* [73] So Hirata and Martin Head-Gordon.  Time-dependent density functional theory within the Tamm–Dancoff approximation.  Chem. Phys. Lett., 314(3-4):291–299, 1999.
  * [74] Florian Buchholz, Iris Theophilou, Soeren EB Nielsen, Michael Ruggenthaler, and Angel Rubio.  Reduced density-matrix approach to strong matter-photon interaction.  ACS Photonics, 6(11):2694–2711, 2019.
  * [75] Tor S Haugland, Enrico Ronca, Eirik F Kjønstad, Angel Rubio, and Henrik Koch.  Coupled cluster theory for molecular polaritons: Changing ground and excited states.  Phys. Rev. X, 10(4):041043, 2020.
  * [76] Florian Buchholz, Iris Theophilou, Klaas J. H. Giesbertz, Michael Ruggenthaler, and Angel Rubio.  Light–matter hybrid-orbital-based first-principles methods: The influence of polariton statistics.  J. Chem. Theory Comput., 16(9):5601–5620, 2020.  PMID: 32692551.
  * [77] Uliana Mordovina, Callum Bungey, Heiko Appel, Peter J. Knowles, Angel Rubio, and Frederick R. Manby.  Polaritonic coupled-cluster theory.  Phys. Rev. Res., 2:023262, Jun 2020.
  * [78] Jonathan McTague and Jonathan J Foley.  Non-Hermitian cavity quantum electrodynamics–configuration interaction singles approach for polaritonic structure with ab initio molecular Hamiltonians.  J. Chem. Phys., 156(15), 2022.
  * [79] Pedro Miguel M. C. de Melo and Andrea Marini.  Unified theory of quantized electrons, phonons, and photons out of equilibrium: A simplified ab initio approach based on the generalized baym-kadanoff ansatz.  Phys. Rev. B, 93:155102, Apr 2016.
  * [80] Pekka Pyykkö.  Relativistic effects in chemistry: more common than you thought.  Annu. Rev. Phys. Chem., 63:45–64, 2012.
  * [81] Pekka Pyykko and Jean Paul Desclaux.  Relativity and the periodic system of elements.  Acc. Chem. Res., 12(8):276–281, 1979.
  * [82] P Romaniello and PL de Boeij.  Relativistic two-component formulation of time-dependent current-density functional theory: Application to the linear response of solids.  J. Chem. Phys., 127(17):174111, 2007.
  * [83] Florent Calvo, Elke Pahl, Michael Wormit, and Peter Schwerdtfeger.  Evidence for low-temperature melting of mercury owing to relativity.  Angew. Chem. Int. Ed., 52(29):7583–7585, 2013.
  * [84] Rajeev Ahuja, Andreas Blomqvist, Peter Larsson, Pekka Pyykkö, and Patryk Zaleski-Ejgierd.  Relativity and the lead-acid battery.  Phys. Rev. Lett., 106(1):018301, 2011.
  * [85] V. P. Kosheleva, A. V. Volotka, D. A. Glazov, D. V. Zinenko, and S. Fritzsche.  $g$ factor of lithiumlike silicon and calcium: Resolving the disagreement between theory and experiment.  Phys. Rev. Lett., 128:103001, Mar 2022.
  * [86] D. A. Glazov, F. Köhler-Langes, A. V. Volotka, K. Blaum, F. Heiße, G. Plunien, W. Quint, S. Rau, V. M. Shabaev, S. Sturm, and G. Werth.  $g$ factor of lithiumlike silicon: New challenge to bound-state qed.  Phys. Rev. Lett., 123:173001, Oct 2019.
  * [87] Robert Laskowski and Peter Blaha.  Understanding the l 2, 3 x-ray absorption spectra of early 3 d transition elements.  Phys. Rev. B, 82(20):205104, 2010.
  * [88] Nanna Holmgaard List, Trond Saue, and Patrick Norman.  Rotationally averaged linear absorption spectra beyond the electric-dipole approximation.  Mol. Phys., 115(1-2):63–74, June 2016.
  * [89] Marta L Vidal, Pavel Pokhilko, Anna I Krylov, and Sonia Coriani.  Equation-of-motion coupled-cluster theory to model L-edge x-ray absorption and photoelectron spectra.  J. Phys. Chem. Lett., 11(19):8314–8321, 2020.
  * [90] Joseph M. Kasper, Torin F. Stetina, Andrew J. Jenkins, and Xiaosong Li.  Ab initio methods for L-edge x-ray absorption spectroscopy.  Chem. Phys. Rev., 1(1):011304, December 2020.
  * [91] Lukas Konecny, Jan Vicha, Stanislav Komorovsky, Kenneth Ruud, and Michal Repisky.  Accurate X-ray Absorption Spectra near L- and M-Edges from Relativistic Four-Component Damped Response Time-Dependent Density Functional Theory.  Inorg. Chem., 61(2):830–846, 2022.
  * [92] Peter Hrobarik, Veronika Hrobarikova, Florian Meier, Michal Repisky, Stanislav Komorovsky, and Martin Kaupp.  Relativistic four-component dft calculations of 1h nmr chemical shifts in transition-metal hydride complexes: unusual high-field shifts beyond the buckingham–stephens model.  J. Phys. Chem. A, 115(22):5654–5659, 2011.
  * [93] Peter Hrobarik, Veronika Hrobarikova, Anja H Greif, and Martin Kaupp.  Giant spin-orbit effects on nmr shifts in diamagnetic actinide complexes: Guiding the search of uranium (vi) hydride complexes in the correct spectral range.  Angew. Chem. Int. Ed., 51(43):10884–10888, 2012.
  * [94] Jan Vicha, Radek Marek, and Michal Straka.  High-frequency^13C and ^29Si NMR chemical shifts in diamagnetic low-valence compounds of Tl^I and Pb^II: Decisive role of relativistic effects.  Inorg. Chem., 55(4):1770–1781, January 2016.
  * [95] Jan Vicha, Jan Novotny, Stanislav Komorovsky, Michal Straka, Martin Kaupp, and Radek Marek.  Relativistic heavy-neighbor-atom effects on nmr shifts: Concepts and trends across the periodic table.  Chem. Rev., 120(15):7065–7103, 2020.
  * [96] Bob Martin and Jochen Autschbach.  Temperature dependence of contact and dipolar NMR chemical shifts in paramagnetic molecules.  J. Chem. Phys., 142(5), February 2015.
  * [97] Jan Novotny, Martin Sojka, Stanislav Komorovsky, Marek Necas, and Radek Marek.  Interpreting the paramagnetic NMR spectra of potential ru(III) metallodrugs: Synergy between experiment and relativistic DFT calculations.  J. Am. Chem. Soc., 138(27):8432–8445, July 2016.
  * [98] Arobendo Mondal, Michael W. Gaultois, Andrew J. Pell, Marcella Iannuzzi, Clare P. Grey, Jürg Hutter, and Martin Kaupp.  Large-scale computation of nuclear magnetic resonance shifts for paramagnetic solids using CP2k.  J. Chem. Theory Comp., 14(1):377–394, December 2017.
  * [99] Peter John Cherry, Syed Awais Rouf, and Juha Vaara.  Paramagnetic enhancement of nuclear spin–spin coupling.  J. Chem. Theory Comput., 13(3):1275–1283, February 2017.
  * [100] Stanislav Komorovsky.  Relativistic theory of EPR and (p)NMR.  In Reference Module in Chemistry, Molecular Sciences and Chemical Engineering. Elsevier, 2023.
  * [101] Irina Malkin, Olga L. Malkina, Vladimir G. Malkin, and Martin Kaupp.  Relativistic two-component calculations of electronic g-tensors that include spin polarization.  J. Chem. Phys., 123(24), December 2005.
  * [102] Peter Hrobarik, Michal Repisky, Stanislav Komorovsky, Veronika Hrobarikova, and Martin Kaupp.  Assessment of higher-order spin–orbit effects on electronic g-tensors of $d^{1}$ transition-metal complexes by relativistic two-and four-component methods.  Theor. Chem. Acc., 129:715–725, 2011.
  * [103] Sebastian Gohr, Peter Hrobarik, Michal Repisky, Stanislav Komorovsky, Kenneth Ruud, and Martin Kaupp.  Four-component relativistic density functional theory calculations of epr g-and hyperfine-coupling tensors using hybrid functionals: validation on transition-metal complexes with large tensor anisotropies and higher-order spin–orbit effects.  J. Phys. Chem. A, 119(51):12892–12905, 2015.
  * [104] Hélène Bolvin and Jochen Autschbach.  Relativistic methods for calculating electron paramagnetic resonance (EPR) parameters.  In Handbook of Relativistic Quantum Chemistry, pages 725–763. Springer Berlin Heidelberg, October 2016.
  * [105] Debora Misenkova, Florian Lemken, Michal Repisky, Jozef Noga, Olga L. Malkina, and Stanislav Komorovsky.  The four-component DFT method for the calculation of the EPR g-tensor using a restricted magnetically balanced basis and london atomic orbitals.  J. Chem. Phys., 157(16), October 2022.
  * [106] Ben Joseph R. Cuyacot, Jan Novotny, Raphael J. F. Berger, Stanislav Komorovsky, and Radek Marek.  Relativistic spin–orbit electronegativity and the chemical bond between a heavy atom and a light atom.  Chem. Eur. J., 28(24), April 2022.
  * [107] Taye B. Demissie, Brady D. Garabato, Kenneth Ruud, and Pawel M. Kozlowski.  Mercury methylation by cobalt corrinoids: Relativistic effects dictate the reaction mechanism.  Angew. Chem. Int. Ed., 55:11503–11506, 2016.
  * [108] Torsha Moitra, Pijush Karak, Sayantani Chakraborty, Kenneth Ruud, and Swapan Chakrabarti.  Behind the scenes of spin-forbidden decay pathways in transition metal complexes.  Phys. Chem. Chem. Phys., 23(1):59–81, 2021.
  * [109] V Pershina.  Electronic structure and properties of superheavy elements.  Nucl. Phys. A, 944:578–613, 2015.
  * [110] S.A. Giuliani, Z. Matheson, W. Nazarewicz, E. Olsen, P.-G. Reinhard, J. Sadhukhan, B. Schuetrumpf, N. Schunck, and P. Schwerdtfeger.  Colloquium: Superheavy elements: Oganesson and beyond.  Rev. Mod. Phys., 91(1), January 2019.
  * [111] Marius Kadek, Baokai Wang, Marc Joosten, Wei-Chi Chiu, Francois Mairesse, Michal Repisky, Kenneth Ruud, and Arun Bansil.  Band structures and Z2 invariants of two-dimensional transition metal dichalcogenide monolayers from fully relativistic dirac-kohn-sham theory using gaussian-type orbitals.  Phys. Rev. Mater., 7(6):064001, 2023.
  * [112] EI Rashba and E Ya Sherman.  Spin-orbital band splitting in symmetric quantum wells.  Phys. Lett. A, 129(3):175–179, 1988.
  * [113] Martin Gmitra, Sergej Konschuh, Christian Ertler, Claudia Ambrosch-Draxl, and Jaroslav Fabian.  Band-structure topologies of graphene: Spin-orbit coupling effects from first principles.  Phys. Rev. B, 80(23):235431, 2009.
  * [114] Zhiyong Y Zhu, Yingchun C Cheng, and Udo Schwingenschlögl.  Giant spin-orbit-induced spin splitting in two-dimensional transition-metal dichalcogenide semiconductors.  Phys. Rev. B, 84(15):153402, 2011.
  * [115] Andreas Hermann, Jürgen Furthmüller, Heinz W Gäggeler, and Peter Schwerdtfeger.  Spin-orbit effects in structural and electronic properties for the solid state of the group-14 elements from carbon to superheavy element 114.  Phys. Rev. B, 82(15):155116, 2010.
  * [116] K. G. Dyall and K. Faegri Jr.  Introduction to Relativistic Quantum Chemistry.  Oxford University Press, New York, 2007.
  * [117] Markus Reiher and Alexander Wolf.  Relativistic Quantum Chemistry: The Fundamental Theory of Molecular Science; 2nd Edition.  Wiley-VCH, 2014.
  * [118] Trond Saue.  Relativistic Hamiltonians for chemistry: A primer.  ChemPhysChem, 12(17):3077–3094, 2011.
  * [119] Matthew S. Kelley and Toru Shiozaki.  Large-scale Dirac–Fock–Breit method using density fitting and 2-spinor basis functions.  J. Chem. Phys., 138(20), May 2013.
  * [120] Shichao Sun, Jordan Ehrman, Tianyuan Zhang, Qiming Sun, Kenneth G Dyall, and Xiaosong Li.  Scalar breit interaction for molecular calculations.  J. Chem. Phys., 158(17), 2023.
  * [121] Chad E Hoyer, Lixin Lu, Hang Hu, Kirill D Shumilov, Shichao Sun, Stefan Knecht, and Xiaosong Li.  Correlated dirac–coulomb–breit multiconfigurational self-consistent-field methods.  J. Chem. Phys., 158(4), 2023.
  * [122] BJ Mackenzie, IP Grant, and PH Norrington.  Program to calculate transverse breit and qed corrections to energy levels in a multiconfiguration dirac-fock environment.  Comput. Phys. Commun.;(Netherlands), 21(2), 1980.
  * [123] VM Shabaev, II Tupitsyn, and VA Yerokhin.  Model operator approach to the lamb shift calculations in relativistic many-electron atoms.  Phys. Rev. A, 88(1):012513, 2013.
  * [124] Gustavo A Aucar.  Toward a QFT-based theory of atomic and molecular properties.  Phys. Chem. Chem. Phys., 16(10):4420–4438, 2014.
  * [125] Peter Schwerdtfeger, Lukas F Pašteka, Andrew Punnett, and Patrick O Bowman.  Relativistic and quantum electrodynamic effects in superheavy elements.  Nucl. Phys. A, 944:551–577, 2015.
  * [126] LF Pašteka, E Eliav, A Borschevsky, U Kaldor, and P Schwerdtfeger.  Relativistic coupled cluster calculations with variational quantum electrodynamics resolve the discrepancy between experiment and theory concerning the electron affinity and ionization potential of gold.  Phys. Rev. Lett., 118(2):023002, 2017.
  * [127] AV Malyshev, DA Glazov, VM Shabaev, II Tupitsyn, VA Yerokhin, and VA Zaytsev.  Model-qed operator for superheavy elements.  Phys. Rev. A, 106(1):012806, 2022.
  * [128] Ayaki Sunaga, Maen Salman, and Trond Saue.  4-component relativistic Hamiltonian with effective QED potentials for molecular calculations.  arXiv:2207.14101v1.
  * [129] Nobuki Inoue, Yoshihiro Watanabe, and Haruyuki Nakano.  Theoretical examination of QED Hamiltonian in relativistic molecular orbital theory.  J. Chem. Phys., 159(5), 2023.
  * [130] Werner Kutzelnigg and Wenjian Liu.  Quasirelativistic theory equivalent to fully relativistic theory.  J. Chem. Phys., 123(24):241102, December 2005.
  * [131] Wenjian Liu and Werner Kutzelnigg.  Quasirelativistic theory. II. theory at matrix level.  J. Chem. Phys., 126(11):114107, March 2007.
  * [132] Miroslav Ilias and Trond Saue.  An infinite-order two-component relativistic hamiltonian by a simple one-step transformation.  J. Chem. Phys., 126(6):064102, February 2007.
  * [133] Stefan Knecht, Michal Repisky, Hans Jørgen Aagaard Jensen, and Trond Saue.  Exact two-component hamiltonians for relativistic quantum chemistry: Two-electron picture-change corrections made simple.  J. Chem. Phys., 157(11):114106, 2022.
  * [134] Can Liao, Joseph M Kasper, Andrew J Jenkins, Ping Yang, Enrique R Batista, Michael J Frisch, and Xiaosong Li.  State interaction linear response time-dependent density functional theory with perturbative spin–orbit coupling: Benchmark and perspectives.  JACS Au, 3(2):358–367, 2023.
  * [135] Thomas Fransson, Trond Saue, and Patrick Norman.  Four-component damped density functional response theory study of UV/Vis absorption spectra and phosphorescence parameters of group 12 metal-substituted porphyrins.  J. Chem. Theory Comput., 12(5):2324–2334, 2016.
  * [136] Hans Ågren, Olav Vahtras, and Boris Minaev.  Response theory and calculations of spin-orbit coupling phenomena in molecules.  In Advances in Quantum Chemistry, volume 27, pages 71–162. Elsevier, 1996.
  * [137] Clàudia Climent, David Casanova, Johannes Feist, and Francisco J Garcia-Vidal.  Not dark yet for strong light-matter coupling to accelerate singlet fission dynamics.  Cell Rep., 3(4), 2022.
  * [138] Elad Eizner, Luis A Martínez-Martínez, Joel Yuen-Zhou, and Stéphane Kéna-Cohen.  Inverting singlet and triplet excited states using strong light-matter coupling.  Sci. Adv., 5(12):eaax4482, 2019.
  * [139] S Kéna-Cohen and SR Forrest.  Green polariton photoluminescence using the red-emitting phosphor PtOEP.  Phys. Rev. B, 76(7):075202, 2007.
  * [140] Manuel Hertzog, Mao Wang, Jürgen Mony, and Karl Börjesson.  Strong light–matter interactions: a new direction within chemistry.  Chem. Soc. Rev., 48(3):937–961, 2019.
  * [141] Xiaodong Xu, Wang Yao, Di Xiao, and Tony F Heinz.  Spin and pseudospins in layered transition metal dichalcogenides.  Nat. Phys., 10(5):343–350, 2014.
  * [142] See-Hun Yang, Ron Naaman, Yossi Paltiel, and Stuart SP Parkin.  Chiral spintronics.  Nat. Rev. Phys., 3(5):328–343, 2021.
  * [143] Herbert Spohn.  Dynamics of charged particles and their radiation field.  Cambridge university press, 2004.
  * [144] Fumio Hiroshima.  Ground states of quantum field models: perturbation of embedded eigenvalues, volume 35.  Springer, 2019.
  * [145] Lewis H Ryder.  Quantum field theory.  Cambridge university press, 1996.
  * [146] Franz Mandl and Graham Shaw.  Quantum field theory.  John Wiley & Sons, 2010.
  * [147] Toshimitsu Takaesu.  On the spectral analysis of quantum electrodynamics with spatial cutoffs. I.  J. Mat. Phys., 50(6), 2009.
  * [148] John C Baez, Irving E Segal, and Zhengfang Zhou.  Introduction to algebraic and constructive quantum field theory, volume 47.  Princeton University Press, 2014.
  * [149] W. Greiner and J. Reinhardt.  Field Quantization.  Springer, 1996.
  * [150] M. E. Peskin and D. V. Schroeder.  An Introduction to Quantum Field Theory.  Westview Press, 1995.
  * [151] JOACHIM Reinhardt and WALTER Greiner.  Quantum electrodynamics of strong fields.  Rep. Prog. Phys., 40(3):219, 1977.
  * [152] Sølve Selstø, Eva Lindroth, and Jakob Bengtsson.  Solution of the dirac equation for hydrogenlike systems exposed to intense electromagnetic pulses.  Phys. Rev. A, 79, 4 2009.
  * [153] D.-A. Deckert, D. Dürr, F. Merkl, and M. Schottenloher.  Time-evolution of the external field problem in Quantum Electrodynamics.  J. Math. Phys., 51(12):122301, 12 2010.
  * [154] V. Rokaj, D. M. Welakuh, M. Ruggenthaler, and A. Rubio.  Light–matter interaction in the long-wavelength limit: No ground-state without dipole self-energy.  J. Phys. B: At., Mol. Opt. Phys., 51:034005, 2018.
  * [155] Mark Kamper Svendsen, Kristian Sommer Thygesen, Angel Rubio, and Johannes Flick.  Molecules in real cavities with quantum electroynamical density functional theory.  arXiv preprint arXiv:2305.02391, 2023.
  * [156] C. Schäfer, M. Ruggenthaler, V. Rokaj, and A. Rubio.  Relevance of the quadratic diamagnetic and self-polarization terms in cavity quantum electrodynamics.  ACS Photonics, 7:975–990, 2020.
  * [157] F Gesztesy, B Thaller, and Harald Grosse.  Efficient method for calculating relativistic corrections for spin-$1/2$ particles.  Phys. Rev. Lett., 50(9):625, 1983.
  * [158] Jürg Fröhlich and Urban M Studer.  Gauge invariance and current algebra in nonrelativistic many-body theory.  Rev. Mod. Phys., 65(3):733, 1993.
  * [159] Dominik Sidler, Thomas Schnappinger, Anatoly Obzhirov, Michael Ruggenthaler, Markus Kowalewski, and Angel Rubio.  Unraveling a cavity induced molecular polarization mechanism from collective vibrational strong coupling.  arXiv preprint arXiv:2306.06004, 2023.
  * [160] Thomas Schnappinger, Dominik Sidler, Michael Ruggenthaler, Angel Rubio, and Markus Kowalewski.  Cavity Born–Oppenheimer Hartree–Fock Ansatz: Light–Matter Properties of Strongly Coupled Molecular Ensembles.  J. Phys. Chem. Lett., 14:8024–8033, 2023.
  * [161] Carsten A Ullrich.  Time-dependent density-functional theory: concepts and applications.  OUP Oxford, R, 2011.
  * [162] Lukas Konecny, Michal Repisky, Kenneth Ruud, and Stanislav Komorovsky.  Relativistic four-component linear damped response TDDFT for electronic absorption and circular dichroism calculations.  J. Chem. Phys., 151(19):194112, 2019.
  * [163] Roi Baer and Leeor Kronik.  Time-dependent generalized Kohn–Sham theory.  Eur. Phys. J. B, 91:1–9, 2018.
  * [164] Mark E Casida.  Time-dependent density functional response theory for molecules.  In Recent Advances In Density Functional Methods: (Part I), pages 155–192. World Scientific, 1995.
  * [165] Mark E Casida.  Time-dependent density-functional theory for molecules and molecular solids.  J. Mol. Struct.: THEOCHEM, 914(1-3):3–18, 2009.
  * [166] P. Norman, K. Ruud, and T. Saue.  Principles and Practices of Molecular Properties.  Wiley-VCH, Chichester, 2018.
  * [167] Dominik Sidler, Christian Schäfer, Michael Ruggenthaler, and Angel Rubio.  Polaritonic chemistry: Collective strong coupling implies strong local modification of chemical properties.  J. Phys. Chem. Lett., 12(1):508–516, 2020.
  * [168] Michal Repisky, Stanislav Komorovsky, Marius Kadek, Lukas Konecny, Ulf Ekström, Elena Malkin, Martin Kaupp, Kenneth Ruud, Olga L. Malkina, and Vladimir G. Malkin.  ReSpect: Relativistic spectroscopy DFT program package.  J. Chem. Phys., 152(18):184101, May 2020.
  * [169] Stanislav Komorovsky, Peter J Cherry, and Michal Repisky.  Four-component relativistic time-dependent density-functional theory using a stable noncollinear DFT ansatz applicable to both closed-and open-shell systems.  J. Chem. Phys., 151(18):184111, 2019.
  * [170] Ernest R. Davidson.  The iterative calculation of a few of the lowest eigenvalues and corresponding eigenvectors of large real-symmetric matrices.  J. Comput. Phys., 17(1):87–94, 1975.
  * [171] Jeppe Olsen, Poul Jørgensen, and Jack Simons.  Passing the one-billion limit in full configuration-interaction (FCI) calculations.  Chem. Phys. Lett., 169(6):463–472, 1990.
  * [172] Trond Saue and HJ Aa Jensen.  Linear response at the 4-component relativistic level: Application to the frequency-dependent dipole polarizabilities of the coinage metal dimers.  J. Chem. Phys., 118(2):522–536, 2003.
  * [173] Radovan Bast, Hans Jørgen Aa Jensen, and Trond Saue.  Relativistic adiabatic time-dependent density functional theory using hybrid functionals and noncollinear spin magnetization.  Int. J. Quantum Chem., 109(10):2091–2112, 2009.
  * [174] Richard E Stanton and Stephen Havriliak.  Kinetic balance: A partial solution to the problem of variational safety in dirac calculations.  J. Chem. Phys., 81(4):1910–1918, 1984.
  * [175] Kenneth G Dyall.  Relativistic double-zeta, triple-zeta, and quadruple-zeta basis sets for the 4d elements y–cd.  Theor. Chem. Acc., 117(4):483–489, 2007.
  * [176] Kenneth G Dyall and Andre S P Gomes.  Revised relativistic basis sets for the 5d elements hf–hg.  Theor. Chem. Acc., 125(1-2):97–100, 2010.
  * [177] Kenneth G Dyall.  Dyall dz, tz, and qz basis sets for relativistic electronic structure calculations, 2023.  Zenodo archive doi:10.5281/zenodo.7574629 https://zenodo.org/records/7574629.
  * [178] Kenneth G Dyall.  Dyall double-zeta, triple-zeta, and quadruple-zeta basis set archive files, 2023.  Zenodo archive doi:10.5281/zenodo.7606547 https://zenodo.org/records/7606547.
  * [179] Thom H Dunning Jr.  Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen.  J. Chem. Phys., 90(2):1007–1023, 1989.
  * [180] John C Slater.  A simplification of the hartree-fock method.  Phys. Rev., 81(3):385, 1951.
  * [181] S H Vosko, L Wilk, and M Nusair.  Accurate spin-dependent electron liquid correlation energies for local spin density calculations: a critical analysis.  Can. J. Phys., 58(8):1200–1211, 1980.
  * [182] Axel D Becke.  Density-functional exchange-energy approximation with correct asymptotic behavior.  Phys. Rev. A, 38(6):3098–3100, 1988.
  * [183] Chengteh Lee, Weitao Yang, and Robert G Parr.  Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density.  Phys. Rev. B, 37(2):785–789, 1988.
  * [184] P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J. Frisch.  Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields.  J. Phys. Chem., 98(45):11623–11627, 1994.
  * [185] Jun Gao, Wenli Zou, Wenjian Liu, Yunlong Xiao, Daoling Peng, Bo Song, and Chengbu Liu.  Time-dependent four-component relativistic density-functional theory for excitation energies. II. The exchange-correlation kernel.  J. Chem. Phys., 123(5), 2005.
  * [186] Michal Repisky, Lukas Konecny, Marius Kadek, Stanislav Komorovsky, Olga L Malkin, Vladimir G Malkin, and Kenneth Ruud.  Excitation energies from real-time propagation of the four-component Dirac–Kohn–Sham equation.  J. Chem. Theory Comput., 11(3):980–991, 2015.
  * [187] J. E. Sansonetti and W. C. Martin.  Handbook of Basic Atomic Spectroscopic Data.  J. Phys. Chem. Ref. Data, 34(4):1559–2259, 2005.
  * [188] Mark Kamper Svendsen, Michael Ruggenthaler, Hannes Hübener, Christian Schäfer, Martin Eckstein, Angel Rubio, and Simone Latini.  Theory of Quantum Light–Matter Interaction in Cavities: Extended Systems and the Long Wavelength Approximation.  arXiv:2312.17374.
  * [189] L Edwards, DH Dolphin, Martin Gouterman, and Alan D Adler.  Porphyrins XVII. Vapor absorption spectra and redox reactions: Tetraphenylporphins and porphin.  J. Mol. Spectrosc., 38(1):16–32, 1971.
  * [190] EJ Baerends, G Ricciardi, A Rosa, and SJA Van Gisbergen.  A DFT/TDDFT interpretation of the ground and excited states of porphyrin and porphyrazine complexes.  Coord. Chem. Rev., 230(1-2):5–27, 2002.
  * [191] Paul N Day, Kiet A Nguyen, and Ruth Pachter.  Calculation of one-photon and two-photon absorption spectra of porphyrins using time-dependent density functional theory.  J. Chem. Theory Comput., 4(7):1094–1106, 2008.
  * [192] Aleksandr G Avramenko and Aaron S Rury.  Cavity polaritons formed from spatially separated quasi-degenerate porphyrin excitons: Structural modulations of bright and dark state energies and compositions.  J. Phys. Chem. C, 126(37):15776–15787, 2022.
  * [193] Shichao Sun, Bing Gu, and Shaul Mukamel.  Polariton ring currents and circular dichroism of Mg-porphyrin in a chiral cavity.  Chem. Sci., 13(4):1037–1048, 2022.


††thanks<EMAIL_ADDRESS>email:junrenshi@pku.edu.cn(Corresponding Author)

# Density functional theory of Composite Fermions

Zhang Yin-Han  Beijing National Laboratory for Condensed Matter Physics and
Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China  Shi
Jun-Ren  International Center for Quantum Materials, Peking University,
Beijing 100871, China  Collaborative Innovation Center of Quantum Matter,
Beijing,100871, China

###### Abstract

We construct a density functional theory for two-dimension electron (hole)
gases subjected to both strong magnetic fields and external potentials. In
particular, we are focused on regimes near even-denominator filling factors,
in which the systems form composite fermion liquids. Our theory provides a
systematic and rigorous approach to determine the properties of ground states
in fractional quantum Hall regime modified by artificial structures. We also
propose a practical way to construct approximated functional.

###### pacs:

71.10.Pm,73.43.-f,75.75.Cd,71.15.Mb

High mobility semiconductor two-dimensional electron (hole) gas (2DEG/2DHG) is
one of the cleanest condensed matter systems and hosts some of most important
discoveries in physics, including quantum Hall effect Klitzing1980  and
fractional quantum Hall effect Tsui1982 . Recently, there is a surge of
interest in designing and imposing artificial structures in these systems. For
instance, it is proposed that a periodic potential modulation could modify
2DEGs to create graphene-like electronic structure, upon which one could
engineer exotic phases such as topological insulator and superconductivity
Polini2013 . This has become practical as semiconductor nano-fabrication
techniques provide versatile and feasible tools to realize these designs in
laboratory Evers2013 . Until recently, most of studies along the line are
focused on non-interacting systems. On the other hand, a recent theoretical
study reveals that the artificial structure could induce novel and interesting
topological phases in strongly-correlated regime of fractional quantum Hall
effect (FQHE), namely, the quantum anomalous Hall insulator of composite
fermions Zhang2014 . The development asks for a systematic approach to treat
the extremely difficult theoretical problem involving external magnetic
fields, potentials, as well as electron-electron interaction.

Density functional theory provides such a systematic and rigorous approach to
tangle the problem. For the usual solid state systems, the density functional
theory has achieved great successes in predicting structures and properties.
In the presence of the external magnetic field, the theory needs to be
modified and extended to the spin and current-density functional theory
(SCDFT) Vignale1987 ; Vignalebook . In principle, the theory is rigorous and
sufficient to predict ground states in arbitrarily strong magnetic field,
provided one knows _exact_ energy functional in terms of electron/spin density
and paramagnetic current density. In practice, however, the energy functional
is determined by the assumption that the paramagnetic current density is
small, effectively limiting the applicability of the theory in the regime of
small magnetic fields. The more ambitious attempt to predict fractional
quantum Hall states within the density functional theory requires substantial
modification of the theoretical framework and yields limited success
Heinonen1995 . It is clear that to treat our problem, a new approach is needed
to circumvent these difficulties.

In this Letter, we develop a density functional theory for composite fermions
(CF). In the theory, instead of using the electron paramagnetic current
density, we express the energy functional in term of the paramagnetic current
density of _composite fermions_ Jains ; Heinonen ; Kalmeyer1992 ; Willett1993
. It is known that 2DEG/2DHG systems with even-denominator filling factors
form composite fermion liquids subjected to zero effective magnetic field, and
the CF paramagnetic current density vanishes HLR . As a result, we are able to
construct the energy functional in the vicinity of these filling factors by
assuming the smallness of the CF paramagnetic current density. To do this, we
re-formulate the SCDFT for CFs, and present details of derivation for the
theoretical framework sketched in an earlier publication Zhang2014 . We
further determine the gauge invariance of the exchange correlation energy
functional and show a practical way to determine it.

Consider a spinless 2DEG system near a magnetic filling factor $\nu=1/2p$ with
$p$ being an integer under a set of the external scalar field $U(\bm{r})$ and
vector field $\bm{A}(\bm{r})$. The Hamiltonian for electrons can be written as

$\displaystyle\hat{H}_{e}=\hat{T}_{e}+\hat{V}_{ext}+\hat{V}_{ee}$ (1)

where the kinetic energy and the external scalar potential operators are

$\displaystyle\hat{T}_{e}=$
$\displaystyle\sum_{j}\left[\frac{1}{2m_{b}}\left(-i\hbar\nabla_{j}-\frac{e}{c}\bm{A}(\bm{r}_{j})\right)^{2}\right],$
(2) $\displaystyle\hat{V}_{ext}$ $\displaystyle=\sum_{j}U(\bm{r}_{j}),$ (3)

and $\hat{V}_{ee}$ is the electron-electron interaction. For simplicity, we
further decompose the vector potential $\bm{A}$ into
$\bm{A}_{0}+\bm{A}^{\prime}$, where $\bm{A}_{0}$ is the vector potential that
gives rise to a filling factor $\nu=1/2p$, i.e.,
$B_{0}\equiv\mathrm{curl}\bm{A}_{0}=2p\phi_{0}n_{0}$, with $\phi_{0}$ being
the magnetic flux quantum and $n_{0}$ being average electron density of the
system.

We can reformulate the problem by making a Chern-Simons transformation:

$\displaystyle\Psi(\bm{r}_{1},..,\bm{r}_{N})\longrightarrow\exp(-2pi\sum_{i<j}\arg(\bm{r}_{i}-\bm{r}_{j}))\Psi(\bm{r}_{1},...,\bm{r}_{N}),$
(4)

where $\arg(\bm{r})$ is the angle between the vector $\bm{r}$ and the positive
real axis. As a result, the kinetic energy operator is transformed to,

$\displaystyle\hat{T}_{CS}=$
$\displaystyle\sum_{j}\left[\frac{1}{2m_{b}}\left(\hat{\bm{p}}_{j}-\frac{e}{c}(\bm{A}_{0}(\bm{r}_{j})+\bm{A}^{\prime}(\bm{r}_{j})+\bm{a}_{CS}(\bm{r}_{j}))\right)^{2}\right],$
(5)

and the Chern-Simons vector potential is,

$\displaystyle\bm{a}_{CS}(\bm{r}_{j})=\frac{p\phi_{0}}{\pi}\sum_{i\neq
j}\frac{\hat{z}\times(\bm{r}_{j}-\bm{r}_{i})}{|\bm{r}_{j}-\bm{r}_{i}|^{2}}.$
(6)

The resulting Hamiltonian can be solved by a mean field approximation
$\mathrm{curl}\bm{a}_{CS}(\bm{r})=-2p\phi_{0}n(\bm{r})$. With the
approximation and in the uniform density limit, $\bm{a}_{CS}$ in Eq. (5) will
exactly cancel $\bm{A}_{0}$. It gives rise to the very successful CF picture
of the fractional quantum Hall effect, in which each electron attaches $2p$
quantum vortices to form a CF subjected to an effective vector potential
$\bm{A}^{\prime}(\bm{r})$ of reduced effective magnetic strength Jains ; HLR ;
Heinonen ; Kalmeyer1992 ; Willett1993 . In particular, systems with $\nu=1/2p$
is equivalent to a CF liquid with zero effective magnetic field HLR .

Starting with the CF representation, we can develop a density functional
theory to go beyond the usual CF mean field theory Zhang2014 ; HLR ;
Kalmeyer1992 ; Heinonen . We introduce the CF-paramagnetic current density:

$\displaystyle\bm{j}_{p}(\bm{r}_{1})$
$\displaystyle=\int\mathrm{d}\bm{r}_{2}\dots\mathrm{d}\bm{r}_{N}\left[-\frac{i\hbar}{2m_{b}}\left(\Psi^{*}(\bm{\nabla}_{1}\Psi)-(\bm{\nabla}_{1}\Psi^{*})\Psi\right)\right.$
$\displaystyle\left.-\frac{e}{c}\left(A_{0}(\bm{r}_{1})+\frac{p\phi_{0}}{m_{b}\pi}\sum_{i\neq
1}\frac{\hat{z}\times(\bm{r}_{1}-\bm{r}_{i})}{|\bm{r}_{1}-\bm{r}_{i}|^{2}}\right)|\Psi|^{2}\right]$
(7)

It is important to note that $\bm{j}_{p}(\bm{r})$ vanishes for a uniform
2DEG/2DHG at $\nu=1/2p$.

We can prove generalized Hohenberg-Kohn (HK) theorem which dictates that
ground state energy is a functional of density $n(\bm{r})$ and CF-paramagnetic
current density $\bm{j}_{p}(\bm{r})$. First, we prove that the ground-state
distributions of $n(\bm{r})$ and $\bm{j}_{p}(\bm{r})$ uniquely determine the
external scale potential $U(\bm{r})$ and vector potential
$\bm{A}^{\prime}(\bm{r})$, and hence the non-degenerate ground state
$|\Psi\rangle$. To see this, we assume that there exist two sets of external
fields $U_{1(2)}(\bm{r})$, $\bm{A}^{\prime}_{1(2)}(\bm{r})$ that yield the
same set of distributions $n(\bm{r})$ and $\bm{j}_{p}(\bm{r})$. The CF
Hamiltonian $\hat{H}_{1(2)}$ have ground-state $|\Psi_{1(2)}\rangle$ with
their ground-state energy $E_{1(2)}$. For $|\Psi_{1}\rangle$ and
$\hat{H}_{1}$, we have

$\displaystyle E_{1}<\langle\Psi_{2}|\hat{H}_{1}|\Psi_{2}\rangle=E_{2}+\int
d^{2}r(U_{2}(\bm{r})-U_{1}(\bm{r}))n(\bm{r})$ $\displaystyle-\frac{e}{c}\int
d^{2}r(\bm{A}_{1}(\bm{r})-\bm{A}_{2}(\bm{r}))\cdot\bm{j}_{p}(\bm{r})$
$\displaystyle+\frac{e^{2}}{2m_{b}c^{2}}\int
d^{2}rn(\bm{r})\left[\bm{A}_{1}^{2}(\bm{r})-\bm{A}_{2}^{2}(\bm{r})\right].$

Similarly, for $|\Psi_{2}\rangle$ and $\hat{H}_{2}$, another inequality is
obtained by exchanging the suffix $(1\leftrightarrow 2)$. Summing up the two
inequalities, one observes contradiction $E_{1}+E_{2}<E_{2}+E_{1}$.

Since the non-degenerate ground state wave function $|\Psi\rangle$ is a
functional of the density distributions $n(\bm{r})$ and $\bm{j}_{p}(\bm{r})$,
the energy expectation value can also be expressed as a functional:

$\displaystyle E_{U,\bm{A}^{\prime}}[n,\bm{j}_{p}]=F[n,\bm{j}_{p}]+\int
d^{2}rn(\bm{r})U(\bm{r})$ $\displaystyle-\frac{e}{c}\int
d^{2}r\bm{A}^{\prime}(\bm{r})\cdot\bm{j}_{p}(\bm{r})+\frac{e^{2}}{2m_{b}c^{2}}\int
d^{2}r\bm{A}^{\prime 2}(\bm{r})n(\bm{r})$ (8)

with a universal functional defined as

$F[n,\bm{j}_{p}]\equiv\min_{\Psi^{\prime}\rightarrow\\{n,\bm{j}_{p}\\}}\langle\Psi^{\prime}|\hat{T}_{0}+\hat{V}_{ee}|\Psi^{\prime}\rangle$
(9)

where $\hat{T}_{0}\equiv\left.\hat{T}_{CS}\right|_{\bm{A}^{\prime}=0}$, and
the search for the minimum is restricted in all possible asymmetry wave
function $\Psi^{\prime}$ that yield the give density distributions. With the
energy functional, the ground state density distributions can be obtained by
minimizing the functional against $n(\bm{r})$ and $\bm{j}_{p}(\bm{r})$
Vignalebook .

We can further map the many-body problem to an effective non-interacting one.
To do this, we assume that for any ground state density distribution yielded
by the interacting system, one can find a non-interacting reference system
that has the same set of ground-state density distributions. The non-
interacting system is obtained by solving Kohn-Sham equation:

$\displaystyle\left[\frac{1}{2m_{b}}(-i\hbar\nabla-\frac{e}{c}\bm{A}_{KS}(\bm{r}))^{2}+U_{KS}(\bm{r})\right]\psi_{j}(\bm{r})=\epsilon_{j}\psi_{j}(\bm{r})$
(10)

And the corresponding density distributions are determined by,

$\displaystyle\bm{j}_{p}(\bm{r})$
$\displaystyle=\sum^{N}_{j=1}\frac{\hbar}{2m_{b}i}\left(\psi^{\ast}_{j}(\bm{r})\nabla\psi_{j}(\bm{r})-c.c\right)$
$\displaystyle-\frac{e}{m_{b}c}\left(\bm{A}_{0}(\bm{r})+\bar{\bm{a}}_{CS}(\bm{r})\right)n(\bm{r}),$
$\displaystyle n(\bm{r})=\sum_{i=1}^{N}|\psi_{i}|^{2}.$ (11)

where the summation is for the $N$-lowest energy states, and

$\bar{\bm{a}}_{CS}(\bm{r})=\frac{p\phi_{0}}{\pi}\int\mathrm{d}\bm{r}^{\prime}\frac{\hat{z}\times(\bm{r}-\bm{r}^{\prime})}{|\bm{r}-\bm{r}^{\prime}|^{2}}n(\bm{r}^{\prime}),$
(12)

is the mean field expectation value of $\bm{a}_{CS}$. In Eq. (10), we
introduce Kohn-Sham potentials $U_{KS}$ and $\bm{A}_{KS}$, which are
functionals of the density distributions.

For the non-interacting system, there is also a HK theorem dictating the
existence of a universal energy functional:

$\displaystyle
T_{s}=\min_{\Psi_{s}\rightarrow\\{n,\bm{j}_{p}\\}}\langle\Psi_{s}|\hat{T}_{0}|\Psi_{s}\rangle,$
(13)

where the search for the minimum is restricted in all possible non-interacting
wave functions of slater determinants. We can decompose the universal energy
functional of the interacting system $F[n,\bm{j}_{p}]$ into:

$\displaystyle F[n,\bm{j}_{p}]=T_{s}[n,\bm{j}_{p}]+\frac{1}{2}\int
d^{2}r_{1}r_{2}$ $\displaystyle n(\bm{r}_{1})V_{ee}(\bm{r})n(\bm{r}_{2})$
$\displaystyle+E_{xc}[n,\bm{j}_{p}],$ (14)

where we introduce exchange-correlation energy functional $E_{xc}$ to register
residual energy correction.

Using definition of $T_{s}$, one can re-express it as:

$\displaystyle T_{s}[n,\bm{j}_{p}]=\sum_{i}\epsilon_{i}-\frac{e}{c}\int
d^{2}r\left[\bm{A}_{0}(\bm{r})+\bm{a}_{CS}(\bm{r})-\bm{A}_{KS}(\bm{r})\right]$
$\displaystyle\cdot\left[\bm{j}_{p}(\bm{r})+\frac{e}{m_{b}c}n(\bm{r})(\bm{A}_{0}(\bm{r})+\bm{a}_{CS}(\bm{r}))\right]-n(\bm{r})U_{KS}(\bm{r})$
$\displaystyle+\int
d^{2}rn(\bm{r})\left[\frac{e^{2}}{2m_{b}c^{2}}\left((\bm{A}_{0}(\bm{r})+\bm{a}_{CS}(\bm{r}))^{2}-\bm{A}^{2}_{KS}(\bm{r})\right)\right],$
(15)

with $\epsilon_{j}$, $\bm{a}_{CS}(\bm{r})$, and KS potentials $U_{KS}$ and
$\bm{A}_{KS}$ being functionals of $n$ and $\bm{j}_{p}$.

Substituting Eq. (15) into Eq. (14) and then Eq. (8), we obtain the total
energy functional $E_{U,\bm{A}^{\prime}}[n,\bm{j}_{p}]$. To determine KS
potentials, we carry out minimization of energy functional
$E_{U,\bm{A}^{\prime}}[n,\bm{j}_{p}]$ with respect of $n(\bm{r})$ and
$\bm{j}_{p}(\bm{r})$:

$\displaystyle\frac{\delta E_{U,\bm{A}^{\prime}}[n,\bm{j}_{p}]}{\delta
n(\bm{r})}=0,\quad\frac{\delta
E_{U,\bm{A}^{\prime}}[n,\bm{j}_{p}]}{\delta\bm{j}_{p}(\bm{r})}=0.$ (16)

The KS potentials can then be determined:

$\displaystyle\bm{A}_{KS}(\bm{r})=\bm{A}(\bm{r})+\bm{a}_{CS}(\bm{r})+\bm{A}_{xc}(\bm{r}),$
$\displaystyle U_{KS}(\bm{r})=U(\bm{r})+U_{xc}(\bm{r})+\int
dr^{\prime}V_{ee}(\bm{r},\bm{r}^{\prime})n(\bm{r}^{\prime})$
$\displaystyle+\frac{e^{2}}{2m_{b}c^{2}}\left(\bm{A}^{\prime
2}(\bm{r})-(\bm{A}^{\prime}(\bm{r})+\bm{A}_{xc}(\bm{r}))^{2}\right)+2p\phi_{0}m_{z}(\bm{r}),$
(17)

where $m_{z}(\bm{r})$ is the orbital magnetization density of the non-
interacting KS system Shi2007 , and the exchange-correlation corrections to
the KS potentials are:

$\displaystyle U_{xc}(\bm{r})\equiv\frac{\delta E_{xc}[n,\bm{j}_{p}]}{\delta
n(\bm{r})},\quad\bm{A}_{xc}(\bm{r})\equiv-\frac{c}{e}\frac{\delta
E_{xc}[n,\bm{j}_{p}]}{\delta\bm{j}_{p}(\bm{r})}.$ (18)

We can show that gauge invariance puts important constraint on the possible
form of exchange-correlation energy functional. We introduce the gauge
transformation:

$\displaystyle\bm{A}(\bm{r})\rightarrow\bm{A}_{G}=\bm{A}(\bm{r})+\nabla
f(\bm{r}),$ (19)

with an arbitrary analytic function $f(\bm{r})$. The density distributions $n$
and $\bm{j}_{p}$, and the wave function $\Psi$ of the system are transformed
as

$\displaystyle
n_{G}(\bm{r})=n(\bm{r}),\quad\Psi_{G}=\exp(\sum_{j}\frac{ief(\bm{r}_{j})}{\hbar
c})\Psi,$
$\displaystyle\bm{j}_{Gp}(\bm{r})=\bm{j}_{p}(\bm{r})+\frac{e}{m_{b}c}\nabla
f(\bm{r})n(\bm{r}).$ (20)

The invariance of the energy functional under the gauge transformation
dictates:

$\displaystyle F$ $\displaystyle[n,\bm{j}_{p}+\frac{e}{m_{b}c}n\nabla
f]=F[n,\bm{j}_{p}]$ $\displaystyle+\frac{e^{2}}{2m_{b}c^{2}}\int
d^{2}rn(\bm{r})(\nabla f(\bm{r}))^{2}+\frac{e}{c}\int d^{2}r\nabla
f(\bm{r})\cdot\bm{j}_{p}(\bm{r}).$ (21)

The kinetic functional $T_{s}$ for the non-interacting KS system transforms
similarly. It follows that the exchange-correlation functional is gauge
invariant:

$\displaystyle
E_{xc}[n(\bm{r}),\bm{j}_{p}(\bm{r})+\frac{e}{m_{b}c}n(\bm{r})\nabla
f(\bm{r})]=E_{xc}[n(\bm{r}),\bm{j}_{p}(\bm{r})].$ (22)

As a result, the exchange-correlation energy functional can be rewritten as,

$\displaystyle E_{xc}[n,\bm{j}_{p}]\equiv\bar{E}_{xc}[n,\bm{v}],$ (23)

which only depends on vorticity defined as:

$\displaystyle\bm{v}(\bm{r})\equiv\nabla\times\frac{\bm{j}_{p}(\bm{r})}{n(\bm{r})}.$
(24)

The exchange-correlation scalar and vector potential can be rewritten as

$\displaystyle U_{xc}(\bm{r})=\frac{\delta\bar{E}_{xc}[n,\bm{\nu}]}{\delta
n(\bm{r})}+$
$\displaystyle\frac{e}{c}\bm{A}_{xc}(\bm{r})\cdot\frac{\bm{j}_{p}(\bm{r})}{n(\bm{r})},$
$\displaystyle-\frac{e}{c}\bm{A}_{xc}(\bm{r})=-\frac{1}{n(\bm{r})}$
$\displaystyle\nabla\times\frac{\delta\bar{E}_{xc}[n,\bm{\nu}]}{\delta\bm{\nu}(\bm{r})}.$
(25)

To determine $E_{xc}$, we employ the local density approximation (LDA) by
assuming that $E_{xc}$ at position $\bm{r}$ is only determined by the local
density $n(\bm{r})$. We further exploit the fact that the CF paramagnetic
current density vanishes at the uniform limit and at $\nu=1/2p$. It suggests
that $\bm{j}_{p}(\bm{r})$ should be small in the vicinity of even-denominator
filling factors. Based on the stiffness theorem Vignalebook , the universal
functional $F$ can be approximated as, in the uniform density limit:

$\displaystyle
F[n,\bm{j}_{p}]=F[n,0]-\frac{1}{2V}\sum_{\bm{q}}\delta\bm{j}_{p}(-\bm{q})\chi^{-1}_{\bm{j}_{p},\bm{j}_{p}}(\bm{q})\delta\bm{j}_{p}(\bm{q})$
(26)

where the static paramagnetic current response function matrix
$\chi_{\bm{j}_{p},\bm{j}_{p}}(\bm{q})$ is defined as

$\displaystyle\chi_{\bm{j}_{p},\bm{j}_{p}}(\bm{q})=-i\int
dt\Theta(t)\langle[\bm{j}_{p}(\bm{q},t),\bm{j}_{p}(-\bm{q},0)]\rangle$ (27)

For an isotropic system, we can decompose

$\displaystyle\chi_{\bm{j}_{p},\bm{j}_{p}}(\bm{q})=\chi_{\parallel}(\bm{q})\frac{q_{i}q_{j}}{q^{2}}+\chi_{\perp}(\bm{q})(\delta_{i,j}-\frac{q_{i}q_{j}}{q^{2}})$
(28)

into a longitudinal component $\chi_{\parallel}(\bm{q})$ and a transverse
component $\chi_{\perp}(\bm{q})$. Because a purely longitudinal and static
vector potential cannot induce any physical current, i.e.,
$\chi_{\parallel}(\bm{q},0)=-n/m$, we can consider transverse current response
$\chi_{\perp}(\bm{q})$ only.

It is straightforward to determine the exchange-correlation energy functional:

$\displaystyle E_{xc}[n,\bm{j}_{p}]=$
$\displaystyle-\frac{1}{2V}\sum_{\bm{q}}\delta\bm{j}^{\perp}_{p}(-\bm{q})(\chi^{-1}_{\perp}(\bm{q})-\chi^{(0)-1}_{\perp}(\bm{q}))\delta\bm{j}^{\perp}_{p}(\bm{q})$
$\displaystyle+E_{xc}[n,0],$ (29)

where $\chi_{\perp}^{0}$ is the transverse response function for the non-
interacting system.

The exchange-correlation energy functional has the general form:

$\displaystyle\bar{E}_{xc}[n,\bm{v}]=E_{xc}[n,0]+\frac{1}{2}\int
dr^{2}_{1}dr^{2}_{2}f(\bm{r}_{1}-\bm{r}_{2})\bm{v}(\bm{r}_{1})\cdot\bm{v}(\bm{r}_{2})$
(30)

where the function $f(\bm{r})$ is related to the transverse current response
function by,

$\displaystyle
f(\bm{r})=\frac{1}{V}\sum_{\bm{q}}\frac{n^{2}}{q^{2}}(\chi^{-1}_{\perp}(\bm{q})-\chi^{(0)-1}_{\perp}(\bm{q}))\exp(i\bm{q}\cdot\bm{r}).$

In general, $f(\bm{r})$ is a function of the density $n$, which can be
replaced with average density $n_{0}$, as the first order of approximation.

We have formulated the current density functional theory for CFs. A gauge
invariant self-consistent Kohn-Sham equation for CFs is constructed. The
general form of the exchange-correlation functional is established. Explicit
determination of the exchange-correlation functional requires numerical
evaluation of ground state energy of uniform 2DEG at $\nu=1/2p$ and
corresponding transverse current response function, which will be investigated
in the future publications.

This work is supported by 973 program of China (2012CB921304) and National
Science Foundation of China (11325416).

## References

  * (1) Klitzing K V, Dorda G and Pepper M 1980 Phys. Rev. Lett. 45 494
  * (2) Tsui D C, Stormer H L and Gossard A C 1982 Phys. Rev. Lett. 48 1559
  * (3) Polini M, Guinea F, Lewenstein M, Manoharan H C and Pellegrini V 2013 Nature Nanotech. 4 625
  * (4) Evers W H, Goris B, Bals S, Casavola M, Graaf J D, Roij R V, Dijkstra M, and Vanmaekelbergh D 2013 Nano Lett. 13, 2317
  * (5) Zhang Y and Shi J 2014 Phys. Rev. Lett. 113, 010806
  * (6) Vignale G and Rasolt M 1987 Phys. Rev. Lett. 59, 2360; 1988 Phys. Rev. B 37, 10685
  * (7) Giuliani G F and Vignale G 2005 Quantum Theory of the Electron Liquid (Cambridge: Cambridge University Press) chap 7 p 390
  * (8) Heinonen O, Lubin M I and Johnson M J 1995 Phys. Rev. Lett. 75, 4110
  * (9) Jain J K 2007 Composite Fermions (Cambridge: Cambridge University Press) chap 5 p105
  * (10) Heinonen O 1998 Composite Fermions (Singapore: World Scientific)
  * (11) Kalmeyer V and Zhang S C 1992 Phys. Rev. B 46, 9889
  * (12) Willett R L, Ruel R R, West K W and Pfeiffer L N 1993 Phys. Rev. Lett. 71, 3846
  * (13) Halperin B I, Lee P A and Read N 1993 Phys. Rev. B 47, 7312
  * (14) Shi J, Vignale G, Xiao D and Niu Q 2007 Phys. Rev. Lett. 99, 197202


#  A Data-Driven Modeling Framework of Time-Dependent Switched Dynamical
Systems via Extreme Learning Machine

Weiming Xiang Weiming Xiang is with the School of Computer and Cyber Sciences,
Augusta University, Augusta GA 30912 USA. Email<EMAIL_ADDRESS>
###### Abstract

In this work, a data-driven modeling framework of switched dynamical systems
under time-dependent switching is proposed. The learning technique utilized to
model system dynamics is Extreme Learning Machine (ELM). First, a method is
developed for the detection of the switching occurrence events in the training
data extracted from system traces. The training data thus can be segmented by
the detected switching instants. Then, ELM is used to learn the system
dynamics of subsystems. The learning process includes segmented trace data
merging and subsystem dynamics modeling. Due to the specific learning
structure of ELM, the modeling process is formulated as an iterative Least-
Squares (LS) optimization problem. Finally, the switching sequence can be
reconstructed based on the switching detection and segmented trace merging
results. An example of the data-driven modeling DC-DC converter is presented
to show the effectiveness of the developed approach.

##  I Introduction

A switched system is composed of a family of continuous or discrete-time
subsystems along with a switching rule orchestrating the activation of them.
Switched systems are frequently used to model and characterize the behaviors
of practical systems that are inherently multi-modal such as power electrical
systems, medical robotic systems, autonomous vehicles, and a variety of Cyber-
Physical Systems (CPS). Relevant analysis, control, and verification methods
and tools for hybrid systems were developed based on explicit models. For
instance, multiple Lyapunov function approaches were developed for stability
analysis for switched systems and extended to a variety of related problems
[1, 2, 3, 4]. After several decades of development of switched and hybrid
system theory and increasingly implementations to various applications,
switched system theory shows its powerful capability in dealing with complex
dynamical systems, and those achievements are made on the basis of explicit
mathematical models. However, accurate models of many practical systems such
as complex CPS are often not easy (or even impossible) to obtain due to
incomplete knowledge of certain processes.

Data-driven methods have been shown to outperform physics-based models in many
disciplines, especially for those physical processes that are not fully
understood by researchers, but for which data of adequate quality and quantity
is available. When applying data-driven modeling methods for switched systems,
there are three unique challenging problems that need to be taken into
account:

  * •

_How does one detect the switching occurrence and segment collected system
traces/data?_

  * •

_How does one classify and model multiple subsystems?_

  * •

_How does one reconstruct switching signals?_

The above three core problems are the main focus of data-driven modeling of
switched systems. In [5], a general convex framework for the identification of
switched linear systems is proposed in the framework of the Least-Squares (LS)
problem with sum-of-norms regularization. A model identification framework is
proposed for hybrid systems consisting of linear subsystems in [6], in which
linear regression is used to estimate the dynamics. In [7], the switched
system identification problem is formulated as a constrained multi-objective
optimization problem. A recursive procedure is developed in [8] for the
identification of switched linear models from input-output data. For more
details of the results on the modeling and identification for switched systems
and hybrid systems, the reader is referred to survey papers such as [9, 10],
and the references cited therein. This paper aims to propose a data-driven
modeling framework of time-dependent switched systems by using a learning
technique called Extreme Learning Machine (ELM) which is with a learning
structure appropriate for modeling switched systems. First, a switching
detection method is developed to detect the occurrence of switchings and
segment the collected system traces. Then, taking advantage of the ELM
learning structure in the form of an LS optimization problem, an algorithm is
developed to merge segmented traces and then model subsystem dynamics.
Finally, the switching sequence can be reconstructed based on the switching
detection and segmented trace merging results.

The rest of this paper is organized as follows. The problem formulation and
preliminaries are given in Section II. A data-driven modeling framework of
switched systems using ELM are presented in Section III. In Section IV, an
evaluation is performed on the modeling of DC-DC converters. Conclusion and
further remarks are given in Section V.

##  II Preliminaries and Problem Formulation

In this paper, the dynamical systems considered for modeling are assumed to be
a class of switched nonlinear systems in the following form

$\mathbf{x}(k+1)=f_{\theta(k)}\left(\mathbf{x}(k),\mathbf{u}(k)\right)$ (1)

where $\mathbf{x}(k)\in\mathbb{R}^{n_{x}}$ is the system state, and
$\mathbf{u}(k)\in\mathbb{R}^{n_{u}}$ is the system input, respectively.
Defining index set $\mathcal{I}\triangleq\\{1,2,\ldots,N\\}$ where $N$ is the
number of subsystems, $\theta(k):\mathbb{N}\to\mathcal{I}$ is a function of
time $k$ denoting the switching signal. $f_{i}$, $i\in\mathcal{I}$ is the
unknown nonlinear discrete-time process for the $i$th subsystem,
characterizing system dynamical behaviors. However, accurate models of many
practical systems are often not easy (or even impossible) to obtain due to
incomplete knowledge of certain processes including switched nonlinear systems
with multiple subsystems along with a switching law. In this paper, we are
going to develop a data-driven modeling framework to construct individual
subsystem dynamics $f_{i}$, $i\in\mathcal{I}$ and switching law $\theta(k)$
using the data collected from system trajectories and inputs, i.e., the input-
state traces $\\{\mathbf{x}(k),\mathbf{u}(k)\\}$, $k=0,1,\ldots,K$.

###### Definition 1

A input-state trace $\phi_{K}$ is a sequence of pairs of
$\psi(k)=\\{\mathbf{x}(k),\mathbf{u}(k)\\}$, $k=0,1,\ldots,K$ generated by
system (1) starting from initial values
$\psi(0)=\\{\mathbf{x}(0),\mathbf{u}(0)\\}$.

In this paper, it is assumed that there exist sufficient measurable input-
state traces available for the data-driven modeling for switched nonlinear
systems in the form of (1).

###### Assumption 1

It is assumed that state $\mathbf{x}(k)$ and input $\mathbf{u}(k)$ are all
measurable for data-driven modeling of switched nonlinear system (1).

###### Remark 1

In practice, the collected data inevitably contains noise, e.g., measured
state $\hat{\mathbf{x}}(k)=\mathbf{x}(k)+\delta(k)$, where $\delta(k)$ is the
noise vector at instant $k$. There are various noise reduction algorithms
available to reduce or remove noise from a signal. These denosing methods
heavily rely on the types of signals and noise, and multiple methods are
successful within the scope of their target signal and noise type. In the rest
of the paper, we assume that collected traces are well pre-denoised using some
noise reduction method that is not in the scope of this work, and we still use
$\mathbf{x}(k)$ and $\mathbf{u}(k)$ to denote denoised measured state and
input, if there is no ambiguity in the context.

Extreme Learning Machine (ELM) is dedicated to train Single-hidden Layer
Feedforward Networks (SLFN). An SLFN with $L$ hidden nodes can be represented
in the form of

$\displaystyle\mathbf{y}=f(\mathbf{x})=\sum\nolimits_{i=1}^{L}\boldsymbol{\upbeta}_{i}h_{i}(\mathbf{x})=\mathbf{h}(\mathbf{x})\boldsymbol{\upbeta}$
(2)

where $\mathbf{x}\in\mathbb{R}^{n}$ and $\mathbf{y}\in\mathbb{R}^{m}$ are the
input and output vectors, respectively. The output nodes are chosen linear and
$\boldsymbol{\upbeta}\triangleq[\boldsymbol{\upbeta}_{1},\ldots,\boldsymbol{\upbeta}_{L}]^{\top}$
is the output weight vector between the hidden layer of $L$ nodes to the
$m\geq 1$ output nodes.
$\mathbf{h}(\mathbf{x})\triangleq[h_{1}(\mathbf{x}),\ldots,h_{L}(\mathbf{x})]$
is the ELM nonlinear feature mapping. In particular, $h_{i}(\mathbf{x})$ can
be defined as follows:

$\displaystyle
h_{i}(\mathbf{x})=g(\mathbf{a}_{i}\cdot\mathbf{x}+b_{i}),~{}\mathbf{a}_{i}\in\mathbb{R}^{d},~{}b_{i}\in\mathbb{R}$
(3)

where $g(\cdot)$ is the activation function for hidden layer neurons including
the sigmoidal functions as well as the radial basis, sine, exponential, and
many other nonregular functions. $\mathbf{a}_{i}$ and $b_{i}$ are the the
input weight vector and the bias weight of the $i$th neuron in hidden layer.

Given $N$ arbitrary distinct input-output samples
$\\{\mathbf{x}_{i},\mathbf{t}_{i}\\}$ with $\mathbf{x}_{i}\in\mathbb{R}^{n}$
and $\mathbf{t}_{i}\in\mathbb{R}^{m}$, ELM aims to find the following minimum
norm Least-Squares (LS) solution of SLFN, which is expressed by

$\displaystyle\min_{\boldsymbol{\upbeta}\in\mathbb{R}^{L\times
m}}\left\|\mathbf{H}\boldsymbol{\upbeta}-\mathbf{T}\right\|$ (4)

where $\mathbf{H}$ is the randomized hidden layer output matrix as

$\displaystyle\mathbf{H}=\begin{bmatrix}\mathbf{h}(\mathbf{x}_{1})\\\
\vdots\\\
\mathbf{h}(\mathbf{x}_{N})\end{bmatrix}=\begin{bmatrix}h_{1}(\mathbf{x}_{1})&\cdots&h_{L}(\mathbf{x}_{N})\\\
\vdots&\ddots&\vdots\\\
h_{1}(\mathbf{x}_{1})&\cdots&h_{L}(\mathbf{x}_{N})\end{bmatrix}$ (5)

and $\mathbf{T}$ is the training data target matrix

$\displaystyle\mathbf{T}=\begin{bmatrix}\mathbf{t}_{1}^{\top}\\\ \vdots\\\
\mathbf{t}_{N}^{\top}\end{bmatrix}=\begin{bmatrix}t_{11}&\cdots&t_{1m}\\\
\vdots&\ddots&\vdots\\\ t_{N1}&\cdots&t_{Nm}\end{bmatrix}.$ (6)

###### Remark 2

Compared with other conventional iterative learning frameworks, this ELM
algorithm can provide a good generalization performance as well as an
extremely fast learning speed. Even with one single randomized hidden layer,
ELM is still a learner with universal approximation and classification
capabilities as proved in [11]. This universal approximation feature make ELM
be able to capture the nonlinear dynamics in switched systems.

###### Remark 3

According to [12, 13], the the smallest norm LS solution
$\boldsymbol{\upbeta}^{*}$ for optimization problem (4) is

$\displaystyle{\boldsymbol{\upbeta}^{*}}=\mathbf{H}^{\dagger}\mathbf{T}$ (7)

where $\mathbf{H}^{\dagger}$ is the Moore-Penrose generalized inverse of
matrix $\mathbf{H}$. Remarkably, the input weight and bias for hidden layer do
not need to be tuned once random values have been assigned at the beginning of
learning, and the learning process only requires the computation on
$\boldsymbol{\upbeta}^{*}$, which is extremely efficiently described by (7).

Under Assumption 1 and assuming that all the collected traces are with a same
initial time and a same length $K$, we assume that there exist $N$ input-state
traces defined by

$\phi_{K}^{(n)}=\\{\psi(0),\psi(1),\ldots,\psi(K)\\},~{}n=1,2,\ldots,N$ (8)

which are collected from the system. The main objective is to extract switched
system model (1) including subsystem dynamics $f_{i}$, $i\in\mathcal{I}$ and
switching law $\theta(k)$. In the framework of ELM, the subsystem dynamics
$f_{i}$, $i\in\mathcal{I}$ are supposed to be modeled in the form of (2). To
achieve the data-driven modeling goal of switched nonlinear system (1), three
core sub-problems need to be solved, i.e., switching instant detection,
subsystem modeling, and switching law reconstruction. Thus, the data-driven
modeling problem of switched systems using ELM is summarized as below:

###### Problem 1

Given $N$ input-state traces $\phi_{K}^{(n)}$, $n=1,2,\ldots,N$ collected from
the system, how does one extract the switched system model (1), including the
following three sub-problems:

  1. 1.

How does one determine the set of switching instants
$\mathcal{S}=\\{k_{1},k_{2},\ldots,k_{\ell}\\}$, $0<k_{1}\leq\ k_{\ell}<K$?

  2. 2.

How does one model subsystem dynamics $f_{i}$, $i\in\mathcal{I}$ by ELM?

  3. 3.

How does one reconstruct time-dependent switching law $\theta(k)$?

Problem 1 with the three core problems is the main concern to be addressed in
the rest of this paper.

##  III Data-Driven Modeling Framework

In this section, a data-driven modeling framework of switched systems will be
presented, including switching detection, subsystem modeling, and switching
reconstruction.

###  III-A Switching Detection

The occurrence of switching is usually reflected by the abrupt change of state
trajectories. First, let’s consider a continuous-time switched nonlinear
system

$\dot{\mathbf{x}}(t)=f_{\theta(t)}(\mathbf{x}(t),\mathbf{u}(t))$ (9)

where $f_{i}$, $i\in\mathcal{I}$ are $\mathcal{C}^{\infty}$ functions. The set
of switching instants is denoted by $\mathcal{S}=\\{t_{1},t_{2},\ldots\\}$. To
detect switching instants in $\mathcal{S}$, the key is to identify the events
that indicate the occurrence of switchings.

###### Definition 2

A switching occurrence event is an outcome that can be only resulted from the
occurrence of switching.

Under the trivial assumption that input signals $\mathbf{u}(t)$ are assumed to
be $\mathcal{C}^{0}$ functions and due to the continuity of $f_{i}$, the
following condition explicitly implies an switching occurrence event at time
instant $t$:

$\displaystyle\lim\nolimits_{t\to
t^{-}}f_{\theta(t)}(\mathbf{x}(t),\mathbf{u}(t))\neq\lim\nolimits_{t\to
t^{+}}f_{\theta(t)}(\mathbf{x}(t),\mathbf{u}(t))$ (10)

which can be reflected by the measured state $\mathbf{x}(t)$ such that

$\displaystyle\left\|\mathbf{x}^{\prime}_{+}(t)-\mathbf{x}^{\prime}_{-}(t)\right\|>0$
(11)

where $\mathbf{x}^{\prime}_{-}(t)\ =\limsup\nolimits_{h\to
0^{+}}\frac{\mathbf{x}(t)-\mathbf{x}(t-h)}{h}$ and
$\mathbf{x}^{\prime}_{+}(t)=\limsup\nolimits_{h\to
0^{+}}\frac{\mathbf{x}(t+h)-\mathbf{x}(t)}{h}$. Thus, the switching set
determined by switching occurrence event (11) is
$\mathcal{S}^{\prime}=\\{t\in\mathbb{R}_{+}\mid\left\|\mathbf{x}^{\prime}_{+}(t)-\mathbf{x}^{\prime}_{-}(t)\right\|>0\\}$,
and we have $\mathcal{S}^{\prime}\subseteq\mathcal{S}$.

It worth mentioning that (10) and (11) may not be able to characterize all
switching occurrence events, e.g., if we consider a switched system with two
modes $f_{i}(\mathbf{x}(t))$, $i=1,2$, and the event (11) does not happen if
the system switches at a particular state in set
$\mathcal{X}=\\{\mathbf{x}\mid f_{1}(\mathbf{x})=f_{2}(\mathbf{x})\\}$. In
this case, we can further exploit the continuity of the derivative of $f_{i}$
in addition to $f_{i}$. Given a trivial assumption
$\mathbf{u}(t)\in\mathcal{C}^{1}$ and due to $f_{i}\in\mathcal{C}^{\infty}$,
the following condition implies an switching occurrence event at time instant
$t$:

$\displaystyle\lim\nolimits_{t\to t^{-}}\frac{\partial
f_{\theta(t)}}{\partial\mathbf{x}}$ $\displaystyle
f_{\theta(t)}+\frac{\partial
f_{\theta(t)}}{\partial\mathbf{u}}\frac{\partial\mathbf{u}}{\partial t}$
$\displaystyle\neq\lim\nolimits_{t\to t^{+}}\frac{\partial
f_{\theta(t)}}{\partial\mathbf{x}}f_{\theta(t)}+\frac{\partial
f_{\theta(t)}}{\partial\mathbf{u}}\frac{\partial\mathbf{u}}{\partial t}$ (12)

The above condition can be reflected by

$\displaystyle\left\|\mathbf{x}^{\prime\prime}_{+}(t)-\mathbf{x}^{\prime\prime}_{-}(t)\right\|>0$
(13)

which generates
$\mathcal{S}^{\prime\prime}=\\{t\in\mathbb{R}_{+}\mid\left\|\mathbf{x}^{\prime\prime}_{+}(t)-\mathbf{x}^{\prime\prime}_{-}(t)\right\|>0\\}$,
and the detected switching set is
$\mathcal{S}^{\prime}\bigcup\mathcal{S}^{\prime\prime}\subseteq\mathcal{S}$.

The above process can be further generalized to $p$th derivative of
$\mathbf{x}$ to generate switching set
$\mathcal{S}^{(p)}=\\{t\in\mathbb{R}_{+}\mid\left\|\mathbf{x}^{(p)}_{+}(t)-\mathbf{x}^{(p)}_{-}(t)\right\|>0\\}$,
and the detected switching set is
$\mathcal{S}=\bigcup\nolimits_{i=1}^{p}\mathcal{S}^{(i)}$.

In summary, switching occurrence events can be represented by the
discontinuity of derivatives of system state $\mathbf{x}(t)$, thus the
switching detection problem is converted to the problem of the discontinuity
detection for derivatives of state $\mathbf{x}(t)$. The discontinuity
detection of derivatives of a signal can be performed by multiple methods
within the scope of signal processing, e.g., the well-known wavelet method.

In the discrete-time case, the above discontinuity detection of derivatives of
state $\mathbf{x}(t)$ can be formulated in a discretized version using
$\mathbf{x}(k)$. For the $p$th derivative case, the switching occurs at
instant $k$ if the following condition is triggered

$\displaystyle\mathbf{s}^{(p)}(k)=\frac{\left\|\Delta^{(p)}\mathbf{x}(k+1)-\Delta^{(p)}\mathbf{x}(k)\right\|}{\left\|\Delta^{(p-1)}\mathbf{x}(k)\right\|}\geq\epsilon_{p}$
(14)

where $\Delta^{(0)}\mathbf{x}(k)=\mathbf{x}(k)$ and
$\Delta^{(p)}\mathbf{x}(k)=\Delta^{(p-1)}\mathbf{x}(k)-\Delta^{(p-1)}\mathbf{x}(k-1)$,
and $\epsilon_{p}>0$ is a prescribed threshold. For instance, (11) can be
rewritten as

$\displaystyle\mathbf{s}^{(1)}(k)=\frac{\left\|\mathbf{x}(k+1)-2\mathbf{x}(k)+\mathbf{x}(k-1)\right\|}{\left\|\mathbf{x}(k)\right\|}\geq\epsilon_{1}$
(15)

where $k=1,2\ldots K-1$.

Based on (14), the induced switching set derived from input-state traces
$\phi_{K}^{(n)}$ can be determined by
$\mathcal{S}_{n}=\bigcup\nolimits_{i=1}^{p}\mathcal{S}^{(i)}$ where
$\mathcal{S}^{(i)}=\\{k\mid\mathbf{s}^{(i)}(k)\geq\epsilon_{i},~{}\epsilon_{i}>0\\}$.

###### Remark 4

The proposed condition (14) provides a data-driven detection framework of
switching occurrence events. There are two factors in the proposed framework
that need some further trade-off discussions for practical applications:

  1. 1.

Using higher-order derivatives, i.e., larger $p$ in (14), can detect more
singular switching occurrence events for low-order derivatives, but the
detection results are more sensitive to the noise in the measured data which
may lead to faulty switching detection. Thus, a trade-off should be considered
for the choice of the order of derivatives. The decision is supposed to
heavily rely on the occurrence of singular switching events and the quality of
collected data.

  2. 2.

The threshold $\epsilon_{i}$ is crucial in the proposed framework. A large
threshold may neglect switching occurrence events, and on the other hand, a
small threshold may cause faulty switching detection. Thus, a trade-off should
be considered to determine a proper threshold in the proposed framework, e.g.,
a threshold can be chosen small enough to include all abrupt changes but also
should be larger than rest of values in $\mathbf{s}(k)$ by observation as in
the DC-DC converter example in Section IV. More advanced methods of threshold
selection will be our future work.

Considering multiple input-state traces $\phi_{K}^{(n)}$, $n=1,2,\ldots,N$, we
can obtain $N$ switching instant sets $\mathcal{S}_{n}$, $n=1,2,\ldots,N$ and
the switching set for the switched system is the union of these sets, which is
described as

$\displaystyle\mathcal{S}=\bigcup\nolimits_{n=1}^{N}\mathcal{S}_{n}$ (16)

where $\mathcal{S}=\\{k_{1},k_{2},\ldots,k_{\ell}\\}$, $0<k_{1}\leq\
k_{\ell}<K$ includes all switching instants detected from the collected
traces.

With the switching instant set $\mathcal{S}$ containing all detected switching
instants, each trace $\phi_{K}^{(n)}$ can be split to $\ell$ segments
$\\{\phi_{K,1}^{(n)},\phi_{K,2}^{(n)},\ldots,\phi_{K,\ell}^{(n)}\\}$ where
$\phi_{K,m}^{(n)}=\\{\mathbf{x}(k),\mathbf{u}(k)\\}$,
$k=k_{m},\ldots,k_{m+1}$, $m=1,2,\ldots,\ell$. Then, by consolidating the
$m$th segments of $N$ traces, the following segmented training data generated
from collected $N$ traces in the $K$-length time window can be obtained

$\displaystyle\Phi_{K,m}=\\{\phi_{K,m}^{(1)},\phi_{K,m}^{(2)},\ldots,\phi_{K,m}^{(N)}\\},~{}m=1,2,\ldots,\ell.$
(17)

These segmented data $\Phi_{K,m}$, $m=1,2,\ldots,\ell$ will be used for
training data-driven models of subsystems $f_{i}$, $i\in\mathcal{I}$.

###  III-B Subsystem Modeling

For each segmented traces $\Phi_{K,m}$, we can extract the training data of
input and output. The input data is $N(k_{m+1}-k_{m})$ pairs of
$\\{\mathbf{x}(k),\mathbf{u}(k)\\}$ and output data is $N(k_{m+1}-k_{m})$
state values $\mathbf{x}(k+1)$. Thus, the input-output training samples are

$\displaystyle\\{\tilde{\mathbf{x}}_{i,m},\mathbf{t}_{i,m}\\},~{}i=1,2,\ldots,\tilde{N}_{m}$
(18)

where $\tilde{\mathbf{x}}_{i,m}=\\{\mathbf{x}(k),\mathbf{u}(k)\\}$,
$\mathbf{t}_{i,m}=\mathbf{x}(k+1)$, $k=k_{m},\ldots,k_{m+1}$,
$\tilde{N}_{m}=N(k_{m+1}-k_{m})$.

Intuitively, we can use an ELM in the form of (2) to approximate system
dynamics in the interval $[k_{m},k_{m+1}]$ via solving the following
optimization problem

$\displaystyle\min_{\boldsymbol{\upbeta}\in\mathbb{R}^{L\times
n_{x}}}\left\|\mathbf{H}_{m}\boldsymbol{\upbeta}-\mathbf{T}_{m}\right\|$ (19)

where $\mathbf{H}_{m}$ is the randomized hidden layer output matrix based on
training input $\tilde{\mathbf{x}}_{i,m}$, and $\mathbf{T}_{m}$ is the
training data target matrix constructed based on $\mathbf{t}_{i,m}$,
respectively.

The above training process will generate $\ell$ optimal output weight vectors
$\boldsymbol{\upbeta}_{m}$ for $\ell$ individual segmented traces
$\Phi_{K,m}$. This will lead to $\ell$ ELMs for $\ell$ segmented traces
$\Phi_{K,m}$, $m=1,2,\ldots,\ell$, which is obviously not acceptable for the
most cases of switched system modeling since the number of subsystems is often
supposed to be much smaller than the switching instants. Therefore, the key
issue to solve the problem of subsystem modeling is how to properly merge
segmented traces for model training.

Considering two segmented traces $\Phi_{K,m_{1}}$ and $\Phi_{K,m_{2}}$, we
merge two segmented data into one trace set
$\\{\Phi_{K,m_{1}},\Phi_{K,m_{2}}\\}$ and the optimization problem is

$\displaystyle\min_{\boldsymbol{\upbeta}\in\mathbb{R}^{L\times
n_{x}}}\left\|\mathbf{H}_{m_{1},m_{2}}\boldsymbol{\upbeta}-\mathbf{T}_{m_{1},m_{2}}\right\|$
(20)

where $\mathbf{H}_{m_{1},m_{2}}$ is the randomized hidden layer output matrix
based on the combined training input $\tilde{\mathbf{x}}_{i,m}$ as

$\displaystyle\mathbf{H}_{m_{1},m_{2}}=\begin{bmatrix}\mathbf{H}_{m_{1}}\\\
\mathbf{H}_{m_{2}}\end{bmatrix}$ (21)

and $\mathbf{T}_{m_{1},m_{2}}$ is the combined training data target matrix

$\displaystyle\mathbf{T}_{m_{1},m_{2}}=\begin{bmatrix}\mathbf{T}_{m_{1}}\\\
\mathbf{T}_{m_{2}}\end{bmatrix}.$ (22)

Solving the above LS optimization problem to obtain the solution is
$\boldsymbol{\upbeta}_{m_{1},m_{2}}^{*}$, and the minimal value is

$\displaystyle\gamma(\boldsymbol{\upbeta}_{m_{1},m_{2}}^{*})=\left\|\mathbf{H}_{m_{1},m_{2}}\boldsymbol{\upbeta}_{m_{1},m_{2}}^{*}-\mathbf{T}_{m_{1},m_{2}}\right\|.$
(23)

If $\gamma(\boldsymbol{\upbeta}^{*}_{m_{1},m_{2}})$ is sufficient small, i.e.,
smaller than a prescribed threshold $\zeta>0$, it implies that there exist one
single common output weight vector $\boldsymbol{\upbeta}^{*}_{m_{1},m_{2}}$
that works for both segmented data $\Phi_{K,m_{1}}$ and $\Phi_{K,m_{2}}$ with
a small approximation error less than the prescribed threshold $\zeta$,
leading to a common ELM model extracted from two segmented data.

The segmented data merging process can be performed recursively from
$\Phi_{K,1}$ to $\Phi_{K,\ell}$ to eventually properly merge all segmented
data and model subsystem dynamics. The process is summarized in Algorithm 1.

###### Remark 5

In Algorithm 1, an inappropriately large threshold would wrongly merge
segmented data generated from different subsystems into one trace set, and on
the other hand, an inappropriately small threshold may fail to merge segmented
data from one same subsystem into one trace set. Thus, the prescribed
threshold $\zeta$ in the segmented data merging and subsystem modeling process
is crucial for the performance of Algorithm 1. Developing advanced methods of
threshold selection will be future work to enrich the modeling framework
proposed in this paper.

Input : Segmented trace data $\Phi_{K,m}$, $n=1,\ldots,\ell$

Output : Subsystem ELM models $f_{i}$, $i\in\tilde{\mathcal{I}}$

1 Function  _$\mathrm{systemModel}(\Phi_{K,m})$_

2

3 $n\leftarrow\ell$, $m_{1}\leftarrow 1$;

/* Segmented data merge */

4

5 while  _$m_{1} <n$_ do

6 $m_{2}\leftarrow 1$;

7 while  _$m_{2}\leq n$_ do

8 $m_{2}\leftarrow m_{2}+1$;

9 $\min_{\boldsymbol{\upbeta}\in\mathbb{R}^{L\times
n_{x}}}\left\|\mathbf{H}_{m_{1},m_{2}}\boldsymbol{\upbeta}-\mathbf{T}_{m_{1},m_{2}}\right\|$;

10 if  _$\gamma(\boldsymbol{\upbeta}^{*}_{m_{1},m_{2}})\leq\zeta$_ then

11 $\Phi_{K,m_{1}}\leftarrow\\{\Phi_{K,m_{1}},\Phi_{K,m_{2}}\\}$;

12 $n\leftarrow n-1$;

13 end if

14

15 end while

16 $m_{1}\leftarrow m_{1}+1$;

17 end while

/* Generate subsystem models */

18 $m\leftarrow 1$;

19 while  _$m\leq m_{1}$_ do

20 $\min_{\boldsymbol{\upbeta}\in\mathbb{R}^{L\times
n_{x}}}\left\|\mathbf{H}_{m}\boldsymbol{\upbeta}-\mathbf{T}_{m}\right\|$;

21 Obtain solution $\boldsymbol{\upbeta}^{*}_{m}$, $m=1,\ldots,m_{1}$ to
construct subsystem models $f_{i}$, $i=1,\ldots,m_{1}$ in the form of (2);

22 return _subsystem models $f_{i}$, $i=1,\ldots,m_{1}$_

23 end while

24

25

Algorithm 1 Subsystem Modeling

###  III-C Switching Reconstruction

As the input-state traces $\phi_{K}^{(n)}$, $n=1,2\ldots,N$ collected in a
finite-time window with a length of $K$ steps, we can infer time-dependent
switching laws valid for the $K$-length time interval. Thus, the time-
dependent switching law reconstruction for switching signal $\theta(k)$,
$k\in\mathbb{N}$ is reduced to reconstruction of time-dependent function
$\theta(k)$, $k\in[0,K]$ where $K$ is the length of training traces.

Using the switching detection results and subsystem modeling results by
Algorithm 1, we can label the time segments obtained by switching detection
with corresponding subsystem index obtained by Algorithm 1 to generate the
following sequence

$\displaystyle\\{\theta(k_{m}-1),\theta(k_{m})),k_{m}\\},~{}m=1,2,\ldots,\ell$
(24)

where $k_{m}$, $m=1,2,\ldots,\ell$ is the detected switching instants,
$\theta(k_{m}-1)$ and $\theta(k_{m})$ indicate the system switching from index
$\theta(k_{m}-1)$ to $\theta(k_{m})$.

If more information is provided, e.g., the switching law is a periodic or
semi-periodic switching, the switching law $\theta(k)$ can be inferred based
on the traces in the $K$-length time interval. A semi-periodic switching law
reconstruction is illustrated by the data-driven modeling of DC-DC converters
in the next section.

##  IV Evaluation by Data-Driven Modeling of DC-DC Converters

The DC-DC converter is a typical switched system. The example presented is
borrowed from [14], which is a boost converter with its topology shown in Fig.
1. The values of these circuit components are listed in Table I. When the
switch is _on_ , there is an increase in the inductor current $i_{L}(t)$. On
the other hand, when the switch is _off_ , the current ramps down through the
flyback diode, the capacitor, and the load. This leads to a transfer of the
energy accumulated in the inductance into the capacitor. The process executes
cyclically, whereas the boost converter operates with the switching period
$T_{s}$ and the duty cycle $u\in[0,1]$, which corresponds to the ratio of
activation duration of an _switch-on_ mode to the period.

Figure 1: DC-DC boost converter. TABLE I: Values of circuit components Component | $E$ | $L$ | $R_{L}$ | $C$ | $R_{C}$ | $T_{s}$
---|---|---|---|---|---|---
Value |  $20$ V |  $1$ mH |  $0.1$ $\Omega$ |  $10$ $\mu$F |  $0.06$ $\Omega$ |  $0.1$ ms

Corresponding to the modes of _switch-on_ and _switch-off_ , the DC-DC model
can be constructed. Letting system state
$\mathbf{x}(t)=[i_{L}(t)~{}~{}v(t)]^{\top}$, the state-space model of the DC-
DC converter is in the following form

$\displaystyle\dot{\mathbf{x}}(t)=\mathbf{A}_{\theta(t)}\mathbf{x}(t)+\mathbf{B}_{\theta(t)}$
(25)

where $\theta(t)={1,2}$ is the switching law defined by

$\displaystyle\theta(t)=\begin{cases}1,~{}t\in[2nT_{s},(2n+1)T_{s}+uT_{s})\\\
2,~{}t\in[(2n+1)T_{s}+uT_{s},2(n+1)T_{s})\end{cases}$ (26)

in which $n=0,1,2,\ldots$. System matrices are given as below:

$\displaystyle\mathbf{A}_{1}=\begin{bmatrix}{-\frac{R_{L}}{L}}&{0}\\\
{0}&{\frac{1}{(R+R_{C})C}}\end{bmatrix},~{}\mathbf{B}_{1}=\begin{bmatrix}\frac{E}{L}\\\
0\end{bmatrix}$
$\displaystyle\mathbf{A}_{2}=\begin{bmatrix}{-\frac{R_{L}}{L}-\frac{R_{C}R}{L(R+R_{C})}}&{-\frac{R}{R+R_{C}}L}\\\
{\frac{R}{R+R_{C}}C}&{-\frac{1}{R+R_{C}}C}\end{bmatrix},~{}\mathbf{B}_{2}=\begin{bmatrix}\frac{E}{L}\\\
0\end{bmatrix}.$

To evaluate the proposed data-driven modeling method for switched systems, we
first collect system traces generated by the DC-DC model, then utilize the
modeling framework proposed in this paper to develop a data-driven ELM model

$\mathbf{x}(k+1)=f_{\theta(k)}\left(\mathbf{x}(k)\right)$ (27)

where $f_{i}$ are ELMs and expected to be able to effectively predict the
system behaviors of the DC-DC converter.

With a sampling period of $0.01$ ms, we collect 20000 training data from 20
system traces in a $10$ ms time window generated from the DC-DC model (25),
then the data-driven modeling process is illustrated below111The source code
is available at: https://github.com/xiangweiming/data-driven-modeling-
switched-systems:

  1. 1.

Switching Detection: By executing switching detection with a prescribed
threshold $\epsilon_{1}=0.002$ for the average value of $\mathbf{s}(k)$ of 20
traces, the switching instants are detected as shown in Fig. 2. Training data
is also segmented for subsystem modeling.

  2. 2.

Subsystem Modeling: Using Algorithm 1 with a prescribed threshold $\zeta=1$
and ELM models with 200 neurons and sigmoid activation functions, two
subsystems of ELM are obtained.

  3. 3.

Switching Law Reconstruction: Based on the obtained switching instants and two
subsystems, the switching sequence
$\\{\theta(k_{m}-1),\theta(k_{m})),k_{m}\\},~{}m=1,2,\ldots,\ell$ is
reconstructed for interval $[0,10]$ ms as shown in Fig. 3, which can
successfully reconstruct the switching signal $\theta(k)$ described by (26).
Given the prior information that the switching signal is a semi-periodic
switching between subsystems, the switching law can be reconstructed as

$\displaystyle\theta(k)=\begin{cases}1,~{}k\in[20n,20n+11]\\\
2,~{}k\in[20n+11,20n+19]\end{cases}$

where $n=0,1,\ldots$.

Given the initial state $\mathbf{x}_{0}=[0.5~{}~{}0.5]^{\top}$, the state
responses and trajectories of data-driven switched system model are shown in
Figs. 4 and 5 which show data-driven model can predict the system trajectories
very well.

Figure 2: Switching detection. Figure 3: Reconstructed switching instants.
Figure 4: State responses of DC-DC model and data-driven model. Figure 5:
State trajectories of DC-DC model and data-driven model.

##  V Conclusion and Further Remarks

A data-driven modeling framework of switched dynamical systems using ELM has
been studied in this paper. An switching detection method is developed to
identify the time instants when switching occurs and generate segmented
training data. ELM is then utilized to learn subsystem dynamics. The modeling
algorithm is able to properly merge segmented data and learn subsystem models.
Finally, the switching sequence with both detected switching instants and
subsystem activation indices is reconstructed. The developed approach is
evaluated by an example of DC-DC converters.

Within the modeling framework of using ELM, there exist several challenging
problems remain open for further study. For instance, as aforementioned, the
selection of thresholds in switching detection and data merging is crucial in
the proposed modeling framework and deserves an in-depth study. In addition,
this paper assumed that the data has been well-denoised before modeling, it is
worthwhile to further explore the potentials of the proposed modeling
framework in presence of innegligible noise in measured data. It is also
meaningful to extend the proposed framework to state-dependent switched
systems in future work.

## References

  * [1] H. Lin and P. J. Antsaklis, “Stability and stabilizability of switched linear systems: a survey of recent results,” IEEE Transactions on Automatic control, vol. 54, no. 2, pp. 308–322, 2009.
  * [2] D. Liberzon and A. S. Morse, “Basic problems in stability and design of switched systems,” IEEE Control Systems Magazine, vol. 19, no. 5, pp. 59–70, 1999.
  * [3] W. Xiang, H.-D. Tran, and T. T. Johnson, “Nonconservative lifted convex conditions for stability of discrete-time switched systems under minimum dwell-time constraint,” IEEE Transactions on Automatic Control, vol. 64, pp. 3407–3414, 2019.
  * [4] W. Xiang, “Necessary and sufficient condition for stability of switched uncertain linear systems under dwell-time constraint,” IEEE Transactions on Automatic Control, vol. 61, no. 11, pp. 3619–3624, 2016.
  * [5] H. Ohlsson and L. Ljung, “Identification of switched linear regression models using sum-of-norms regularization,” Automatica, vol. 49, no. 4, pp. 1045–1050, 2013.
  * [6] R. Medhat, S. Ramesh, B. Bonakdarpour, and S. Fischmeister, “A framework for mining hybrid automata from input/output traces,” in 2015 International Conference on Embedded Software, pp. 177–186, IEEE, 2015.
  * [7] Q. Zhang, Q. Wang, and G. Li, “Switched system identification based on the constrained multi-objective optimization problem with application to the servo turntable,” International Journal of Control, Automation and Systems, vol. 14, no. 5, pp. 1153–1159, 2016.
  * [8] L. Bako, K. Boukharouba, E. Duviella, and S. Lecoeuche, “A recursive identification algorithm for switched linear/affine models,” Nonlinear Analysis: Hybrid Systems, vol. 5, no. 2, pp. 242–253, 2011.
  * [9] A. Garulli, S. Paoletti, and A. Vicino, “A survey on switched and piecewise affine system identification,” IFAC Proceedings Volumes, vol. 45, no. 16, pp. 344–355, 2012.
  * [10] S. Paoletti, A. L. Juloski, G. Ferrari-Trecate, and R. Vidal, “Identification of hybrid systems a tutorial,” European Journal of Control, vol. 13, no. 2-3, pp. 242–260, 2007.
  * [11] G.-B. Huang, L. Chen, and C. K. Siew, “Universal approximation using incremental constructive feedforward networks with random hidden nodes,” IEEE Transactions on Neural Networks, vol. 17, no. 4, pp. 879–892, 2006\.
  * [12] G.-B. Huang, Q.-Y. Zhu, and C.-K. Siew, “Extreme learning machine: theory and applications,” Neurocomputing, vol. 70, no. 1-3, pp. 489–501, 2006.
  * [13] K. Banerjee, “Generalized inverse of matrices and its applications,” 1973.
  * [14] W. Heemels, D. Lehmann, J. Lunze, and B. De Schutter, “Introduction to hybrid systems,” Handbook of Hybrid Systems Control–Theory, Tools, Applications, vol. 2, 2009.