* Mniszewski  _et al._ (2020) S. M. Mniszewski, J. Belak, J.-L. Fattebert, C. F. A. Negre, S. R. Slattery, A. A. Adedoyin, R. F. Bird, C. S. Chang, G. Chen, S. Ethier, S. Fogerty, S. Habib, C. Junghans, D. Lebrun-Grandi, J. Mohd-Yusof, S. G. Moore, D. Osei-Kuffuor, S. J. Plimpton, A. Pope, S. T. Reeve, L. Ricketson, A. Scheinberg, A. Y. Sharma, and M. E. Wall, To be submitted to: Special Journal Issue: ECP Co-design and computational motifs (2020).
  * Goedecker and Colombo (1994) S. Goedecker and L. Colombo, Phys. Rev. Lett. 73, 122 (1994).
  * Galli (1996) G. Galli, Cur. Op. Sol. State Mat. Sci. 1, 864 (1996).
  * Tsuchida (2008) E. Tsuchida, J. Phys.: Condens. Matter 20, 294212 (2008).
  * Shimojo  _et al._ (2008) F. Shimojo, R. K. Kalia, A. Nakano, and P. Vashista, Phys. Rev. B 77, 085103 (2008).
  * Cawkwell and Niklasson (2012) M. J. Cawkwell and A. M. N. Niklasson, J. Chem. Phys. 137, 134105 (2012).
  * VandeVondele  _et al._ (2012) J. VandeVondele, U. Borstnik, and J. Hutter, J. Chem. Theory Comput. 8, 3565 (2012).
  * Arita  _et al._ (2014) M. Arita, D. R. Bowler, and T. Miyazaki, J. Chem. Theory Comput. 10, 5419 (2014).
  * Shimojo  _et al._ (2014) F. Shimojo, S. Hattori, R. K. Kalia, M. Kusaneth, W. W. Mou, A. Nakano, K. Nomura, S. Ohmura, P. Rajak, K. Shimamura, and P. Vashista, J. Chem. Phys. 140, 18529 (2014).
  * Osei-Kuffuor  _et al._ (2014) D. Osei-Kuffuor, J. L. Fattebert, and F. Gygi, Phys. Rev. Lett. 112, 046401 (2014).
  * Otsuka  _et al._ (2016) T. Otsuka, M. Taiji, D. R. Bowler, and T. Miyazaki, Japanese Journal of Applied Physics 55, 1102B1 (2016).
  * Niklasson  _et al._ (2016a) A. M. N. Niklasson, S. M. Mnizsewski, C. F. A. Negre, M. J. Cawkwell, P. J. Swart, J. Mohd-Yusof, T. C. Germann, M. E. Wall, N. Bock, E. H. Rubensson, and H. N. Djidjev, J. Chem. Phys. 144, 234101 (2016a).
  * (23) H. N. Djidjev, G. Hahn, S. M. Mniszewski, C. F. Negre, A. M. Niklasson, and V. B. Sardeshmukh, “Graph partitioning methods for fast parallel quantum molecular dynamics,” in  _2016 Proceedings of the Seventh SIAM Workshop on Combinatorial Scientific Computing_, pp. 42–51, http://epubs.siam.org/doi/pdf/10.1137/1.9781611974690.ch5 .
  * Lass  _et al._ (2018) M. Lass, S. Mohr, H. Wiebeler, T. Kühne, and C. Plessl, Proc. Platform for Advanced Scientific Computing (PASC) Conference, ACM (2018).
  * Djidjev  _et al._ (2019a) H. N. Djidjev, G. Hahn, S. M. Mniszewski, C. F. A. Negre, and A. M. N. Niklasson, Algorithms 12 (2019a), 10.3390/a12090187.
  * Lass  _et al._ (2020) M. Lass, R. Schade, T. Kühne, and C. Plessl, Proc. International Conference for High Performance Computing, Networking, Storage and Analysis (SC), IEEE Computer Society , 1127 (2020).
  * Yang (1991) W. Yang, Phys. Rev. Lett. 66, 1438 (1991).
  * Walker and Mezey (1993) P. D. Walker and P. G. Mezey, J. Am. Chem. Soc. 115, 12423 (1993).
  * Yang and Lee (1995) W. T. Yang and T. S. Lee, J. Chem. Phys. 103, 5674 (1995).
  * Abrikosov  _et al._ (1996) I. A. Abrikosov, A. M. N. Niklasson, S. I. Simak, B. Johansson, A. V. Ruban, and H. L. Skriver, Phys. Rev. Lett. 76, 4203 (1996).
  * Kitaura  _et al._ (1999) K. Kitaura, E. Ikeo, T. Nakano, and M. Uebayasi, Chem. Phys. Lett. 313, 701 (1999).
  * Li  _et al._ (2005) S. Li, W. Li, and T. Fang, Journal of the American Chemical Society 127, 7215 (2005), pMID: 15884963, https://doi.org/10.1021/ja0427247 .
  * Ozaki (2006) T. Ozaki, Phys. Rev. B 74, 245101 (2006).
  * Nishimoto  _et al._ (2014) Y. Nishimoto, D. G. Fedorov, and S. Irle, J. Chem. Theory Comput. 10, 4801 (2014).
  * Vuong  _et al._ (2019) V. Q. Vuong, Y. Nishimoto, D. G. Fedorov, B. G. Sumpter, T. A. Niehaus, and S. Irle, Journal of Chemical Theory and Computation 15, 3008 (2019), pMID: 30998360, https://doi.org/10.1021/acs.jctc.9b00108 .
  * Nishimoto and Irle (2021) Y. Nishimoto and S. Irle, “The fmo-dftb method,” in  _Recent Advances of the Fragment Molecular Orbital Method: Enhanced Performance and Applicability_, edited by Y. Mochizuki, S. Tanaka, and K. Fukuzawa (Springer Singapore, Singapore, 2021) pp. 459–485.
  * Stewart (1976) G. W. Stewart, in  _Sparse Matrix Computations_ , edited by J. R. Banch and D. J. Rose (Academic Press, New York, 1976) pp. 113–130.
  * Pissanetzky (1984) S. Pissanetzky,  _Sparse Matrix Technology_ (Academic Press, London, 1984).
  * Brussino and Sonnad (1989) G. Brussino and V. Sonnad, International Journal for Numerical Methods in Engineering 28, 801 (1989).
  * Saad (1996) Y. Saad,  _Iterative methods for Sparse Linear Systems_ (PWS Publishing, Boston, 1996).
  * Schwegler and Challacombe (1996) E. Schwegler and M. Challacombe, J. Chem. Phys. 105, 2726 (1996).
  * Challacombe and Schwegler (1997) M. Challacombe and E. Schwegler, J. Chem. Phys. 106, 5526 (1997).
  * Daniels and Scuseria (1999) A. D. Daniels and G. E. Scuseria, J. Chem. Phys. 110, 1321 (1999).
  * Niklasson (2002) A. M. N. Niklasson, Phys. Rev. B 66, 155115 (2002).
  * Niklasson  _et al._ (2003) A. M. N. Niklasson, C. J. Tymczak, and M. Challacombe, J. Chem. Phys. 118, 8611 (2003).
  * Rubensson and Salek (2005) E. H. Rubensson and P. Salek, J. Comput. Chem 26, 1628 (2005).
  * Jordan and Mazziotti (2005) D. K. Jordan and D. A. Mazziotti, J. Chem. Phys. 122, 084114 (2005).
  * Aradi  _et al._ (2007) B. Aradi, B. Hourahine, and T. Frauenheim, The Journal of Physical Chemistry A 111, 5678 (2007).
  * Rubensson  _et al._ (2008a) E. H. Rubensson, E. Rudberg, and P. Salek, J. Chem. Phys. 128, 074106 (2008a).
  * Niklasson (2011a) A. M. N. Niklasson, “Density matrix methods in linear scaling electronic structure theory,” in  _Linear-Scaling Techniques in Computational Chemistry and Physics: Methods and Applications_, edited by R. Zalesny, M. G. Papadopoulos, P. G. Mezey, and J. Leszczynski (Springer Netherlands, Dordrecht, 2011) pp. 439–473.
  * Buluc and Gilbert (2012) A. Buluc and J. R. Gilbert, SIAM J. Sci. Comput. 34, 170 (2012).
  * Bock and Challacombe (2013) N. Bock and M. Challacombe, SIAM J. Sci. Comput. 35, 72 (2013).
  * Borstnik  _et al._ (2014) U. Borstnik, J. VandeVondele, V. Weber, and J. Hutter, Parallel Computing 40, 47 (2014).
  * Weber  _et al._ (2015) V. Weber, T. Latino, A. Pozdeev, I. Feduova, and A. Curioni, J. Chem. Theory Comput. 11, 3145 (2015).
  * Pinski  _et al._ (2015) P. Pinski, C. Riplinger, E. F. Valeev, and F. Neese, The Journal of Chemical Physics 143, 034108 (2015), https://doi.org/10.1063/1.4926879 .
  * Truflandier  _et al._ (2016) L. A. Truflandier, R. M. Dianzinga, and D. R. Bowler, J. Chem. Phys. 144 (2016), 10.1063/1.4943213.
  * Kruchinina  _et al._ (2016) A. Kruchinina, E. Rudberg, and E. H. Rubensson, Journal of Chemical Theory and Computation 12, 5788 (2016), pMID: 27783507, https://doi.org/10.1021/acs.jctc.6b00626 .
  * Schade  _et al._ (2022a) R. Schade, T. Kenter, H. Elgabarty, M. Lass, O. Schütt, A. Lazzaro, H. Pabst, S. Mohr, J. Hutter, T. D. Kühne, and C. Plessl, Parallel Computing 111, 102920 (2022a).
  * Note (1) Graph-based linear scaling electronic structure theory was recently used by Schade and co-workers, who call it the submatrix method or the non-orthogonal local submatrix method (NOLSM) Schade  _et al._ (2022a), but without referring to the original and equivalent method in Ref. Niklasson  _et al._ (2016a).
  * Niklasson (2017) A. M. N. Niklasson, J. Chem. Phys. 147, 054103 (2017).
  * Niklasson (2020a) A. M. N. Niklasson, J. Chem. Phys. 152, 104103 (2020a).
  * Niklasson (2021) A. M. N. Niklasson, Eur. Phys. J. B 94, 164 (2021).
  * Voter (1997) A. F. Voter, Phys. Rev. Lett. 78, 3908 (1997).
  * Voter and ad T. C. Germann (2002) A. F. Voter and F. M. ad T. C. Germann, Annu. Rev. Mater. Res. 32, 321 (2002).
  * Perez  _et al._ (2009) D. Perez, B. P. Uberuaga, Y. Shim, J. G. Amar, and A. F. Voter (Elsevier, 2009) pp. 79 – 98\.
  * Perez  _et al._ (2016) D. Perez, E. D. Cubuk, A. Waterland, E. Kaxiras, and A. F. Voter, Journal of Chemical Theory and Computation 12, 18 (2016), pMID: 26605853, https://doi.org/10.1021/acs.jctc.5b00916 .
  * Baroni  _et al._ (2001) S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
  * Weber  _et al._ (2004) V. Weber, A. M. N. Niklasson, and M. Challacombe, Phys. Rev. Lett. 92, 193002 (2004).
  * Ochsenfeld  _et al._ (2004) C. Ochsenfeld, J. Kussmann, and F. Koziol, Angewandte Chemie 43, 4485 (2004).
  * Kussmann and Ochsenfeld (2007) J. Kussmann and C. Ochsenfeld, J. Chem. Phys. 127, 204103 (2007).
  * Kussmann  _et al._ (2015) J. Kussmann, A. Luenser, M. Beer, and C. Ochsenfeld, The Journal of Chemical Physics 142, 094101 (2015), https://doi.org/10.1063/1.4908131 .
  * Elstner  _et al._ (1998) M. Elstner, D. Poresag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai, and G. Seifert, Phys. Rev. B 58, 7260 (1998).
  * Finnis  _et al._ (1998) M. W. Finnis, A. T. Paxton, M. Methfessel, and M. van Schilfgarde, Phys. Rev. Lett. 81, 5149 (1998).
  * et al. (2020) B. H. et al., J. Chem. Phys. 152, 124101 (2020).
  * Dewar and Thiel (1977) M. J. S. Dewar and W. Thiel, Theoret.Chim. Acta 46, 89 (1977).
  * Dewar  _et al._ (1985) M. J. S. Dewar, E. G. Zoebisch, E. F. Healy, and J. J. P. Stewart, J. Am. Chem. Soc. 107, 3902 (1985).
  * Stewart (2013) J. J. P. Stewart, J. Mol. Model. 19, 1 (2013).
  * Bannwarth  _et al._ (2018) C. Bannwarth, S. Ehlert, and S. Grimme, J. Chem. Theory Comput. 15, 1652 (2018).
  * Dral  _et al._ (2019) P. O. Dral, X. Wu, and W. Thiel, J. Chem. Theory Comput. 15, 1743 (2019).
  * Malone  _et al._ (2020) W. Malone, B. Nebgen, A. White, Y. Zhang, H. Song, J. A. Bjorgaard, A. E. Sifain, B. Rodriguez-Hernandez, V. M. Freixas, S. Fernandez-Alberti, A. E. Roitberg, T. R. Nelson, and S. Tretiak, Journal of Chemical Theory and Computation 16, 5771 (2020), pMID: 32635739, https://doi.org/10.1021/acs.jctc.0c00248 .
  * Zhou  _et al._ (2020) G. Zhou, B. Nebgen, N. Lubbers, W. Malone, A. M. N. Niklasson, and S. Tretiak, Journal of Chemical Theory and Computation 16, 4951 (2020), pMID: 32609513, https://doi.org/10.1021/acs.jctc.0c00243 .
  * Bannwarth  _et al._ (2020) C. Bannwarth, E. Caldeweyher, S. Ehlert, A. H. ans P. Pracht, J. Seibert, S. Spicher, and S. Grimme, WIREs Comput. Lol. Sci. 11, 1 (2020).
  * Golub and van Loan (1996) G. Golub and C. F. van Loan,  _Matrix Computations_ (Johns Hopkins University Press, Baltimore, 1996).
  * Niklasson (2020b) A. M. N. Niklasson, J. Chem. Theory Comput. 16, 3628 (2020b).
  * Niklasson (2011b) A. M. N. Niklasson, “Density matrix methods in linear scaling electronic structure theory,” in  _Linear-Scaling Techniques in Computational Chemistry and Physics: Methods and Applications_, edited by R. Zalesny, M. G. Papadopoulos, P. G. Mezey, and J. Leszczynski (Springer Netherlands, Dordrecht, 2011) pp. 439–473.
  * Rubensson  _et al._ (2008b) E. H. Rubensson, E. Rudberg, and P. Salek, J. Chem. Phys. 128, 74109 (2008b).
  * Negre  _et al._ (2016) C. F. A. Negre, S. M. Mnizsewski, M. J. Cawkwell, N. Bock, M. E. Wall, and A. M. N. Niklasson, J. Chem. Theory Comput. 12, 3063 (2016).
  * Niklasson  _et al._ (2005) A. M. N. Niklasson, V. Weber, and M. Challacombe, J. Chem. Phys. 123, 44107 (2005).
  * Niklasson  _et al._ (2016b) A. M. N. Niklasson, S. M. Mniszewski, C. F. A. Negre, M. E. Wall, M. J. Cawkwell, and N. Bock, “PROGRESS version 1.0,” (2016b), https://github.com/lanl/qmd-progress.
  * Cawkwell and et al. (2010) M. J. Cawkwell and et al., “LATTE,” (2010), Los Alamos National Laboratory (LA- CC-10004), http://www.github.com/lanl/latte.
  * McWeeny (1956) R. McWeeny, Proc. R. Soc. London Ser. A-Math 235, 496 (1956).
  * Palser and Manolopoulos (1998) A. H. R. Palser and D. E. Manolopoulos, Phys. Rev. B 58, 12704 (1998).
  * Nemeth and Scuseria (2000) K. Nemeth and G. E. Scuseria, J. Chem. Phys. 113, 6035 (2000).
  * Niklasson (2003) A. M. N. Niklasson, Phys. Rev. B 68, 233104 (2003).
  * Rudberg and Rubensson (2011) E. Rudberg and E. H. Rubensson, J. Phys.: Condens. Matter 23, 075502 (2011).
  * Rubensson (2011) E. H. Rubensson, J. Chem. Theory and Comput. 7, 1233 (2011).
  * Suryanarayana (2013) P. Suryanarayana, Chem. Phys. Lett. 555, 291 (2013).
  * Rubensson and Niklasson (2014) E. H. Rubensson and A. M. N. Niklasson, SIAM J. Sci. Comput. 36, 148 (2014).
  * Shang  _et al._ (2021b) H. Shang, W. Liang, Y. Zhang, and J. Yang, Comput. Phys. Commun. 258, 107613 (2021b).
  * Finkelstein  _et al._ (2022) J. Finkelstein, E. H. Rubensson, S. M. Mniszewski, C. F. A. Negre, and A. M. N. Niklasson, Journal of Chemical Theory and Computation 18, 4255 (2022), pMID: 35670603, https://doi.org/10.1021/acs.jctc.2c00274 .
  * Pederson  _et al._ (2022) R. Pederson, J. Kozlowski, R. Song, J. Beall, M. Ganahl, M. Hauru, A. G. M. Lewis, S. B. Mallick, V. Blum, and G. Vidal, “Tensor processing units as quantum chemistry supercomputers,”  (2022).
  * Karypis  _et al._ (1999) G. Karypis, , and V. Kumar, SIAM J. Sci. Comput. 20, 359 (1999).
  * Lab (2013) K. Lab, http://glaros.dtc.umn.edu/gkhome/views/metis (2013).
  * Neidigh  _et al._ (2002) J. W. Neidigh, R. M. Fesinmeyer, and N. H. Andersen, “NMR structure of Trp-Cage miniprotein construct TC5b,” (2002).
  * Djidjev  _et al._ (2019b) H. N. Djidjev, G. Hahn, S. M. Mniszewski, C. F. A. Negre, and A. M. N. Niklasson, Algorithms 12, 187 (2019b).
  * Ghale  _et al._ (2017) P. Ghale, M. P. Kroonblawd, S. Mniszewski, C. F. A. Negre, R. Pavel, S. Pino, V. Sardeshmukh, G. Shi, and G. Hahn, SIAM J. Sci. Comput. 39, C466 (2017).
  * Mermin (1963) N. D. Mermin, Annals of Physics 21, 99 (1963).
  * Mermin (1965) N. D. Mermin, Phys. Rev. B 137, A1441 (1965).
  * Parr and Yang (1989) R. G. Parr and W. Yang,  _Density-functional theory of atoms and molecules_ (Oxford University Press, Oxford, 1989).
  * Niklasson (2008a) A. M. N. Niklasson, J. Chem. Phys. 129, 244107 (2008a).
  * Niklasson (2008b) A. M. N. Niklasson, Phys. Rev. Lett. 100, 123004 (2008b).
  * Souvatzis and Niklasson (2014) P. Souvatzis and A. M. N. Niklasson, J. Chem. Phys. 140, 044117 (2014).
  * Kolafa (2003) J. Kolafa, J. Comput. Chem. 25, 335 (2003).
  * Pulay and Fogarasi (2004) P. Pulay and G. Fogarasi, Chem. Phys. Lett. 386, 272 (2004).
  * Herbert and Head-Gordon (2005) J. Herbert and M. Head-Gordon, Phys. Chem. Chem. Phys. 7, 3269 (2005).
  * Niklasson  _et al._ (2006) A. M. N. Niklasson, C. J. Tymczak, and M. Challacombe, Phys. Rev. Lett. 97, 123001 (2006).
  * Kühne  _et al._ (2007) T. D. Kühne, M. Krack, F. R. Mohamed, and M. Parrinello, Phys. Rev. Lett. 98, 066401 (2007).
  * Yoshida (1990) H. Yoshida, Phys. Lett. A 150, 262 (1990).
  * Grebogi  _et al._ (1990) C. Grebogi, S. M. Hammel, J. A. Yorke, and T. Saur, Phys. Rev. Lett. 65, 1527 (1990).
  * Toxvaerd (1994) S. Toxvaerd, Phys. Rev. E 50, 2271 (1994).
  * Gans and Shalloway (2000) J. Gans and D. Shalloway, Phys. Rev. E 61, 4587 (2000).
  * Engel  _et al._ (2005) R. D. Engel, R. D. Skeel, and M. Drees, J. Comput. Phys. 206, 432 (2005).
  * Bond and Leimkuhler (2007) S. D. Bond and B. J. Leimkuhler,  _Molecular dynamics and the accuracy of numerically computed averages_ (Cambride University Press, United Kingdom, 2007).
  * Toxvaerd  _et al._ (2012) S. Toxvaerd, O. J. Heilmann, and J. C. Dyre, J. Chem. Phys. 136, 224106 (2012).
  * Hammonds and Heyes (2020) K. D. Hammonds and D. M. Heyes, J. Chem. Phys. 152, 024114 (2020).
  * Forest and Ruth (1990) E. Forest and R. Ruth, Physica D 43, 105 (1990).
  * Channel and Scovel (1990) J. P. Channel and C. Scovel, Nonlinearity 3, 231 (1990).
  * McLachlan and Atela (1992) R. McLachlan and P. Atela, Nonlinearity 5, 541 (1992).
  * Leimkuhler and Skeel (1994) B. J. Leimkuhler and R. D. Skeel, J. Comput. Phys. 112, 117 (1994).
  * Finkelstein  _et al._ (2020) J. Finkelstein, C. Cheng, G. Fiorin, B. Seibold, and N. Grønbech-Jensen, The Journal of Chemical Physics 153, 134101 (2020).
  * Car and Parrinello (1985) R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).
  * Bonella  _et al._ (2020) S. Bonella, A. Coretti, R. Vuilleumier, and G. Ciccotti, Phys. Chem. Chem. Phys. 22, 10775 (2020).
  * Coretti  _et al._ (2020) A. Coretti, L. Scalfi, C. Bacon, B. Rotenberg, R. Vuilleumier, G. Ciccotti, M. Salanne, and S. Bonella, The Journal of Chemical Physics 152, 194701 (2020), https://doi.org/10.1063/5.0007192 .
  * Coretti  _et al._ (0) A. Coretti, T. J. Baird, R. Vuilleumier, and S. Bonella, The Journal of Chemical Physics 0, null (0), https://doi.org/10.1063/5.0130117 .
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B:864 (1964).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, 1133 (1965).
  * Dreizler and Gross (1990) R. Dreizler and K. Gross,  _Density-functional theory_ (Springer Verlag, Berlin Heidelberg, 1990).
  * Weinert and Davenport (1992) M. Weinert and J. W. Davenport, Phys. Rev. B 45, R13709 (1992).
  * Wentzcovitch  _et al._ (1992) R. M. Wentzcovitch, J. L. Martins, and P. B. Allen, Phys. Rev. B 45, R11372 (1992).
  * Roothaan (1951) C. C. J. Roothaan, Rev. Mod. Phys. 23, 69 (1951).
  * McWeeny (1959) R. McWeeny, Phys. Rev. 114, 1528 (1959).
  * Hirakawa  _et al._ (2017) T. Hirakawa, T. suzuki, D. R. Bowler, and T. Myazaki, J. Phys.: Condens. Matter 29, 405901 (2017).
  * Note (2) When $n=\rho_{0}[n]$ we have an exact self-consistent solution and then $\rho_{0}[n]=\rho_{\rm min}$.
  * Niklasson  _et al._ (2009) A. M. N. Niklasson, P. Steneteg, A. Odell, N. Bock, M. Challacombe, C. J. Tymczak, E. Holmstrom, G. Zheng, and V. Weber, J. Chem. Phys. 130, 214109 (2009).
  * Steneteg  _et al._ (2010) P. Steneteg, I. A. Abrikosov, V. Weber, and A. M. N. Niklasson, Phys. Rev. B 82, 075110 (2010).
  * Zheng  _et al._ (2011) G. Zheng, A. M. N. Niklasson, and M. Karplus, J. Chem. Phys. 135, 044122 (2011).
  * Niklasson  _et al._ (2015) A. M. N. Niklasson, M. J. Cawkwell, E. H. Rubensson, and E. Rudberg, Phys. Rev. E 92, 063301 (2015).
  * Nishimoto (2017) Y. Nishimoto, J. Chem. Phys. 146, 084101 (2017).
  * McWeeny (1962) R. McWeeny, Phys. Rev. 126, 1028 (1962).
  * Niklasson and Challacombe (2004) A. M. N. Niklasson and M. Challacombe, Phys. Rev. Lett. 92, 193001 (2004).
  * Truflandier  _et al._ (2020) L. A. Truflandier, R. M. Dianzinga, and D. R. Bowler, The Journal of Chemical Physics 153, 164105 (2020), https://doi.org/10.1063/5.0022244 .
  * Weber  _et al._ (2005) V. Weber, A. M. N. Niklasson, and M. Challacombe, J. Chem. Phys. 123, 44106 (2005).
  * Cooley and Tukey (1965) J. W. Cooley and J. W. Tukey, Math. Comput. 19, 297 (1965).
  * Darden  _et al._ (1993) T. Darden, D. York, and L. Pedersen, The Journal of Chemical Physics 98, 10089 (1993), https://doi.org/10.1063/1.464397 .
  * Note (3) Notice, this approach needs to be adapted to work with other representations or the electronic degrees of freedom besides some coarse grained charge density, e.g. of partial atomic charges.
  * Broyden (1965) C. G. Broyden, Math. Comput. 19, 577 (1965).
  * Anderson (1965) D. G. Anderson, J. Assoc. Comput. Mach. 12, 547 (1965).
  * Pulay (1980) P. Pulay, Chem. Phys. Let. 73, 393 (1980).
  * Porezag  _et al._ (1995) D. Porezag, T. Frauenheim, T. Köhler, G. Seifert, and R. Kaschner, Phys. Rev. B 51, 12947 (1995).
  * Seifert  _et al._ (1996) G. Seifert, D. Porezag, and T. Frauenheim, International Journal of Quantum Chemistry 58, 185 (1996).
  * Frauenheim  _et al._ (2000) T. Frauenheim, G. Seifert, M. Elstner, Z. Hajnal, G. Jungnickel, D. Poresag, S. Suhai, and R. Scholz, Phys. Stat. sol. 217, 41 (2000).
  * Gaus  _et al._ (2011) M. Gaus, Q. Cui, and M. Elstner, J, Chem. Theory Comput. 7, 931 (2011).
  * Krishnapryian  _et al._ (2017) A. Krishnapryian, P. Yang, A. M. N. Niklasson, and M. J. Cawkwell, J. Chem. Theory Comput. 13, 6191 (2017).
  * Pulay  _et al._ (1979) P. Pulay, , Fogarasi, F. Pang, and J. E. Boggs, J. Am. Chem. Soc. 101, 2550 (1979).
  * Berendsen  _et al._ (1995) H. Berendsen, D. Vanderspoel, and R. Vandrunen, Computer Phys. Comm. 91, 43 (1995).
  * (166) “CHARMM-GUI,” https://www.charmm-gui.org/, accessed: 2022-11-27.
  * (167) R. P. D. Bank, “1L2Y,” https://www.rcsb.org/structure/1L2Y, accessed: 2022-11-27.
  * (168) “AmberTools22,” https://ambermd.org/AmberTools.php, accessed: 2022-11-27.
  * Schade  _et al._ (2022b) R. Schade, T. Kenter, H. Elgabarty, M. Lass, T. D. Kühne, and C. Plessl, “Breaking the exascale barrier for the electronic structure problem in ab-initio molecular dynamics,”  (2022b).
  * Dral  _et al._ (2015) P. O. Dral, O. A. von Lilienfeld, and W. Thiel, Journal of Chemical Theory and Computation 11, 2120 (2015), pMID: 26146493, https://doi.org/10.1021/acs.jctc.5b00141 .
  * Ramakrishnan  _et al._ (2015) R. Ramakrishnan, P. O. Dral, M. Rupp, and O. A. von Lilienfeld, Journal of Chemical Theory and Computation 11, 2087 (2015), pMID: 26574412, https://doi.org/10.1021/acs.jctc.5b00099 .
  * Li  _et al._ (2018) H. Li, C. Collins, M. Tanha, G. J. Gordon, and D. J. Yaron, Journal of Chemical Theory and Computation 14, 5764 (2018), pMID: 30351008, https://doi.org/10.1021/acs.jctc.8b00873 .
  * Kranz  _et al._ (2018) J. J. Kranz, M. Kubillus, R. Ramakrishnan, O. A. von Lilienfeld, and M. Elstner, Journal of Chemical Theory and Computation 14, 2341 (2018), pMID: 29579387, https://doi.org/10.1021/acs.jctc.7b00933 .
  * Goldman  _et al._ (2018) N. Goldman, B. Aradi, R. K. Lindsey, and L. E. Fried, Journal of Chemical Theory and Computation 14, 2652 (2018), pMID: 29614217, https://doi.org/10.1021/acs.jctc.8b00165 .
  * Zheng  _et al._ (2021) P. Zheng, R. Zubatyuk, W. Wu, O. Isayev, and P. O. Dral, Nature Communications 12, 7022 (2021).
  * Zhou  _et al._ (2022) G. Zhou, N. Lubbers, K. Barros, S. Tretiak, and B. Nebgen, Proc. Nat. Ac. Sci. 119, 2120333119 (2022).
  * Hu  _et al._ (2022) F. Hu, F. He, and D. J. Yaron, “https://arxiv.org/abs/2210.11682: semiempirical hamiltonians learned from data can have accuracy comparable to density functional theory,”  (2022).


# Non-Debye relaxations: The characteristic exponent in the excess wings model

K. Górska<EMAIL_ADDRESS>A. Horzela<EMAIL_ADDRESS>Institute of Nuclear Physics, Polish Academy of Sciences,
ul. Radzikowskiego 152, PL-31342 Kraków, Poland  T. K. Pogány  Faculty of
Maritime Studies, University of Rijeka, 51000 Rijeka, Croatia
Institute of Applied Mathematics, Óbuda University, H-1034 Budapest, Hungary
<EMAIL_ADDRESS>
###### Abstract

The characteristic (Laplace or Lévy) exponents uniquely characterize
infinitely divisible probability distributions. Although of purely
mathematical origin they appear to be uniquely associated with the memory
functions present in evolution equations which govern the course of such
physical phenomena like non-Debye relaxations or anomalous diffusion. Commonly
accepted procedure to mimic memory effects is to make basic equations time
smeared, i.e., nonlocal in time. This is modeled either through the
convolution of memory functions with those describing relaxation/diffusion or,
alternatively, through the time smearing of time derivatives. Intuitive
expectations say that such introduced time smearings should be physically
equivalent. This leads to the conclusion that both kinds of so far introduced
memory functions form a “twin” structure familiar to mathematicians for a long
time and known as the Sonine pair. As an illustration of the proposed scheme
we consider the excess wings model of non-Debye relaxations, determine its
evolution equations and discuss properties of the solutions.

excess wings relaxation, time smeared evolution equations, memory functions,
characteristic exponents

###### pacs:

02.50.Ey, 02.30.Uu, 02.30.Gp

##  I Introduction

Typical example of dielectric relaxation is provided by a dipolar system which
approaches the equilibrium being earlier driven out of it by a step or
alternating external electric field. The phenomenon is usually described in
terms of the relaxation function $n(t)$ which counts dipoles surviving
depolarization during the time $(0,t)\subset(0,\infty)$ and, if normalized,
evolves form $n(0+)=1$ to $n(\infty)=0$. The function $n(t)$ comes out as the
solution of macroscopic differential equation

$\dot{n}(t)=-r(t,\tau)n(t).$ (1)

The non-negative quantity $r(t,\tau)$ is the transition rate of the system and
besides of the time depends on properties characterizing the medium among
which a material constant called the relaxation, or characteristic, time
$\tau$ is the most important. Solution to Eq. (1) is easily got as
$n(t)=\exp[-\int_{0}^{\infty}r(\xi,\tau)\operatorname{d\\!}\xi]$ but it
remains of very limited physical utility because the knowledge of $r(t,\tau)$,
especially for short and long times $t$, is far insufficient except of the
Debye case for which $r(t,\tau)={\tau}^{-1}={\rm const}$. Data provided by the
broadband dielectric spectroscopy (encoded in the so–called spectral
functions) extrapolated to the full frequency range and next transformed to
the time domain, do not help very much - using them to calculate the ratio
$\dot{n}(t)/n(t)$ usually leads to cumbersome formulae Gloeckle ; RHilfer02 ;
RHilfer02a ; KGorska18 , in addition singular at the origin, which makes their
experimental verification rather impossible for the time laps close to the
origin. These difficulties have prompted efforts to look for mesoscopic
description of relaxation phenomena based on dynamical rules being non-local
in time and leading to the evolution equations which from the very beginning
take into account the memory effects. The simplest way to mimic the memory is
to introduce the time smearing which may proceed two-fold: either one smears
the left hand side of Eq. (1), i.e., the time derivative in
$\dot{n}(t)=-r(t,\tau)n(t)$ or rewrites Eq. (1) in the integral form
${n}(t)=1-\int_{0}^{t}r(\xi,\tau)n(\xi)\operatorname{d\\!}\xi$ and uses the
smearing $r(\xi,\tau)\to r(t-\xi,\tau)$. Keeping the relaxation time $\tau$
explicitly separated out from the other material depending parameter’s vector
${p}$, say, this leads to

$\int_{0}^{t}k(t-\xi)~{}\dot{n}(\xi)\,\operatorname{d\\!}\xi=-B(\tau;{p})~{}n(t)\,,$
(2)

where $B(\tau,{p})$ denotes the universal, time independent, transition rate
and $k(t)$ stands for the memory kernel responsible for smearing the time
derivative. In turn, for $r(t-\xi,\tau)=B(\tau;{p})M(t-\xi)$, we get

$n(t)=1-B(\tau,{p})\int_{0}^{t}M(t-\xi)~{}n(\xi)\,\operatorname{d\\!}\xi,$ (3)

with $M(t-\xi)$ being another memory kernel, a priori not connected to
$k(t-\xi)$. Mathematically Eqs. (2) and (3) are both the Volterra type
equations GGripenberg  which utility goes beyond more popular fractional
differential equations introduced in the framework of fractional calculus
approach to the relaxation phenomena KGorska20 . If we require physically
justified equivalence of Eqs. (2) and (3) then the memory kernels become
mutually related and form a coupled (Sonine) pair which appearance and
properties we shall discuss a bit later.

The integro-differential equation (2) is mathematically very well understood
ANKochubei11 . In fact it is the equation which for special choices of
$k(t-\xi)$ reduces to equations with fractional derivatives (more precisely,
various types of them) for a long time proposed to investigate the relaxation
phenomena. Simultaneously, both Eqs. (2) and (3) have the form of kinetic
equations which are the starting point to describe relaxation in the
subordination framework AStanislavsky15 , AStanislavsky17 ,
(AStanislavsky19_1, , Chs. 4.1, 4.3) developed as a general scheme within the
stochastic processes approach to the relaxation phenomena and anomalous
diffusion.

The cornerstone of the stochastic processes based approach to relaxation (as
well as to anomalous diffusion if one adopts a suitably reinterpreted
language) is the assumption that the transition rate $r(t,\tau)$ introduced in
Eq. (1) takes on the meaning of a non-negative stochastic quantity
parametrized by the randomized characteristic time $\tau$. The latter does not
denote any longer one among material properties of the relaxing medium and
becomes physically meaningful variable which shape of postulated randomization
strongly influences, or even determines, modeling the relaxation. The choice
of stochastic processes proposed to investigate relaxation phenomena is domi-
nated by choosing those which are non-negative, non-decreasing and have
distributions which are infinitely divisible. The last means that relavant
distributions functions $f(\widetilde{\beta})$ are representable as
$N\to\infty$ limit of distributions obeyed by random variables
$\widetilde{\beta}^{(N)}=\frac{1}{N}\sum_{i=1}^{N}\lambda_{i}$ where all
$\lambda_{i}$ are independent identically distributed random variables L1 .
Randomization of the characteristic time $\tau$ (in the Debye systems assumed
to be the same and fixed for all dipoles forming the system) means that we are
going to change description of the system - instead of looking for
deterministic evolution in the time $t$ measured by a laboratory clock we
search for stochastic evolution in terms of the “internal” time $\tau(t)$
whose dependence on $t$ is hidden in some probability distribution
$f(\tau,t)$.

Any non-negative stochastic process whose distribution is infinitely
divisible, herewith denoted as $U(\tau)$, satisfies the relation

$\left\langle\exp{(-sU(\tau))}\right\rangle=\exp{(-\tau\widehat{\Psi}(s))},$
(4)

where $\widehat{\Psi}(s)$ bears the name of characteristic (either Laplace or
Lévy) exponent and is uniquely given by the Lévy–Khintchine formula
(Schilling, , Eq. (1.3))

$\widehat{\Psi}(s)=\lambda
s+\int_{0}^{\infty}(1-\operatorname{e}^{-sx})\mu(\operatorname{d\\!}x)$ (5)

where $\mu(\operatorname{d\\!}x)$ (subject to some additional conditions) is
called the Lévy measure while $\lambda$ is named the drift parameter. For
$s>0$ the relation (5) places all functions $\widehat{\Psi}(s)$ in the class
of Bernstein functions (BFs), i.e. non-negative functions on $\mathbb{R}_{+}$,
differentiable infinitely many times and satisfying for $s>0$ and
$n\in\mathbb{N}_{0}$ the conditions $(-1)^{n}f^{(n+1)}(s)\geq 0$ everywhere in
their domain L2 . We remark that the BFs are close relatives to the completely
monotone functions (CMFs), also being non-negative on $\mathbb{R}_{+}$,
differentiable there infinitely many times and satisfying
$(-1)^{n}f^{(n)}(s)\geq 0,s>0,n\in\mathbb{N}_{0}$. To make these notions more
intuitive one may understand BFs as “maximally regularly” increasing positive
functions while CMFs as “maximally regularly” decreasing ones L3 . The deep
mutual relation between infinitely divisible distributions, BFs and CMFs is
encoded as follows: for $h:(0,\infty)\mapsto(0,\infty)$ the following
statements are equivalent: (i) $h$ is CMF and it is infinitely divisible with
$h(0+)\leq 1$; and (ii) $h=\exp(-{\widehat{\Psi}})$ with ${\widehat{\Psi}}$
being BF (RSchilling10, , p. 52, Lemma 5.8). Coming back to the relaxation
phenomena we remind that the relaxation function $n(t)$ (which provides us
with the information on the number of relaxation centers which did not decay
during the time $(0,t)$) if calculated from the spectroscopic data appears to
be CMF for a vast majority of commonly used phenomenological models Hanyga1 ;
EDeOliveira11 ; RGarrappa11 ; Tomovski ; RGarrappa16 ; KGorska20a ; KGorska21
. This fact merged with the just mentioned theorem strongly suggests that
characteristic exponents are inextricably linked with investigation of the
relaxation processes. Research which sheds light on this problem is the
leitmotif of our paper.

We present and discuss a number of arguments which clarify the role played by
characteristic exponents in description of the relaxation phenomena, in
particular provide the reader their interpretation as memory functions. The
methods which we advocate are general and, as recently demonstrated in
KGorska21a ; SW21 , applicable to various phenomenological models of
relaxation. In what follows we focus our attention on the excess wings model
RHilfer02b ; RHilfer02c  which goes beyond the Jonscher universal relaxation
law (URL) AJonscher92  and is less popular among experimentalists if compared
with models of the Havriliak-Negami family. General considerations of Sec. II
show how the characteristic exponent enters the spectral, relaxation, and
memory functions. Also we explain the physical interpretation of
${\widehat{\Psi}}(s)$ and demonstrate that required equivalence of Eqs. (2)
and (3) inevitably leads to the concept of the Sonine pair which
non–negligible role in theoretical studies of relaxation, viscoelasticity and
anomalous diffusion was recently noticed, analysed and developed AHanyga20 ;
AGiusti20 ; AGiusti20a . Starting from Sec. III we investigate the excess
wings model. Using its spectral function we recover suitable characteristic
exponent whose knowledge enables us to find the appropriate relaxation
function. In Sec. IV we use the characteristic exponent to introduce two
coupled memory kernel functions which form the Sonine pair. Thus we arrive at
a pair of evolution equations which involve either the smearing of the
relaxation function or its time derivative and should give the same excess
wings relaxation function. Both equations are solved in Sec. V where also
requirements demanded from their solutions are checked. The paper is concluded
in Sec. VI.

##  II Characteristic exponents as constitutive elements of the relaxation
theory

As sygnalized in the Introduction the first step in the construction of
stochastic approach to relaxation phenomena, see AStanislavsky15  and
AStanislavsky17 ; AStanislavsky19_1  for recent exhaustive reviews, is to
assume that they are underpinned by randomization of the characteristic time
$\tau$ and that the stochastic processes $U(\tau)$ emerging from such a
randomization have non–negative infinitely divisible distributions. In the
majority of physically meaningful applications these distributions are
realized as heavy tailed $\alpha$-stable Lévy ones related to various variants
of random walks. The second step, essential for making the method effective in
modeling physical applications, is to use the subordination formalism Bochner
within which the parent process, usually the Debye law dependent on
operational time $\tau$, is subordinated by a directing process $\tau(t)$
which links $\tau$ and $t$ in a random relation encoded in a probability
density (pdf) $f(\tau,t)$. Intuitively, employing the subordination scheme
means to replace a process described in terms of the laboratory clock measured
time $t$ by a composed random process governed by an irregular non-decreasing
flow of randomized time $\tau$ given by a stochastic process
$t\to\tau=\tau(t)$. Physically it is expected that properties of $\tau(t)$ may
shed light on the internal structure of the system or provide us with some
hints how its macroscopic behaviour is influenced by many-body effects.

According to Eq. (4) the probability theory introduces the characteristic
exponent $\widehat{\Psi}(s)$ in terms of the mean value of the exponentiated
non-negative stochastic process $\exp{(-sU(\tau))}$. Suppose that $f(\tau,t)$
(which says how to find the system in the operational time $\tau$ if it is in
the laboratory time $t$) is also the infinitely divisible pdf of $U(\tau)$.
Then,

$\displaystyle\begin{split}\mathcal{L}[f(\tau,t);s]&=\int_{0}^{\infty}\operatorname{e}^{-st}f(\tau,t)\operatorname{d\\!}t\\\
&=\left\langle\operatorname{e}^{-sU(\tau)}\right\rangle=\operatorname{e}^{-\tau{\widehat{\Psi}}(s)},\quad\tau>0.\end{split}$
(6)

Assumption that the process $U(\tau)$ results from $t\to\tau:U(\tau)\leq t$
with the pdf $f(\tau,t)$ opens the possibility to ask for the “inverse”
process $S(t)={\rm inf}\\{\tau:U(\tau)>t\\}$ and its pdf $g(t,\tau)$. The
latter may be calculated from the cumulant distribution functions of $U(\tau)$
and $S(t)$ (see e.g. AStanislavsky15 ; ChechkinSokolov21 )

$g(t,\tau)=-\frac{\partial}{\partial\tau}\int_{0}^{\,t}f(\tau,\xi)\operatorname{d\\!}\xi.$
(7)

Taking the Laplace transform of Eq. (7) and using Eq. (6) leads to L4

$\widehat{g}(s,\tau)=\frac{\widehat{\Psi}(s)}{s}\operatorname{e}^{-\tau\widehat{\Psi}(s)}.$
(8)

Alternatively, any non-Debye relaxation process may be seen as summing up
effects of multichannel exponential decays with each channel characterized by
some randomly distributed relaxation time $\theta$. Under this assumption the
relaxation function $n(t)$ counting the fraction of objects which have
survived the decay in the laboratory time interval $(0,t)$ boils down to the
weighted average of exponential decays

$n(t)=\int_{0}^{\infty}e^{-t\theta}\mu(\operatorname{d\\!}\theta),$ (9)

where $\mu(\operatorname{d\\!}\theta)$ denotes the probability with which the
random relaxation time $\theta$ occurs. In the framework of the subordination
approach the same quantity $n(t)$ comes from weighted average of the Debye law
expressed in the operational time $\tau$ and the pdf $g(t,\tau)$. Thus Eq. (9)
may be rewritten as the integral decomposition Fogedby94

$n(t)=\int_{0}^{\infty}e^{-B({p})\,\tau}g(t,\tau)\operatorname{d\\!}\tau.$
(10)

Using Eq. (10) enables us to calculate the response (called also spectral)
function defined in the frequency domain as
$\widehat{\phi}(\operatorname{i}\\!\omega)={\cal
L}^{-1}[-\dot{n}(t);\operatorname{i}\\!\omega]$. Because of Eqns. (8) and (10)
it is uniquely expressed in terms of the characteristic exponent

$\widehat{\phi}(\operatorname{i}\\!\omega)=\frac{1}{1+\widehat{\Psi}(\operatorname{i}\\!\omega)/B({p})}.$
(11)

Here we point out that Eq. (11) explicitly determines the relation between
purely phenomenological object which is the spectral function $\widehat{\phi}$
obtained as a fit to experimental data and $\widehat{\Psi}$, a mathematical
quantity one to one related to the stochastic process being assumed to
underlie physical phenomenon under consideration but of origin rather loosely
supported by specific physical properties of the system. To look for physical
justification of so far presented construction notice that the relation

$\widehat{n}(s)=\frac{1-\widehat{\phi}(s)}{s}.$ (12)

and Eq. (11) implies

$\widehat{n}(s)=\frac{s^{-1}}{1+B(p)/\widehat{\Psi}(s)}.$ (13)

As recalled in the Introduction the time evolution equations involving memory
effects may be obtained by modeling memory effects through the time smearing,
either of $\dot{n}(t)$ like it has taken place in Eq. (2) or of
$r(t,\tau)n(t)$ like has been done in Eq. (3). Doing that we arrive at linear
integro-differential equations which without difficulties may be solved in the
Laplace domain. The relaxation function which solves Eq. (2) in the Laplace
domain reads

$\widehat{n}_{k}(s)=\frac{s^{-1}}{1+B({p})/[s~{}\widehat{k}(s)]}.$ (14)

while for Eq. (3) we get

$\widehat{n}_{M}(s)=\frac{s^{-1}}{1+B({p})~{}\widehat{M}(s)},$ (15)

where $\widehat{k}(s)=\mathcal{L}[k(t);s]$ and
$\widehat{M}(s)=\mathcal{L}[M(t);s]$. Physical equivalence of the above
approaches requires that Eqs. (14) and (15) describe the same situation, i.e.,
the memory effects influencing the behaviour of $\dot{n}(t)$ and
$r(t,\tau)n(t)$ should yield the same results for observed properties of
$n(t)$. The equality of Eqs. (14) and (15), i.e.,
$\widehat{n}_{k}(s)=\,\widehat{n}_{M}\,=\widehat{n}(s)$, if compared with Eq.
(13), gives

$\widehat{M}(s)=[s~{}\widehat{k}(s)]^{-1}=[\widehat{\Psi}(s)]^{-1},$ (16)

which merges the deterministic, i.e. evolution equations stemmed, description
of the relaxation with its stochastic roots. Consequently, the stochastic
nature of relaxations puts rigid restrictions on properties of admissible
memory functions, in fact deeply reaching for their analyticity structure.
This is because the memory functions form not only the Sonine pair written
down in the Laplace domain as $\widehat{M}(s)\widehat{k}(s)=s^{-1}$ but being
directly related to the characteristic exponents and Bernstein functions (both
living on the positive semiaxis) may be consistently extended to the complex
domain where they fall into special classes of analytic functions, namely the
Stieltjes and Nevanlinna-Pick functions GGripenberg ; Berg .

##  III The excess wings model

The spectral function which corresponds to the simplest version of the excess
wings model is

$\widehat{\phi}_{\alpha}(\operatorname{i}\\!\omega)=\frac{1+(\operatorname{i}\\!\omega\tau_{2})^{\alpha}}{1+\operatorname{i}\\!\omega\tau_{1}+(\operatorname{i}\\!\omega\tau_{2})^{\alpha}},\quad\alpha\in(0,1).$
(17)

It depends on two characteristic times $\tau_{1}>0$ and $\tau_{2}>0$ and so
does not fit to the Jonscher’s URL RGarrappa16 ; AJonscher92 ; SW16  involving
only a single characteristic time $\tau$. Despite this reservation the excess
wings model appears useful in analysis of experimental data as it successfully
describes the relaxation phenomena in the high frequency regime when the
frequency of applied electric field is of the order $10^{5}-10^{10}$ Hz
RHilfer02b ; RHilfer02c ; PDixon90 ; PDixon90a ; RBrand00 . For $\alpha=1$ the
spectral function $\widehat{\phi}_{\alpha}(\operatorname{i}\\!\omega)$ is
proportional to the Debye spectral function
$\widehat{\phi}_{D}(\operatorname{i}\\!\omega)=[1+(\tau_{1}+\tau_{2})\operatorname{i}\\!\omega]^{-1}$
with the characteristic time $\tau_{1}+\tau_{2}$, i.e.,
$\widehat{\phi}_{1}(\operatorname{i}\\!\omega)=\widehat{\phi}_{D}(\operatorname{i}\\!\omega)+\operatorname{i}\\!\omega\tau_{2}\widehat{\phi}_{D}(\operatorname{i}\\!\omega)$.

Comparing the spectral function (17) with (11) we find that the characteristic
exponent $\widehat{\Psi}$ formally reads

$\widehat{\Psi}(\operatorname{i}\\!\omega)=\frac{\operatorname{i}\\!\omega}{\tau_{2}^{-\alpha}+(\operatorname{i}\\!\omega)^{\alpha}},\quad\alpha\in(0,1)$
(18)

if we set
$B(\tau,{p})=B(\tau_{1},\tau_{2},\alpha)=\tau_{2}^{\alpha}/\tau_{1}=\widetilde{\tau}$.
But some doubt arises: is the construction described in Sec. II legitimate if
we have two characteristic times - which of them, and how, is randomized? To
find out properties of $\widehat{\Psi}$ without referring to the Lévy-
Kchintchine formula consider the function
$\widehat{\Psi}(\operatorname{i}\\!\omega)$ given by Eq. (18) as the function
$\widehat{\Psi}(z)$ of a complex variable $z\in\mathbb{C}$. As shown in
(EDeOliveira11, , Eq. (2.22) et seq.) this function satisfies all conditions
of (GGripenberg, , Theorem 2.6) or (EDeOliveira11, , Theorem). It leads to the
crucially important result - namely enables us to represent in an unique way
$\widehat{\Psi}(z)$ as the Laplace transform of a non-negative function.
Furthermore, restricting the argument $z$ of $\widehat{\Psi}(z)$ to the
positive semiaxis, i.e., $z=s>0$, we can identify
$\widehat{\Psi}(s)|_{s\in\mathbb{R}_{+}}$ as a BF and make use of a plethora
of results concerning CMFs and BFs. In the first step notice that for $s>0$
the function $\widehat{\Psi}(s)$ is non-negative while its first derivative

$\frac{\operatorname{d\\!}\widehat{\Psi}(s)}{\operatorname{d\\!}s}=\frac{1-\alpha}{\tau_{2}^{-\alpha}+s^{\alpha}}+\frac{\alpha\tau_{2}^{-\alpha}}{(\tau_{2}^{-\alpha}+s^{\alpha})^{2}}$

is CMF.

$f(s)$ |  CMF |  BF
---|---|---
$(s+b)^{\,\mu}$, $b\geq 0$ | $\mu\leq 0$ | $\mu\in(0,1)$
$(s^{\nu}+b)^{\,\mu}$, $b>0$ |  $\mu\leq 0$ and $\nu\in(0,1]$ |  $\mu\in(0,1)$ and $\nu\in(0,1)$
Table 1:  Examples of CMF and BF.

Indeed, from Tab. 1 we see that for non-negative $\tau_{2}$ and
$\alpha\in(0,1)$ this expression is a convex sum of CMFs and hence it is CMF
as well. Thus, the characteristic exponent $\widehat{\Psi}(s)$ itself is BF.
This is the result which we do need and which for the case under consideration
is by no means obvious from the stochastic point of view since we lack the
information concerning the infinite divisibility of underlying stochastic
process. Needed result, which obviously confirms infinite divisibility, is
obtained from the completely different sources, namely from the phenomenology
merged with mathematical analysis. We would also like to remark that within
the stochastic approach we deal with functions of real variables exemplified
by those being CMFs and BFs. Starting from the spectral function treated as a
complex function of the complex variable we avoid the path marked out by
principles of the stochastic approach. Equipped with tools of the complex
analysis we can leave aside the probability rooted description of relaxation
phenomena and may understand much better results not once or twice hidden
behind paradigms of the real functions approach KGorska21 .

The relaxation function $n(t)$ Eq. (10) with substituted Eq. (8) reads

$\displaystyle n(t)$
$\displaystyle=\int_{0}^{\infty}\\!\\!\operatorname{e}^{-\xi\widetilde{\tau}}\mathcal{L}^{-1}\\!\left[\frac{1}{\tau_{2}^{-\alpha}+s^{\alpha}}\exp\Big{(}-\frac{\xi
s}{\tau_{2}^{-\alpha}+s^{\alpha}}\Big{)};t\right]\operatorname{d\\!}\xi$
$\displaystyle=\mathcal{L}^{-1}[(s+s^{\alpha}\widetilde{\tau}+\tau_{1}^{-1})^{-1};t],$
(19)

where $\widetilde{\tau}=\tau_{2}^{\alpha}/\tau_{1}$ was introduced a few lines
above. To get Eq. (III) we changed the order of integration over
$\xi\in[0,\infty)$ in the inverse Laplace transform which reduced the integral
over $\xi$ to the elementary one. The inverse Laplace transform in the lower
line of Eq. (III) can be calculated by virtue of the formula (TSandev19, , p.
10, Eq. (1.38))

$\mathcal{L}^{-1}\Big{[}\dfrac{s^{-\beta}}{1+\lambda_{1}s^{-\alpha_{1}}+\lambda_{2}s^{-\alpha_{2}}};t\Big{]}\\\
=t^{\,\beta-1}\,E_{(\alpha_{1},\alpha_{2}),\,\beta}\big{(}-\lambda_{1}t^{\alpha_{1}},-\lambda_{2}t^{\alpha_{2}}\big{)},$

from which we get

$n(t)=E_{(1,1-\alpha),1}(-t/\tau_{1},-\widetilde{\tau}t^{1-\alpha}).$ (20)

As shown in RRNigmatulin16  this function, known as the binomial
(multivariable) Mittag-Leffler function, is defined by the double power series

$E_{(\alpha_{1},\alpha_{2}),\,\beta}(x,y)=\sum_{k\geq
0}\sum_{\underset{l_{1}+l_{2}=k}{l_{1},l_{2}\geq
0}}\frac{k!}{l_{1}!\,l_{2}!}\,\frac{x^{l_{1}}y^{l_{2}}}{\Gamma(\beta+\alpha_{1}l_{1}+\alpha_{2}l_{2})},$
(21)

,$x,y\in\mathbb{R}$, and is non-negative for $\lambda_{1},\lambda_{2}\geq 0$,
$\beta\in(0,1)$, and $\alpha_{1},\alpha_{2}\geq\beta-1$. Thus, $n(t)$ given by
Eq. (20) is non-negative for $\alpha\in(0,1)$. Notice that in Eq. (21) the
infinite sum over $k$ is followed by sums over $l_{1}$ and $l_{2}$ constrained
by $l_{1}+l_{2}=k$. As a consequence the double sum in $l_{1}$ and $l_{2}$ can
be represented two-fold: (a) $l_{1}=0,1,\ldots,k$ and $l_{2}=k-l_{1}$ or (b)
$l_{2}=0,1,\ldots k$ and $l_{1}=k-l_{2}$. Without loss of generality we
consider the case (a). In such a case Eq. (20) becomes

$\displaystyle n(t)$ $\displaystyle=\sum_{k\geq
0}\sum_{l_{1}=0}^{k}\binom{k}{l_{1}}\frac{(-t/\tau_{1})^{k}(\tau_{2}/t)^{\alpha(k-l_{1})}}{\Gamma[1+k-\alpha(k-l_{1})]}$
(22) $\displaystyle=\sum_{k\geq
0}(-\widetilde{\tau}t^{1-\alpha})^{k}\sum_{l_{1}=0}^{k}\binom{k}{l_{1}}\frac{(t/\tau_{2})^{\alpha
l_{1}}}{\Gamma[1+(1-\alpha)k+\alpha l_{1}]}.$

Using the definition of Mittag-Leffler polynomials (42) we get

$n(t)=\sum_{k\geq
0}(-\widetilde{\tau}t^{1-\alpha})^{k}E_{\alpha,(1-\alpha)k+1}^{-k}[-(t/\tau_{2})^{\alpha}].$
(23)

The same expression as in Eq. (22) will be obtained if we change $\sum_{k\geq
0}\sum_{l_{1}=0}^{k}$ into $\sum_{l_{1}\geq 0}\sum_{k\geq l_{1}}$. Making this
change and setting $r=k-l_{1}$ we transform the series and the sum sitting
inside Eq. (22) into two independent series

$n(t)=\sum_{l_{1}\geq 0}\sum_{r\geq
0}\frac{(l_{1}+r)!}{r!l_{1}!}\frac{(-\widetilde{\tau}t^{1-\alpha})^{r}(-t/\tau_{1})^{l_{1}}}{\Gamma[1+(1-\alpha)r+l_{1}]}.$
(24)

Treating once the series over $l_{1}$ and another time the series over $r$ as
the definition (37) of the three parameter Mittag–Leffler (or Prabhakar)
function L5  we express Eq. (24) in two equivalent forms, namely

$\displaystyle n(t)$ $\displaystyle=\sum_{l_{1}\geq
0}(-t/\tau_{1})^{l_{1}}E^{1+l_{1}}_{1-\alpha,1+l_{1}}(-\widetilde{\tau}t^{1-\alpha})$
(25) $\displaystyle=\sum_{r\geq
0}(-\widetilde{\tau}t^{1-\alpha})^{r}E^{1+r}_{1,(1-\alpha)r+1}(-t/\tau_{1})\,.$
(26)

Calculations made for (a) can be repeated for (b) with $l_{2}$ written instead
of $l_{1}$; thus, $r=k-l_{2}$. The formulae (25) and (26) reproduce the
relations (RGarrappa16, , Eqs. (3.73), (3.71)) up to the multiplicative
constant $1/\tau_{1}$. Moreover, we conclude that $n(0+)=1$.

##  IV Evolution equation

###  IV.1 Smearing of $r(t,\tau)n(t)$

First we check what equation is satisfied by $n(t)$. For that purpose we take
Eq. (23). In (AMMathai08, , Theorem 2.3.1. on p. 93), i.e.
$zE_{\mu,\nu}^{\gamma}(z)=E_{\mu,\nu-\mu}^{\gamma}(z)-E_{\mu,\nu-\mu}^{\gamma-1}(z)$,
we set $\mu=\alpha$, $\nu=(1-\alpha)k+1$, and $\gamma=-k$. That allows us to
rewrite Eq. (23) in the form

$\displaystyle
n(t)=-\tau_{1}\sum_{k=0}^{\infty}(-\widetilde{\tau}\,)^{\,k+1}t^{(1-\alpha)(k+1)-1}$
$\displaystyle\times\left\\{E^{-(1+k)}_{\alpha,(1-\alpha)(1+k)}[-(t/\tau_{2})^{\alpha}]-E^{-k}_{\alpha,(1-\alpha)(1+k)}[-(t/\tau_{2})^{\alpha}]\right\\}$
$\displaystyle=-\tau_{1}\left\\{\sum_{r=1}^{\infty}(-\widetilde{\tau}\,)^{\,r}t^{(1-\alpha)r-1}E^{\,-r}_{\alpha,(1-\alpha)r}[-(t/\tau_{2})^{\alpha}]\right.$
$\displaystyle\left.+\widetilde{\tau}\sum_{k=0}^{\infty}(-\widetilde{\tau}\,)^{k}t^{(1-\alpha)k-\alpha}E^{-k}_{\alpha,(1-\alpha)k+1-\alpha}[-(t/\tau_{2})^{\alpha}]\right\\}$
(27)

where we change the summation index in the first series by setting $k+1=r$.
From Eqs. (39) and (40) it comes out that Eq. (27) can be expressed as

$-\frac{1}{\tau_{1}}\,n(t)=\Big{[}\frac{\operatorname{d\\!}}{\operatorname{d\\!}t}+\widetilde{\tau}D_{t}^{\alpha}\Big{]}n(t),\quad\alpha\in(0,1),$
(28)

where
$({D_{x}^{\alpha}}f)(x)=(\frac{\operatorname{d\\!}}{\operatorname{d\\!}x}I_{0}^{1-\alpha}f)(x)$
is the fractional derivative in the Riemann-Liouville sense for
$\alpha\in(0,1)$ whereas $(I_{0}^{\,\nu}f)(x)$ (given by Eq. (41)) is the
Riemann-Liouville fractional integral for $\nu\in(0,1)$. We point out that the
time operator in square bracket of (28) is equivalent to (RHilfer17, , Eq. (5)
for $\beta=0$). Acting with $I_{0}^{\,1}$ on both sides of Eq. (28) we get

$-\frac{1}{\tau_{1}}(I_{0}^{1}n)(t)=n(t)-1+\widetilde{\tau}\,(I^{1-\alpha}_{0}n)(t)$

represented also in the form

$n(t)=1-\dfrac{1}{\tau_{1}}\,(I^{1}_{0}\,n)(t)-\widetilde{\tau}\,(I^{1-\alpha}_{0}n)(t).$

That leads to Eq. (3) with
$B({\tau}_{1},{\tau}_{2},{\alpha})=\widetilde{\tau}=\tau_{2}^{\alpha}/\tau_{1}$
and

$M(t)=\tau_{2}^{-\alpha}+\frac{t^{-\alpha}}{\Gamma(1-\alpha)},$ (29)

which is interpreted as power-like smearing of $r(t,\tau)n(t)$ in Eq. (1). The
related Laplace transform becomes

$\widehat{M}(s)=\frac{\tau_{2}^{-\alpha}+s^{\alpha}}{s}.$ (30)

From the above and Eq. (16) we restore the characteristic exponent described
by Eq. (18).

###  IV.2 Coupled memories

The explicit form of the characteristic exponent $\widehat{\Psi}(s)$ enables
us to find memories $M(t)$ and $k(t)$ responsible for the time smearing of Eq.
(1). Recall that the memory $M(t)$ reflects the smearing of $n(t)$ whereas
$k(t)$ is related to the smearing of the time derivative $\dot{n}(t)$ and that
the memory $M(t)$ and its Laplace form are given by Eqs. (29) and (30). Using
the coupled pair $\widehat{M}(s)\widehat{k}(s)=1/s$ we find that
$\widehat{k}(s)$ and its Laplace form $k(t)$ yield

$\widehat{k}(s)=(\tau_{2}^{-\alpha}+s^{\alpha})^{-1}\quad\text{and}\quad
k(t)=t^{\alpha-1}E_{\alpha,\alpha}[-(t/\tau_{2})^{\alpha}].$ (31)

The singularity of $M(t)$ and $k(t)$ at $t=0$ is controlled by the parameter
$\alpha$. In the example quoted just below Corollary 4.1 in AHanyga20  it is
pointed out that $M(t)$ and $k(t)$ are the so–called Sonine functions and the
coupled pair $(k,M)$ is the Sonine pair (AHanyga20, , pp. 213–4); at a moment
we conclude that they are only Sonine functions $k$ and $M$. Such functions
are locally integrable non–decreasing functions which satisfy

$\sigma(t)\to\infty,\quad t\sigma(t)\to 0,\quad\text{for}\quad t\to
0;\quad\sigma\in\\{M,k\\}.$

Thus, (AHanyga20, , Theorem 3.1) is revealed. According to the philosophy of
the coupled me-
mories $M(t)$ is linked to Eq. (28) and $k(t)$ to

$\int_{0}^{t}(t-\xi)^{\alpha-1}E_{\alpha,\alpha}\Big{[}-\Big{(}\frac{t-\xi}{\tau_{2}}\Big{)}^{\alpha}\Big{]}~{}\dot{n}(\xi)\operatorname{d\\!}\xi=-\widetilde{\tau}\,n(t).$
(32)

Hence, the smearing of the relaxation function $n(t)$ can be changed into the
smearing of its first time derivative $\dot{n}(t)$ like it is done in Eq.
(32).

##  V The series form of solutions to (32)

General conditions of solvability Eq.(32) are precised in (ANKochubei11, ,
Theorem 2). It guarantees the uniqueness of the solution, its continuity,
differentiability, and completely monotone character on $(0,\infty)$. From Eq.
(31) the asymptotics of $\widehat{k}(s)$ turns out to be

$\displaystyle\begin{split}\widehat{k}(s)\to\tau_{2}^{\alpha},\quad
s~{}\widehat{k}(s)\to 0,\quad\text{for}\quad s\to 0,\\\ \widehat{k}(s)\to
0,\quad s~{}\widehat{k}(s)\to\infty,\quad\text{for}\quad
s\to\infty,\end{split}$

so we reconstruct the conditions listed in (ANKochubei11, , Theorem 2) except
of the first of them: $\widehat{k}(s)$ does not tend to infinity with $s\to 0$
but to the constant $\tau_{2}^{\alpha}$ instead. This clearly suggests the
existence of a solution to (32) which differs from (20).

Looking for the solution of (32) we apply (KGorska20, , Eq. (5)):

$n(t)=\mathcal{L}^{-1}[\widehat{k}(s)/(s~{}\widehat{k}(s)+\widetilde{\tau});t],$

in which we extract from denominator either $s~{}\widehat{k}(s)$ or
$\widetilde{\tau}$. This extraction procedure enables us to infer two kinds of
formulae. To derive them we employ the series form of $(1+x)^{-1}=\sum_{r\geq
0}(-x)^{r}$ for $|x|<1$, where we take either
$x=\widetilde{\tau}[s~{}\widehat{k}(s)]^{-1}$ or
$x=s~{}\widehat{k}(s)/\widetilde{\tau}$, getting the series form of solutions
$n_{\alpha}(\tau_{1},\tau_{2};t)$ and
$\widetilde{n}_{\alpha}(\tau_{1},\tau_{2};t)$, respectively. Clearly

$\displaystyle n_{\alpha}(\tau_{1},\tau_{2};t)$ $\displaystyle=\sum_{r\geq
0}(-\widetilde{\tau})^{r}\mathcal{L}^{-1}\Big{[}\dfrac{(\tau_{2}^{-\alpha}+s^{\alpha})^{r}}{s^{1+r}};t\Big{]},$
$\displaystyle\widetilde{n}_{\alpha}(\tau_{1},\tau_{2};t)$
$\displaystyle=\frac{1}{\widetilde{\tau}}\sum_{r\geq
0}(-\widetilde{\tau})^{-r}\mathcal{L}^{-1}\Big{[}\dfrac{s^{r}}{(\tau_{2}^{-\alpha}+s^{\alpha})^{1+r}};t\Big{]}.$

We point out that these formulae are equivalent to those obtained in
(KGorska19, , Eqs. (1.4), (3.1)) or (KGorska20, , Eqs. (6), (7))). The inverse
Laplace transforms present in $n_{\alpha}(\tau_{1},\tau_{2};t)$ and
$\widetilde{n}_{\alpha}(\tau_{1},\tau_{2};t)$ are calculated applying the
technique of (TRPrabhakar69, , Eq. (2.5)) exhibited in Eq. (38). For
$\alpha\in(0,1)$ we have

$n_{\alpha}(\tau_{1},\tau_{2};t)=\sum_{r\geq
0}(-\widetilde{\tau})^{r}\,t^{(1-\alpha)r}E_{\alpha,(1-\alpha)r+1}^{-r}[-(t/\tau_{2})^{\alpha}],$
(33)

while for $\alpha>1$

$\widetilde{n}_{\alpha}(\tau_{1},\tau_{2};t)=\frac{1}{\widetilde{\tau}}\sum_{r\geq
0}(-\widetilde{\tau})^{-r}\,t^{(\alpha-1)r}E_{\alpha,(\alpha-1)r+1}^{\,r}[-(t/\tau_{2})^{\alpha}].$
(34)

(Notice that Eq. (33) is the same as Eq. (23).) Both calculation procedures
are legitimate because the Mittag–Leffler functions setting in the series,
either (33) or (34), are well defined for all $r\in\mathbb{N}_{0}$ as
depending on the parameters $(1-\alpha)r+1>0;\alpha\in(0,1)$, and
$(\alpha-1)r+1>0;\alpha>1$, respectively. We point out that the Laplace
transforms in both solutions yield to
$(s+s^{\alpha}\,\widetilde{\tau}+\tau_{1}^{-1})^{-1}$, once for
$\alpha\in(0,1)$ and in turn for $\alpha>1$. In the excess wings model we have
$\alpha\in(0,1)$ so $n_{\alpha}(\tau_{1},\tau_{2};t)$ is the correct solution
for that range of $\alpha$. In the case of non-negative integer
$n\in\mathbb{N}_{0}$, the expression $E_{\alpha,\nu}^{-n}(z)$ becomes the
Mittag–Leffler polynomial of degree $\deg\big{(}E_{\alpha,\nu}^{-n}\big{)}=n$
which basic properties are quoted in A.

For $\alpha=1$ the solutions coincide,
$\widetilde{n}_{1}(\tau_{1},\tau_{2};t)=n_{1}(\tau_{1},\tau_{2};t)$, taking
the exponential decay form

$n_{1}(\tau_{1},\tau_{2};t)=\frac{\tau_{1}}{\tau_{1}+\tau_{2}}\exp\Big{(}-\frac{t}{\tau_{1}+\tau_{2}}\Big{)},$

viz., B. The equality of $n_{\alpha}(\tau_{1},\tau_{2};t)$ and
$\widetilde{n}_{\alpha}(\tau_{1},\tau_{2};t)$ for $\alpha\neq 1$ can be
established in the limit case of large $\tau_{2}$ for which the uniqueness
conditions of the initial Cauchy problem (32) and $n(0+)=1$, given in
ANKochubei11 , are satisfied. The limit of large $\tau_{2}$ means, assuming
$t$ be fixed, that the three parameter Mittag-Leffler function and the Mittag-
Leffler polynomial are considered for small values of their arguments being in
both cases equal to $t/\tau_{2}$. The asymptotic behaviour either of the
Mittag–Leffler function or of the associated $\deg(E^{-n}_{\alpha,\,\beta})$
Mittag-Leffler polynomial coincide for a small values of argument:

$E^{\gamma}_{\alpha,\,\beta}(x)\quad\text{or}\quad
E^{-n}_{\alpha,\,\beta}(x)\sim[\Gamma(\beta)]^{-1},\qquad x\to 0\,.$

Accordingly, when $\tau_{2}$ is growing and $t$ remains fixed, we deduce

$n_{\alpha}(\tau_{1},\tau_{2};t)\propto\sum_{r\geq
0}\frac{(-\widetilde{\tau}t^{1-\alpha})^{r}}{\Gamma[1+(1-\alpha)r]}=E_{1-\alpha}(-\widetilde{\tau}t^{1-\alpha}).$

Using the reciprocal arguments property (RGorenflo14, , Eq. (4.8.5))

$E_{-\nu}(z)+E_{\nu}(z^{-1})=1,\qquad\nu>0,\,z\in\mathbb{C}\setminus\\{0\\},$

we have

$n_{\alpha}(\tau_{1},\tau_{2};t)\propto
1-E_{\alpha-1}[-(\widetilde{\tau}t^{1-\alpha})^{-1}]=\sum_{r\geq
1}\frac{(-\widetilde{\tau}t^{1-\alpha})^{-r}}{\Gamma\big{(}1-(1-\alpha)r\big{)}},$

which is the asymptotics of $\widetilde{n}_{\alpha}(\tau_{1},\tau_{2};t)$ for
$\tau_{2}\to\infty$. Thus, we infer that $\tau_{2}\to\infty$ means
$\widehat{k}(s)\to\infty$ which confirms the equality
$n_{\alpha}(\tau_{1},\tau_{2};t)=\widetilde{n}_{\alpha}(\tau_{1},\tau_{2};t)$,
by bearing in mind the uniqueness of solution guaranteed by (ANKochubei11, ,
Theorem 2).

After some routine, but long and a little boring calculations employing
definitions of the Mittag-Leffler polynomials and the three parameter Mittag-
Leffler function, we get that the solutions $n_{\alpha}(\tau_{1},\tau_{2};t)$
and $\widetilde{n}_{\alpha}(\tau_{1},\tau_{2};t)$ can be presented in the form
of binomial Mittag-Leffler functions:

$\displaystyle n_{\alpha}(\tau_{1},\tau_{2};t)$
$\displaystyle=E_{(1,1-\alpha),1}(-t/\tau_{1},-\widetilde{\tau}t^{1-\alpha}),$
(35) $\displaystyle\widetilde{n}_{\alpha}(\tau_{1},\tau_{2};t)$
$\displaystyle=t^{\alpha-1}\,E_{(\alpha,\alpha-1),\alpha}[-(t/\tau_{2})^{\alpha},-t^{\,\alpha-1}/\,\widetilde{\tau}\,]/\,\widetilde{\tau}.$
(36)

We remind that the first of these results holds for $\alpha\in(0,1)$ while the
second one when $\alpha>1$. In turn, the equality
$n_{1}(\tau_{1},\tau_{2};t)=\widetilde{n}_{1}(\tau_{1},\tau_{2};t)$ yields
$E_{(1,1),1}(-t/\tau_{1},-\widetilde{\tau})=1/\widetilde{\tau}\,E_{(1,1),1}(-t/\tau_{2},-1/\,\widetilde{\tau}\,)$.
Reformulation of the Eqs. (33) and (34) by Eqs. (35) and (36) involves plenty
of technical details, first of all concerning transformations of finite and
infinite sums. All this goes beyond the presented exposition and is shifted to
C. Finally, we remark that $n_{\alpha}(\tau_{1},\tau_{2};t)$ in (35) can be
expressed as (25) and (26) or the formulae (RGarrappa16, , Eqs. (3.71),
(3.73)), whereas (36) coincides with (RGarrappa16, , Eq. (3.72)).

##  VI Conclusions

We have shown that the kinetic equations (2) and (3) assumed to govern the
relaxation phenomena and stemmed from the time smearing of either LHS or RHS
in non-Debye evolution equation $\dot{n}(t)=-r(t,\tau)n(t)$ determine their
stochastic interpretation. The crucial role in the presented approach is
played by the characteristic exponent $\widehat{\Psi}$ which provides us with
a bridge connecting kinetic equations and stochastic methods. Moreover, for a
large set of relaxing systems $\widehat{\Psi}$ obeys well-defined properties
which put it in the class of Bernstein functions and open new ways to push
forward mathematical and physical understanding of the relaxation phenomena.
To illustrate our methods we went beyond the family of the Havriliak-Negami
models and considered the excess wings model of relaxation. We identified the
characteristic exponent related to it and derived and solved kinetic equations
which reflect two ways of introducing the memory effects - the time smearing
of $\dot{n}(t)$ or $r(t,\tau)n(t)$ reflected in Eqns. (2) and (3),
respectively. Natural assumption that both approaches lead to the same
physical results allowed us to claim that the memory functions, $M(t)$ and
$k(t)$, responsible for both variants of smearing, form the Sonine pair, i.e.,
their transforms to the Laplace domain satisfy
$\widehat{M}(s)\widehat{k}(s)=1/s$. Results of the paper complete and, in a
sense, unify so-called deterministic and stochastic processes based
investigations of the non-Debye relaxation phenomena. We show that both these
approaches are not only mutually related but realize a correspondence
principle which joins different, but in fact equivalent views on the same
physical problem.

## Acknowledgments

K.G. and A.H. have been supported by the Polish National Center for Science
(NCN) research grant OPUS12 no. UMO-2016/23/B/ST3/01714. K. G. acknowledges
also support under the project Preludium Bis 2 no. UMO-2020/39/O/ST2/01563
awarded by the NCN and NAWA (Polish National Agency For Academic Exchange).
The research of T.K.P. has been supported in part by the University of Rijeka,
Croatia, under the project uniri-pr-prirod-19-16.

## Conflict of interest

The authors declare that they have no conflict of interest.

##  Appendix A Three parameter Mittag-Leffler function and Mittag-Leffler
polynomials

The three parameter Mittag-Leffler function is defined through the power
series (RGorenflo14, , p. 97, Eq. (5.1.1))

$E_{\alpha,\nu}^{\mu}(x)=\sum_{r\geq
0}\frac{(\mu)_{r}x^{r}}{r!\Gamma(\nu+\alpha r)},$ (37)

where
$\operatorname{\mathfrak{Re}}(\alpha),\operatorname{\mathfrak{Re}}(\nu),\operatorname{\mathfrak{Re}}(\mu)>0$
and $x\in\mathbb{R}$, $(\mu)_{r}$ denotes the familar Pochhammer symbol
(raising factorial) equal to
$\Gamma(\mu+r)/\Gamma(\mu)=\mu(\mu+1)\ldots(\mu+r-1);r\in\mathbb{N}_{0}$. The
Pochhammer symbol for $\mu=1$ is equal to $r!$ and Eq. (37) depends on two
parameters $\alpha$ and $\nu$ only. This case is named the two parameter
Mittag-Leffler (Wiman) function and it is quoted as
$E_{\alpha,\nu}(x)=E^{1}_{\alpha,\nu}(x)$. For $\mu=\nu=1$ Eq. (37) reduces to
the one parameter (standard) Mittag-Leffler function
$E_{\alpha}(x)=E_{\alpha,1}^{1}(x)$. The Laplace transform of
$t^{\nu-1}E_{\alpha,\nu}^{\mu}(\lambda t^{\alpha})$ equals

$\mathcal{L}[t^{\nu-1}E_{\alpha,\nu}^{\mu}(\lambda
t^{\alpha});s]=s^{-\nu}(1-\lambda s^{-\alpha})^{-\mu}\,$ (38)

for $\operatorname{\mathfrak{Re}}(\nu),\operatorname{\mathfrak{Re}}(s)>0$,
$|s|>|\lambda|^{1/\operatorname{\mathfrak{Re}}(\alpha)}$ TRPrabhakar69 .
Derivatives of the three parameter Mittag-Leffler function read

$x^{\,\nu-1}E_{\mu,\nu}^{\gamma}(ax^{\,\mu})=\frac{\operatorname{d\\!}}{\operatorname{d\\!}x}[x^{\nu}E_{\mu,\nu-1}^{\gamma}(ax^{\,\mu})]$
(39)

and

$x^{\,\nu-\alpha}E_{\mu,\nu-\alpha}^{\gamma}(ax^{\mu})=D_{x}^{\,\alpha}[x^{\nu}E_{\mu,\nu-1}^{\gamma}(ax^{\,\mu})],$
(40)

where
$({D_{x}^{\alpha}}f)(x)=(\frac{\operatorname{d\\!}}{\operatorname{d\\!}x}I_{0}^{1-\alpha}f)(x)$
is the fractional derivative in the Riemann-Liouville sense for
$\alpha\in(0,1)$ and

$(I_{0}^{\,\nu}f)(x)=\dfrac{1}{\Gamma(\nu)}\,\int_{0}^{x}(x-\xi)^{\nu-1}f(\xi)\operatorname{d\\!}\xi\,,\quad\nu\in(0,1],$
(41)

stands for the Riemann-Liouville fractional integral.

The Mittag-Leffler polynomials occur when the upper parameter in Eq. (37) is a
negative integer, i.e., $\gamma=-n,n\in\mathbb{N}_{0}$. From the definition of
the Pochhammer symbol all terms in Eq. (37) vanish when the upper parameter
$\gamma<-n$ and the series terminates leading to

$E_{\alpha,1+c}^{-n}(x)=\sum_{r=0}^{n}\binom{n}{r}\frac{(-x)^{r}}{\Gamma(1+c+\alpha
r)},\qquad\alpha,c>0.$ (42)

These objects are related to the Konhauser polynomials $Z^{c}_{n}(x;k)$
(JDEKonhauser67, , p. 304, Eq. (5)) defined via the formula

$Z^{c}_{n}(x;k)=\dfrac{\Gamma(kn+c+1)}{n!}\sum_{j=0}^{n}\binom{n}{j}\dfrac{(-x^{k})^{j}}{\Gamma(kn+c+1)},$

where $c>-1$. The latter extend the generalized (associated) Laguerre
polynomials $L_{n}^{(c)}(x^{k})=Z^{c}_{n}(x;k)$ (for the latter see below in
Appendix B), KGorska20 ; HMSrivastava82 . The connection formula between
Mittag–Leffler and Konhauser polynomials reads (Ozar, , p. 633, Eq. (7))

$E_{\alpha,c+1}^{-n}(x^{\alpha})=\dfrac{\Gamma(\alpha
n+c+1)}{n!}~{}Z^{c}_{n}(x;\alpha)\,.$

##  Appendix B The proof of
$n_{1}(\tau_{1},\tau_{1};t)=\widetilde{n}_{1}(\tau_{1},\tau_{1};t)$

The equality of $n_{1}(\tau_{1},\tau_{1};t)$ and
$\widetilde{n}_{1}(\tau_{1},\tau_{1};t)$ we can established using the
following three facts:

  1. 1.

$E_{1,1}^{-n}(x)=L_{n}(x)$ where $L_{n}(x)$ signifies the $n$th Laguerre
polynomial;

  2. 2.

the generating function for generalized (associated) Laguerre polynomials
$L_{n}^{(\alpha)}(x);\quad L_{n}^{(0)}(x)\equiv L_{n}(x)$, which Laplace
transform we use, reads (APPrudnikov-v2, , Eq. (5.11.2.1))

$\sum_{k\geq
0}t^{k}L_{k}^{\alpha}(x)=\dfrac{1}{(1-t)^{1+\alpha}}\,\exp\Big{(}\frac{tx}{t-1}\Big{)}$

for all $|t|<1$. For another generating functions see for instance Ref. Chat ;

  3. 3.

$E_{1,1}^{\,r}(-x)=\operatorname{e}^{-x}E_{1,1}^{-(r-1)}(x)$ which is Kummer’s
first transformation formula for the confluent hypergeometric function
${}_{1}F_{1}$, namely ${}_{1}F_{1}(r;1;-x)=E_{1,1}^{\,r}(x)$.

##  Appendix C Derivation of equations (35) and (36)

Substituting the Mittag-Leffler polynomial’s expression in Eq. (33) after some
algebra we conclude that

$n_{1}(\tau_{1},\tau_{2};t)=\sum_{r\geq
0}\sum_{j=0}^{r}\binom{r}{j}\frac{\Big{(}-\dfrac{\widetilde{\tau}}{t^{\alpha-1}}\Big{)}^{r-j}\Big{(}-\dfrac{t}{\tau_{1}}\Big{)}^{j}}{\Gamma\big{(}j+(1-\alpha)(r-j)+1\big{)}}.$

Setting $j=l_{1}$ and $r-j=l_{1}$ we can rewrite the right-hand side above as

$n_{1}(\tau_{1},\tau_{2};t)=\sum_{r\geq
0}{\sum_{\underset{l_{1}+l_{2}=r}{l_{1},l_{2}\geq
0}}}\frac{r!}{l_{1}!l_{2}!}\frac{\Big{(}-\dfrac{\widetilde{\tau}}{t^{\alpha-1}}\Big{)}^{l_{2}}\Big{(}-\dfrac{t}{\tau_{1}}\Big{)}^{l_{1}}}{\Gamma\big{(}1+l_{1}+(1-\alpha)l_{2}\big{)}}.$

Comparison with Eq. (21) gives Eq. (35). Analogous calculation can be done for
Eq. (36); during these computations we use the series form of the three
parameter Mittag-Leffler function.

## References

  * (1) Anderssen RS, Loy RJ. Completely monotone fading memory relaxation moduli. Bull Austral Math Soc 2002; 65:449
  * (2) Anderssen RS, Loy RJ. Rheological implications of completely monotone fading memory. J Rheol 2002; 46:1459
  * (3) Berg C. Stieltjes-Pick-Bernstein-Schoenberg and their connection to completely monotonicity. In: Mateu J and Porcu E, editors. Positive Define Functions: From Schoenberg to Space-Time challenges. Dep. Math. of Univ. Jaume I, Castellon; 2008
  * (4) Bochner S. Harmonic Analysis and the Theory of Probability. Univ. of California Press, Berkeley/Los Angeles; 1955
  * (5) Brand R, Lunkenheimer P, Schneider U, Loidl A. Excess wing in the dielectric loss of glass-forming ethanol: A relaxation process. Phys Rev B 2000; 62:8878
  * (6) Capelas de Oliveira E, Mainardi F, Vaz Jr J. Models based on Mittag-Leffler functions for anomalous relaxation in dielectrics. Eur Phys J Special Topics 2011; 193:161
  * (7) Chatterjea SK. On a generating function of Laguerre polynomials. Boll Un Mat Ital 1962; 17:179
  * (8) Chechkin AV, Sokolov IM. On relation between generalized diffusion and subordination schemes. Phys Rev E 2021; 103:032133
  * (9) Dixon PK. Specific-heat spectroscopy and dielectric susceptibility measurements of $salol$ at the glass transition. Phys Rev B 1990; 42:8179
  * (10) Dixon PK, Wu L, Nagel SR, Williams BD, Carini JP. Scaling in the relaxation of supercooled liquids. Phys Rev Lett 1990; 65:1108
  * (11) Fogedby HC, Langevin equations for continuous time Lévy flights. Phys Rev E 1994; 50:1657
  * (12) Garrappa R, Mainardi F, Maione G. Models of dielectric relaxation based on completely monotone functions. Frac Calc Appl Anal 2016; 19:1105 ; corrected version available in arXiv: 1611.04028
  * (13) Giusti A, Colombaro I, Garra R, Garrappa R, Polito F, Popolizio M, Mainardi F. A practical guide to Prabhakar fractional calculus. Frac Calc Appl Anal 2020; 23:9
  * (14) Giusti A. General fractional calculus and Prabhakara’s theory. Comm Nonlinear Sci Numer Simulat 2020; 83:105114
  * (15) Glöckle WG, Nonnenmacher TF. Fox function representation of non-Debye relaxation processes. J Stat Phys 1993; 71:741
  * (16) Gorenflo R, Kilbas AA, Mainardi F, Rogosin SV. Mittag-Leffler Functions, Related Topics and Applications. Springer, New York; 2014
  * (17) Górska K, Horzela A, Bratek Ł, Penson KA, Dattoli G. The Havriliak-Negami relaxation and its relatives: the response, relaxation and probability density functions. J Phys A: Theor Math 2018; 51:135202
  * (18) Górska K, Horzela A, Pogány TK. A note on the article ”Anomalous relaxation model based on the fractional derivative with a Prabhakar-like kernel” $[$Z. Angew. Math. Phys. $(2019)~{}{70}:42]$. Z Angew Math Phys 2019; 70:141
  * (19) Górska K, Horzela A. The Volterra type equation related to the non-Debye relaxation. Comm Nonlinear Sci Numer Simulat 2020; 85:105246
  * (20) Górska K, Horzela A, Lattanzi A, Pogány TK. On the complete monotonicity of the three parameter generalized Mittag-Leffler function $E_{\alpha,\beta}^{\gamma}(-x)$. Appl Anal Discret Math 2021; 15:118
  * (21) Górska K, Horzela A. Non-Debye Relaxations: Two types of memories and their Stieltjes character, Mathematics 2021; 9:477
  * (22) Górska K, Horzela A, Pogány TK. Non-Debye relaxations: smeared time evolution, memory effects, and the Laplace exponents, Comm Nonlinear Sci Numer Simulat 2021; 99:105837
  * (23) Grippenberg G, Londen SO, Staffans OJ. Volterra Integral and Functional Equations. Cambridge University Press, Cambridge; 1990
  * (24) Hanyga A, Seredyńska M. On a Mathematical Framework for the Constitutive Equations of Anisotropic Dielectric Relaxation. J Stat Phys 2008; 131:269
  * (25) Hanyga A. A comment on a controversial issue: A generalized fractional derivative cannot have a regular kernel. Frac Calc Appl Anal 2020; 23:211
  * (26) Hilfer R. Fitting the excess wing in the dielectric $\alpha$-relaxation of propylene carbonate. J Phys: Condens Matter 2002; 14:2297
  * (27) Hilfer R. Experimental evidence for fractional time evolution in glass forming materials. Chem Phys 2002; 284:399
  * (28) Hilfer R. Analytical representations for relaxation functions of glasses. J Non-Cryst Solids 2002; 305:122
  * (29) Hilfer R. $H$-function representations for stretched exponential relaxation and non-Debye susceptibilities in glassy systems. Phys Rev E 2002; 65:061510
  * (30) Hilfer R. Composite continuous time random walks. Eur Phys J B 2017; 90:233
  * (31) Jonscher AK. The universal dielectric response and its physical significance. IEEE Transactions on Electrical Insulation 1992; 27:407
  * (32) Kochubei AN. General fractional calculus, evolution equations, and renewal processes. Integr Equ Oper Theory 2011; 71:583
  * (33) Konhauser JDE. Biorthogonal polynomials suggested by the Laguerre polynomials. Pacific J Math 1967; 21:303
  * (34) Nigmatulin RR, Khamzin AA, Baleanu D. On the Laplace integral representation of multivariable Mittag-Leffler functions in anomalous relaxation. Math Meth Appl Sci 2016; 39:2983
  * (35) Mainardi F, Garrappa R. On complete monotonicity of the Prabhakar function and non-Debye relaxation in dielectrics. J Comp Phys 2015; 293:70
  * (36) Mathai AM, Haubold HJ. Special functions for applied scientists. Springer, New York; 2008
  * (37) Özarslan MA, Kürt C. Bivariate Mittag–Leffler functions arising in the solutions of convolution integral equation with $2D$–Laguerre–Konhauser polynomials in the kernel. Appl Math Comput 2019; 347:631
  * (38) Prabhakar TR. A singular integral equation with a generalized Mittag Leffler function in the kernel. Yokohoma Math J 1971; 19:7
  * (39) Prudnikov AP, Brychkov YuA, Marichev OI. Integrals and Series. Special Functions. Vol. 2. Gordon and Breach, Amsterdam; 1998
  * (40) Sandev T, Tomovski Ž. Fractional Equations and Models. Theory and Applications. Springer, New York; 2019
  * (41) Schilling RL. An introduction to Lévy and Feller processes. In: From Lévy–type processes to parabolic SPDEs. Adv. Courses Math., pp. 1–126. Birkhäuser – Springer, Cham (2016)
  * (42) Schilling RL, Song R, Vondraček Z. Bernstein Functions. De Gruyter, Berlin; 2010
  * (43) Srivastava HM. Some biorthogonal polynomials suggested by the Laguerre polynomials. Pacific J Math 1982; 98:235
  * (44) Stanislavsky A, Weron K, Weron A. Anomalous diffusion approach to non-exponential relaxation in complex physical systems. Comm Nonlinear Sci Numer Simulat 2015; 24:117
  * (45) Stanislavsky A, Weron K. Atypical case of the dielectric relaxation responses and its fractional kinetic equation. Frac Calc Appl Math 2016; 19:212
  * (46) Stanislavsky A, Weron K. Stochastic tools hidden behind the empirical dielectric relaxation laws. Rep Prog Phys 2017; 80:036001
  * (47) Stanislavsky A, Weron K. Fractional-calculus tools applied to study the nonexponential relaxation in dielectrics. In: Tarasov VE, editor. Handbook of fractional calculus with applications in physics, Part B, Vol. 5. De Gruyter, Berlin; 2019. p. 53–70.
  * (48) Stanislavsky A, Weron K. Duality in fractional systems, Comm Nonlinear Sci Numer Simulat 2021; 101:105861
  * (49) Tomovski Ž, Pogány TK, Srivastava HM. Laplace type integral expression for a certain three-parameter family of generalized Mittag-Leffler functions with application involving complete monotonicity. J Franklin Inst 2014; 351:5437
  * (50) Details of these distributions are irrelevant, however for applications to the relaxation phenomena it is usually assumed that we deal with $\alpha$-stable distributions; final results come from the generalized limit theorems.
  * (51) From the probabilistic point of view the Bernstein functions may be identified as subordinators and thus it cannot be strange that they play the crucial role in stochastic analysis of relaxation and anomalous diffusion.
  * (52) CMFs provide us also with an example of the fading memory concept proposed by L. Boltzmann and reintroduced to physics through applications in rheology and elasticity theory, see e.g. RSAnderssen02 ; RSAnderssen02a .
  * (53) Throughout the paper the superscript $\,\widehat{}\,$ denotes the Laplace transform: $\widehat{h}(s)=\mathcal{L}[h(t);s]=\int_{0}^{\infty}\operatorname{e}^{-st}h(t)\operatorname{d\\!}t$; $h(t)$ is the inverse Laplace transform given as $h(t)=\mathcal{L}^{-1}[\widehat{h}(s);t]=\int_{L}\operatorname{e}^{st}\widehat{h}(s)\operatorname{d\\!}s/(2\operatorname{i}\\!\pi)$ where $L$ is a Bromwich contour which leaves all singularities of $\widehat{h}(s)$ left to it.
  * (54) Main properties of the three parameter Mittag–Leffler function $E_{\mu,\nu}^{\gamma}(z)$ are listed in A.