# Octopus, a computational framework for exploring light-driven phenomena and quantum dynamics in extended and finite systems

Nicolas Tancogne-Dejean, Micael J. T. Oliveira, Xavier Andrade, Heiko

Appel, Carlos H. Borca, Guillaume Le Breton, Florian Buchholz, Alberto

Castro, Stefano Corni, Alfredo A. Correa, et al.

► To cite this version:

Nicolas Tancogne-Dejean, Micael J. T. Oliveira, Xavier Andrade, Heiko Appel, Carlos H. Borca, et al. Octopus, a computational framework for exploring light-driven phenomena and quantum dynamics in extended and finite systems. The Journal of Chemical Physics, 2020, Journal of Chemical Physics, 152, pp.124119. 10.1063/1.5142502. hal-04442664

HAL Id: hal-04442664

https://hal.univ-lille.fr/hal-04442664

Submitted on 14 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Octopus, a computational framework for exploring light-driven phenomena and quantum dynamics in extended and finite systems

Nicolas Tancogne-Dejean, Micael J. T. Oliveira, Xavier Andrade, Heiko Appel, Carlos H. Borca, Guillaume Le Breton, Florian Buchholz, Alberto Castro, Stefano Corni, Alfredo A. Correa, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Nicolas Tancogne-Dejean, Micael J. T. Oliveira, Xavier Andrade, Heiko Appel, Carlos H. Borca, et al. Octopus, a computational framework for exploring light-driven phenomena and quantum dynamics in extended and finite systems. The Journal of Chemical Physics, 2020, Journal of Chemical Physics, 152, pp.124119. 10.1063/1.5142502. hal-04442664

HAL Id: hal-04442664

https://lilloa.hal.science/hal-04442664v1

Submitted on 14 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Octopus, a computational framework for exploring light-driven phenomena and quantum dynamics in extended and finite systems

Special Collection: Chemical Physics Software Collection, Electronic Structure Software

Nicolas Tancogne-Dejean Ⓔ ; Micael J. T. Oliveira Ⓔ ; Xavier Andrade; Heiko Appel; Carlos H. Borca Ⓔ ;

Guillaume Le Breton; Florian Buchholz; Alberto Castro Ⓔ ; Stefano Corni; Alfredo A. Correa;

Umberto De Giovannini; Alain Delgado; Florian G. Eich; Johannes Flick; Gabriel Gil; Adrián Gomez Ⓔ ;

Nicole Helbig; Hannes Hübener; René Jestädt; Joaquim Jornet-Somoza; Ask H. Larsen Ⓔ ; Irina V. Lebedeva;

Martin Lüders; Miguel A. L. Marques Ⓔ ; Sebastian T. Ohlmann; Silvio Pipolo; Markus Rampp Ⓔ ;

Carlo A. Rozzi Ⓔ ; David A. Strubbe; Shunsuke A. Sato; Christian Schäfer Ⓔ ; Iris Theophilou Ⓔ ; Alicia Welden;

Angel Rubio ✉

<figure><img src="image_2.png" /></figure>

Check for updates

J. Chem. Phys. 152, 124119 (2020)

https://doi.org/10.1063/1.5142502

<figure><img src="image_3.png" /></figure>

View
Online

<figure><img src="image_4.png" /></figure>

Export
Citation

CrossMark

<figure><img src="image_5.png" /></figure>

AIP
Publishing

<figure><img src="image_6.png" /></figure>

Chemical Physics Reviews

Special Topic: Molecular Approaches
for Spin-based Technologies

Submit Today!

<figure><img src="image_7.png" /></figure>

AIP
Publishing
# Octopus, a computational framework for exploring light-driven phenomena and quantum dynamics in extended and finite systems

Special Collection: Chemical Physics Software Collection, Electronic Structure Software

Nicolas Tancogne-Dejean Ⓔ; Micael J. T. Oliveira Ⓔ; Xavier Andrade; Heiko Appel; Carlos H. Borca Ⓔ;

Guillaume Le Breton; Florian Buchholz; Alberto Castro Ⓔ; Stefano Corni; Alfredo A. Correa;

Umberto De Giovannini; Alain Delgado; Florian G. Eich; Johannes Flick; Gabriel Gil; Adrián Gomez Ⓔ;

Nicole Helbig; Hannes Hübener; René Jestädt; Joaquim Jornet-Somoza; Ask H. Larsen Ⓔ; Irina V. Lebedeva;

Martin Lüders; Miguel A. L. Marques Ⓔ; Sebastian T. Ohlmann; Silvio Pipolo; Markus Rampp Ⓔ;

Carlo A. Rozzi Ⓔ; David A. Strubbe; Shunsuke A. Sato; Christian Schäfer Ⓔ; Iris Theophilou Ⓔ; Alicia Welden;

Angel Rubio 📉

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

Check for updates

J. Chem. Phys. 152, 124119 (2020)

https://doi.org/10.1063/1.5142502

<figure><img src="image_5.png" /></figure>

View
Online

<figure><img src="image_6.png" /></figure>

Export
Citation

CrossMark

<figure><img src="image_7.png" /></figure>

AIP
Publishing

<figure><img src="image_8.png" /></figure>

Chemical Physics Reviews

Special Topic: Molecular Approaches
for Spin-based Technologies

Submit Today!

<figure><img src="image_9.png" /></figure>

AIP
Publishing
# Octopus, a computational framework for exploring light-driven phenomena and quantum dynamics in extended and finite systems

Cite as: J. Chem. Phys. 152, 124119 (2020); doi: 10.1063/1.5142502

Submitted: 17 December 2019 • Accepted: 1 March 2020 •

Published Online: 31 March 2020

<figure><img src="image_8.png" /></figure>

View Online

<figure><img src="image_9.png" /></figure>

Export Citation

<figure><img src="image_10.png" /></figure>

CrossMark

Nicolas Tancogne-Dejean,<fnref n="1" />,ª ID Micael J. T. Oliveira,<fnref n="1" />,ᵇ ID Xavier Andrade,<fnref n="2" /> Heiko Appel,<fnref n="1" /> Carlos H. Borca,<fnref n="2" />,<fnref n="3" /> ID

Guillaume Le Breton,<fnref n="6" /> Florian Buchholz,<fnref n="1" /> Alberto Castro,<fnref n="5" />,<fnref n="6" /> ID Stefano Corni,<fnref n="7" />,<fnref n="8" /> Alfredo A. Correa,<fnref n="2" />

Umberto De Giovannini,<fnref n="1" /> Alain Delgado,<fnref n="9" /> Florian G. Eich,<fnref n="1" /> Johannes Flick,<fnref n="10" />,<fnref n="11" /> Gabriel Gil,<fnref n="7" />,<fnref n="12" />

Adrián Gomez,<fnref n="5" /> ID Nicole Helbig,<fnref n="13" /> Hannes Hübener,<fnref n="1" /> René Jestädt,<fnref n="1" /> Joaquim Jornet-Somoza,<fnref n="1" />

Ask H. Larsen,<fnref n="14" /> ID Irina V. Lebedeva,<fnref n="14" /> Martin Lüders,<fnref n="1" /> Miguel A. L. Marques,<fnref n="15" /> ID Sebastian T. Ohlmann,<fnref n="16" />

Silvio Pipolo,<fnref n="17" /> Markus Rampp,<fnref n="16" /> ID Carlo A. Rozzi,<fnref n="8" /> ID David A. Strubbe,<fnref n="18" /> Shunsuke A. Sato,<fnref n="1" />,<fnref n="19" />

Christian Schäfer,<fnref n="1" /> ID Iris Theophilou,<fnref n="1" /> ID Alicia Welden,<fnref n="2" /> and Angel Rubio<fnref n="1" />,<fnref n="11" />,<fnref n="14" />,ᶜ ID

## AFFILIATIONS

<footnote n="1">¹ Max Planck Institute for the Structure and Dynamics of Matter, Luruper Chaussee 149, D-22761 Hamburg, Germany</footnote>

<footnote n="2">² Quantum Simulations Group, Lawrence Livermore National Laboratory, Livermore, California 94551, USA</footnote>

<footnote n="3">³ School of Chemistry and Biochemistry, Georgia Institute of Technology, Atlanta, Georgia 30332, USA</footnote>

<footnote n="4">⁴ Département de Physique, École Normale Supérieure de Lyon, 46 Allée d'Italie, Lyon Cedex 07, France</footnote>

<footnote n="5">⁵ Institute for Biocomputation and Physics of Complex Systems, University of Zaragoza, Calle Mariano Esquillor, 50018 Zaragoza, Spain</footnote>

<footnote n="6">⁶ ARAID Foundation, Avda. de Ranillas 1-D, 50018 Zaragoza, Spain</footnote>

<footnote n="7">⁷ Dipartimento di Scienze Chimiche, Università degli studi di Padova, via F. Marzolo 1, 35131 Padova, Italy</footnote>

<footnote n="8">⁸ CNR - Istituto Nanoscienze, via Campi 213a, 41125 Modena, Italy</footnote>

<footnote n="9">⁹ Xanadu, 777 Bay Street, Toronto, Ontario M5G 2C8, Canada</footnote>

<footnote n="10">¹⁰ John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, Massachusetts 02138, USA</footnote>

<footnote n="11">¹¹ Center for Computational Quantum Physics, Flatiron Institute, 162 5th Avenue, New York, New York 10010, USA</footnote>

<footnote n="12">¹² Instituto de Cibernética, Matemática y Física, Calle E 309, 10400 La Habana, Cuba</footnote>

<footnote n="13">¹³ Nanomat/Qmat/CESAM and ETSF, Université de Liège, B-4000 Sart-Tilman, Belgium</footnote>

<footnote n="14">¹⁴ Nano-Bio Spectroscopy Group and ETSF, Universidad del País Vasco, 20018 San Sebastián, Spain</footnote>

<footnote n="15">¹⁵ Institut für Physik, Martin-Luther-Universität Halle-Wittenberg, 06120 Halle (Saale), Germany</footnote>

<footnote n="16">¹⁶ Max Planck Computing and Data Facility, Gießenbachstraße 2, 85741 Garching, Germany</footnote>

<footnote n="17">¹⁷ Université de Lille, CNRS, Centrale Lille, ENSCL, Université d'Artois UMR 8181—UCCS Unité de Catalyse et Chimie du Solide, F-59000 Lille, France</footnote>

<footnote n="18">¹⁸ Department of Physics, School of Natural Sciences, University of California, Merced, California 95343, USA</footnote>

<footnote n="19">¹⁹ Center for Computational Sciences, University of Tsukuba, Tsukuba 305-8577, Japan</footnote>

<footnote>**Note:** This article is part of the JCP Special Topic on Electronic Structure Software.</footnote>

<footnote>a) Electronic mail: nicolas.tancogne-dejean@mpsd.mpg.de</footnote>

<footnote>b) Electronic mail: micael.oliveira@mpsd.mpg.de</footnote>

<footnote>c) Author to whom correspondence should be addressed: angel.rubio@mpsd.mpg.de</footnote>
# Octopus, a computational framework for exploring light-driven phenomena and quantum dynamics in extended and finite systems

Cite as: J. Chem. Phys. 152, 124119 (2020); doi: 10.1063/1.5142502

Submitted: 17 December 2019 • Accepted: 1 March 2020 •

Published Online: 31 March 2020

<figure><img src="image_10.png" /></figure>

View Online

<figure><img src="image_11.png" /></figure>

Export Citation

<figure><img src="image_12.png" /></figure>

CrossMark

Nicolas Tancogne-Dejean,<fnref n="1" />,<fnref n="4" /> Micael J. T. Oliveira,<fnref n="1" />,<fnref n="5" /> Xavier Andrade,<fnref n="2" /> Heiko Appel,<fnref n="1" /> Carlos H. Borca,<fnref n="2" />,<fnref n="3" />

Guillaume Le Breton,<fnref n="6" /> Florian Buchholz,<fnref n="1" /> Alberto Castro,<fnref n="5" />,<fnref n="6" /> Stefano Corni,<fnref n="7" />,<fnref n="8" /> Alfredo A. Correa,<fnref n="2" />

Umberto De Giovannini,<fnref n="1" /> Alain Delgado,<fnref n="9" /> Florian G. Eich,<fnref n="1" /> Johannes Flick,<fnref n="10" />,<fnref n="11" /> Gabriel Gil,<fnref n="7" />,<fnref n="12" />

Adrián Gomez,<fnref n="5" /> Nicole Helbig,<fnref n="13" /> Hannes Hübener,<fnref n="1" /> René Jestädt,<fnref n="1" /> Joaquim Jornet-Somoza,<fnref n="1" />

Ask H. Larsen,<fnref n="14" /> Irina V. Lebedeva,<fnref n="14" /> Martin Lüders,<fnref n="1" /> Miguel A. L. Marques,<fnref n="15" /> Sebastian T. Ohlmann,<fnref n="16" />

Silvio Pipolo,<fnref n="17" /> Markus Rampp,<fnref n="16" /> Carlo A. Rozzi,<fnref n="8" /> David A. Strubbe,<fnref n="18" /> Shunsuke A. Sato,<fnref n="1" />,<fnref n="19" />

Christian Schäfer,<fnref n="1" /> Iris Theophilou,<fnref n="1" /> Alicia Welden,<fnref n="2" /> and Angel Rubio<fnref n="1" />,<fnref n="11" />,<fnref n="14" />,ᶜ

## AFFILIATIONS

<footnote n="1">¹ Max Planck Institute for the Structure and Dynamics of Matter, Luruper Chaussee 149, D-22761 Hamburg, Germany</footnote>

<footnote n="2">² Quantum Simulations Group, Lawrence Livermore National Laboratory, Livermore, California 94551, USA</footnote>

<footnote n="3">³ School of Chemistry and Biochemistry, Georgia Institute of Technology, Atlanta, Georgia 30332, USA</footnote>

<footnote n="4">⁴ Département de Physique, École Normale Supérieure de Lyon, 46 Allée d'Italie, Lyon Cedex 07, France</footnote>

<footnote n="5">⁵ Institute for Biocomputation and Physics of Complex Systems, University of Zaragoza, Calle Mariano Esquillor, 50018 Zaragoza, Spain</footnote>

<footnote n="6">⁶ ARAID Foundation, Avda. de Ranillas 1-D, 50018 Zaragoza, Spain</footnote>

<footnote n="7">⁷ Dipartimento di Scienze Chimiche, Università degli studi di Padova, via F. Marzolo 1, 35131 Padova, Italy</footnote>

<footnote n="8">⁸ CNR – Istituto Nanoscienze, via Campi 213a, 41125 Modena, Italy</footnote>

<footnote n="9">⁹ Xanadu, 777 Bay Street, Toronto, Ontario M5G 2C8, Canada</footnote>

<footnote n="10">¹⁰ John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, Massachusetts 02138, USA</footnote>

<footnote n="11">¹¹ Center for Computational Quantum Physics, Flatiron Institute, 162 5th Avenue, New York, New York 10010, USA</footnote>

<footnote n="12">¹² Instituto de Cibernética, Matemática y Física, Calle E 309, 10400 La Habana, Cuba</footnote>

<footnote n="13">¹³ Nanomat/Qmat/CESAM and ETSF, Université de Liège, B-4000 Sart-Tilman, Belgium</footnote>

<footnote n="14">¹⁴ Nano-Bio Spectroscopy Group and ETSF, Universidad del País Vasco, 20018 San Sebastián, Spain</footnote>

<footnote n="15">¹⁵ Institut für Physik, Martin-Luther-Universität Halle-Wittenberg, 06120 Halle (Saale), Germany</footnote>

<footnote n="16">¹⁶ Max Planck Computing and Data Facility, Gießenbachstraße 2, 85741 Garching, Germany</footnote>

<footnote n="17">¹⁷ Université de Lille, CNRS, Centrale Lille, ENSCL, Université d'Artois UMR 8181—UCCS Unité de Catalyse et Chimie du Solide, F-59000 Lille, France</footnote>

<footnote n="18">¹⁸ Department of Physics, School of Natural Sciences, University of California, Merced, California 95343, USA</footnote>

<footnote n="19">¹⁹ Center for Computational Sciences, University of Tsukuba, Tsukuba 305-8577, Japan</footnote>

<footnote>Note: This article is part of the JCP Special Topic on Electronic Structure Software.</footnote>

<footnote>a) Electronic mail: nicolas.tancogne-dejean@mpsd.mpg.de</footnote>

<footnote>b) Electronic mail: micael.oliveira@mpsd.mpg.de</footnote>

<footnote>c) Author to whom correspondence should be addressed: angel.rubio@mpsd.mpg.de</footnote>
## ABSTRACT

Over the last few years, extraordinary advances in experimental and theoretical tools have allowed us to monitor and control matter at short time and atomic scales with a high degree of precision. An appealing and challenging route toward engineering materials with tailored properties is to find ways to design or selectively manipulate materials, especially at the quantum level. To this end, having a state-of-the-art *ab initio* computer simulation tool that enables a reliable and accurate simulation of light-induced changes in the physical and chemical properties of complex systems is of utmost importance. The first principles real-space-based Octopus project was born with that idea in mind, i.e., to provide a unique framework that allows us to describe non-equilibrium phenomena in molecular complexes, low dimensional materials, and extended systems by accounting for electronic, ionic, and photon quantum mechanical effects within a generalized time-dependent density functional theory. This article aims to present the new features that have been implemented over the last few years, including technical developments related to performance and massive parallelism. We also describe the major theoretical developments to address ultrafast light-driven processes, such as the new theoretical framework of quantum electrodynamics density-functional formalism for the description of novel light-matter hybrid states. Those advances, and others being released soon as part of the Octopus package, will allow the scientific community to simulate and characterize spatial and time-resolved spectroscopies, ultrafast phenomena in molecules and materials, and new emergent states of matter (quantum electrodynamical-materials).

© 2020 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/1.5142502

## I. INTRODUCTION

It is a general challenge in the electronic structure community to develop accurate and efficient methods for modeling materials of ever increasing complexity in order to predict their properties. In this respect, time-dependent density functional theory (TDDFT) and related methods have become a natural choice for modeling systems in and out of their equilibrium.

During the last few years, novel directions of research emerged in the fields of chemistry, physics, and materials science that required the development of novel simulation tools needed to face the new challenges posed by those experimental advances and emerging phenomena. One can cite some examples among these new fields of research: the strong coupling between light and matter (including materials embedded in cavities), new states of matter (hidden phases and topological solids and molecules), and strong field dynamics in periodic systems. Whereas the strong-field dynamics of atoms and molecules is now well understood, strong-field dynamics in solids is an active field of research. Real-time (rt) TDDFT¹,² represents a natural tool to study highly non-linear phenomena in solids and low-dimensional materials, and the development of efficient numerical methods to perform real-time TDDFT in such periodic or semi-periodic systems is crucial to explore this new phenomenon (including the quantum nature of light and phonons). Indeed, it allows us to describe highly nonlinear processes without having to resort to perturbation theory on top of equilibrium DFT calculations.

Recent years have seen tremendous experimental progress in the field of strong light-matter interactions,³⁻⁵ where the strong coupling of light to chemical systems, quantum materials, or nanoplasmonic systems, among others, has been demonstrated. In this regime, light and matter meet on the same footing and the electron-photon interaction has to be explicitly considered.⁴⁻⁹ A novel theoretical approach accelerating this field is quantum-electrodynamical density-functional theory (QEDFT),⁴,⁵,¹⁰⁻¹⁴ which complements TDDFT with the photonic degrees of freedom and provides reliable and predictive simulations in this emerging field of research.

In this paper, we explore the recent advances in the Octopus project.¹⁵⁻¹⁹ We focus our attention on the recently added features and, particularly, on the ones that have not been described in the previous papers.¹⁷⁻¹⁹ These new features include the implementation of new levels of self-consistent and microscopic couplings of light and matter, the treatment of solvent effects through the polarizable continuum model (PCM), the implementation of various methods to treat van der Waals (vdW) interactions, new methods to calculate magnons, conductivities, and photoelectron spectroscopy from real-time TDDFT, and the calculation of orbital magneto-optical responses. Advances in numerical algorithms and methods, such as new propagators and the use of iterative eigensolvers in the context of reduced density matrix (RDM) functional theory, are also discussed. Finally, recent improvements in the treatment of periodic systems, as well as more technical code improvements, are also presented. For a more detailed explanation about how to use these newly introduced features in practice, we refer the readers to the Octopus webpage,²⁰ as new tutorials and examples are regularly being added to it.

This paper is organized as follows: First, we present in Secs. II–XII the new implementations of physical theories and algorithms that allow us to deal with non-equilibrium phenomena in materials and nanostructures. This is followed by Secs. XIII–XVI dealing with technical developments that are fundamental in improving the code performance and the stability of the algorithms. Finally, we draw our conclusions in Sec. XVII. Unless otherwise stated, atomic units are used throughout the paper.

## II. COUPLED MAXWELL-KOHN-SHAM EQUATIONS

In most cases when light-matter interactions are considered, a decoupling of light and matter is performed at the outset. Either the electromagnetic fields are prescribed and then the properties of the matter subsystem are determined, as frequently found in, e.g., quantum chemistry or solid-state physics, or the properties of matter are prescribed and then the properties of the photon subsystem are determined, as done in, e.g., quantum optics or photonics.¹²,²¹
## ABSTRACT

Over the last few years, extraordinary advances in experimental and theoretical tools have allowed us to monitor and control matter at short time and atomic scales with a high degree of precision. An appealing and challenging route toward engineering materials with tailored properties is to find ways to design or selectively manipulate materials, especially at the quantum level. To this end, having a state-of-the-art *ab initio* computer simulation tool that enables a reliable and accurate simulation of light-induced changes in the physical and chemical properties of complex systems is of utmost importance. The first principles real-space-based Octopus project was born with that idea in mind, i.e., to provide a unique framework that allows us to describe non-equilibrium phenomena in molecular complexes, low dimensional materials, and extended systems by accounting for electronic, ionic, and photon quantum mechanical effects within a generalized time-dependent density functional theory. This article aims to present the new features that have been implemented over the last few years, including technical developments related to performance and massive parallelism. We also describe the major theoretical developments to address ultrafast light-driven processes, such as the new theoretical framework of quantum electrodynamics density-functional formalism for the description of novel light-matter hybrid states. Those advances, and others being released soon as part of the Octopus package, will allow the scientific community to simulate and characterize spatial and time-resolved spectroscopies, ultrafast phenomena in molecules and materials, and new emergent states of matter (quantum electrodynamical-materials).

© 2020 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/1.5142502

## I. INTRODUCTION

It is a general challenge in the electronic structure community to develop accurate and efficient methods for modeling materials of ever increasing complexity in order to predict their properties. In this respect, time-dependent density functional theory (TDDFT) and related methods have become a natural choice for modeling systems in and out of their equilibrium.

During the last few years, novel directions of research emerged in the fields of chemistry, physics, and materials science that required the development of novel simulation tools needed to face the new challenges posed by those experimental advances and emerging phenomena. One can cite some examples among these new fields of research: the strong coupling between light and matter (including materials embedded in cavities), new states of matter (hidden phases and topological solids and molecules), and strong field dynamics in periodic systems. Whereas the strong-field dynamics of atoms and molecules is now well understood, strong-field dynamics in solids is an active field of research. Real-time (rt) TDDFT¹,² represents a natural tool to study highly non-linear phenomena in solids and low-dimensional materials, and the development of efficient numerical methods to perform real-time TDDFT in such periodic or semi-periodic systems is crucial to explore this new phenomenon (including the quantum nature of light and phonons). Indeed, it allows us to describe highly nonlinear processes without having to resort to perturbation theory on top of equilibrium DFT calculations.

Recent years have seen tremendous experimental progress in the field of strong light-matter interactions,³⁻⁵ where the strong coupling of light to chemical systems, quantum materials, or nanoplasmonic systems, among others, has been demonstrated. In this regime, light and matter meet on the same footing and the electron-photon interaction has to be explicitly considered.⁴⁻⁹ A novel theoretical approach accelerating this field is quantum-electrodynamical density-functional theory (QEDFT),⁴,⁵,¹⁰⁻¹⁴ which complements TDDFT with the photonic degrees of freedom and provides reliable and predictive simulations in this emerging field of research.

In this paper, we explore the recent advances in the Octopus project.¹⁵⁻¹⁹ We focus our attention on the recently added features and, particularly, on the ones that have not been described in the previous papers.¹⁷⁻¹⁹ These new features include the implementation of new levels of self-consistent and microscopic couplings of light and matter, the treatment of solvent effects through the polarizable continuum model (PCM), the implementation of various methods to treat van der Waals (vdW) interactions, new methods to calculate magnons, conductivities, and photoelectron spectroscopy from real-time TDDFT, and the calculation of orbital magneto-optical responses. Advances in numerical algorithms and methods, such as new propagators and the use of iterative eigensolvers in the context of reduced density matrix (RDM) functional theory, are also discussed. Finally, recent improvements in the treatment of periodic systems, as well as more technical code improvements, are also presented. For a more detailed explanation about how to use these newly introduced features in practice, we refer the readers to the Octopus webpage,²⁰ as new tutorials and examples are regularly being added to it.

This paper is organized as follows: First, we present in Secs. II–XII the new implementations of physical theories and algorithms that allow us to deal with non-equilibrium phenomena in materials and nanostructures. This is followed by Secs. XIII–XVI dealing with technical developments that are fundamental in improving the code performance and the stability of the algorithms. Finally, we draw our conclusions in Sec. XVII. Unless otherwise stated, atomic units are used throughout the paper.

## II. COUPLED MAXWELL-KOHN-SHAM EQUATIONS

In most cases when light-matter interactions are considered, a decoupling of light and matter is performed at the outset. Either the electromagnetic fields are prescribed and then the properties of the matter subsystem are determined, as frequently found in, e.g., quantum chemistry or solid-state physics, or the properties of matter are prescribed and then the properties of the photon subsystem are determined, as done in, e.g., quantum optics or photonics.¹²,²¹
The microscopic interaction of light and matter in Octopus has followed this decoupling strategy and has been treated so far only in the forward-coupling direction. In this approximation, an external classical laser pulse or kick is prescribed and the response of the system is computed by evolving the Kohn–Sham orbitals.²² The back-action of the matter subsystem on the laser pulse, the subsequent effect of this modified pulse on the matter, and so on are ignored. This forward-coupling approximation is highly accurate when the total current generated in the system is comparatively small, such as in atoms or small molecules. Small refers here to currents that induce transverse electric fields that are about an order of magnitude smaller than the longitudinal electric field that originates from the Hartree potential of the system. This has been exploited in Octopus and other codes, and many different types of spectroscopies were computed successfully in the past using this approach.

In contrast, in classical electromagnetic modeling, the opposite view is taken. Here, the material properties are prescribed and the resulting electromagnetic fields are computed. In practice, the material properties are routinely approximated by a local continuum model or by the dielectric function of the system (such as the Debye, Lorentz, or Drude models) and then Maxwell's equations are solved for linear dielectric media arranged in appropriate geometries.²³

It is clear that both perspectives, a focus on matter dynamics alone or a focus on electromagnetic field dynamics alone, break down when the total induced currents become large and when electromagnetic near-field effects on the scale of the material system are not negligible anymore. Prime examples for such cases are nanoplasmonic systems, surface plasmon-polaritons, or tip-enhanced spectroscopies. In these cases, the back-action of the material response on the system itself has to be taken into account, leading to screening and retardation effects. The proper theoretical framework, which encompasses all these effects, is quantum electrodynamics. Starting with a generalized Pauli-Fierz field theory for the combined system of electrons, nuclei, and photons, we have recently derived different levels of self-consistent and microscopic couplings of light and matter. In the classical limit, this results in a coupled set of Ehrenfest–Maxwell–Pauli–Kohn–Sham equations.²¹ To implement these equations, we added a Maxwell solver to the Octopus code, which we couple self-consistently to the dynamics of the electrons and nuclei. In the following, we briefly summarize the basic ingredients for this implementation and show an example of self-consistent light-matter interactions for a nanoplasmonic system. Further details of the implementation and nano-optical applications can be found in Ref. 21.

Since over the years, Octopus has been optimized heavily to solve time-dependent Schrödinger and Kohn–Sham equations, we have exploited the fact that Maxwell's equations can be formulated in Schrödinger form²⁴ to benefit from the efficient time-evolution in the code. This reformulation is based on the Riemann–Silberstein vector,²⁵ which is a combination of the electric E(r, t) and magnetic field B(r, t),

$$F_{\pm}(\mathbf{r}, t) = \sqrt{\frac{\epsilon_0}{2}} E(\mathbf{r}, t) \pm i \sqrt{\frac{1}{2\mu_0}} \mathbf{B}(\mathbf{r}, t). \quad (1)$$

The sign of the imaginary part of the Riemann–Silberstein vector corresponds to different helicities. The reformulation of Maxwell's

equations in Schrödinger form is purely algebraic and starts out with the microscopic Maxwell's equations

$$\nabla \cdot \mathbf{E} = \frac{\rho}{\epsilon_0}, \qquad \nabla \cdot \mathbf{B} = 0, \qquad (2)$$

$$\nabla \times \mathbf{B} = \frac{1}{c^2} \frac{\partial \mathbf{E}}{\partial t} + \mu_0 \mathbf{J}, \qquad \nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t}, \qquad (3)$$

where E and B are the classical electric and magnetic fields, respectively, ρ and J are the charge and current densities, ε₀ and μ₀ are the vacuum permittivity and permeability, respectively, and c = (ε₀μ₀)⁻¹/² is the speed of light. Using the Riemann–Silberstein vector, the electric and magnetic Gauss laws may now be combined,

$$\nabla \cdot \mathbf{F} = \frac{1}{\sqrt{2\epsilon_0}} \rho, \qquad (4)$$

and likewise, the Faraday and Ampere laws can be combined into one evolution equation for the Riemann–Silberstein vector,

$$i\hbar \frac{\partial \mathbf{F}}{\partial t} = c \left( \mathbf{S} \cdot \frac{\hbar}{i} \nabla \right) \mathbf{F} - \frac{i\hbar}{\sqrt{2\epsilon_0}} \mathbf{J}. \qquad (5)$$

Here, S = (Sₓ, Sᵧ, S_z) denotes a vector of spin one matrices

$$S_x = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & -i \\ 0 & i & 0 \end{pmatrix}, \quad S_y = \begin{pmatrix} 0 & 0 & i \\ 0 & 0 & 0 \\ -i & 0 & 0 \end{pmatrix}, \quad S_z = \begin{pmatrix} 0 & -i & 0 \\ i & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad (6)$$

which are analogous to the Pauli matrices and show the spin-one character of the photon. Having cast Maxwell's equations as an inhomogeneous Schrödinger equation, it is now straightforward to use the time-evolution algorithms in Octopus to time-evolve the Riemann–Silberstein vector. The only difference to the matter propagation is that we are now dealing with the “Maxwell Hamiltonian”

$$H_{\text{EM}} = c \left( \mathbf{S} \cdot \frac{\hbar}{i} \nabla \right), \qquad (7)$$

which acts upon six orbitals of the Riemann–Silberstein vector corresponding to the three components of the electric and magnetic field vectors. As in the matter case in Octopus, the discretization of the gradient in the Maxwell Hamiltonian is performed with finite-difference stencils, and the domain parallelization of Octopus can be used seamlessly for the Maxwell case as well. The difference to finite-difference time-domain (FDTD) codes based on the Yee algorithm is that we do not employ two shifted grids for the electric and magnetic fields,²³ but rather a single grid for the Riemann–Silberstein vector. This simplifies the coupling to matter and allows us to use higher-order finite-difference discretizations for the gradient. Since the spatial discretization is connected to the temporal discretization through the Courant condition, this, in turn, allows us to use larger time steps. Furthermore, from our experience, a unified grid also improves the stability compared to FDTD.

Instead of using the constitutive relations, we couple Maxwell's equations directly to the microscopic current density of the matter subsystem, consisting of the usual paramagnetic current term, the diamagnetic current term, and the magnetization current term. For the coupling of the electromagnetic fields to the matter subsystem,
The microscopic interaction of light and matter in Octopus has followed this decoupling strategy and has been treated so far only in the forward-coupling direction. In this approximation, an external classical laser pulse or kick is prescribed and the response of the system is computed by evolving the Kohn–Sham orbitals.²² The back-action of the matter subsystem on the laser pulse, the subsequent effect of this modified pulse on the matter, and so on are ignored. This forward-coupling approximation is highly accurate when the total current generated in the system is comparatively small, such as in atoms or small molecules. Small refers here to currents that induce transverse electric fields that are about an order of magnitude smaller than the longitudinal electric field that originates from the Hartree potential of the system. This has been exploited in Octopus and other codes, and many different types of spectroscopies were computed successfully in the past using this approach.

In contrast, in classical electromagnetic modeling, the opposite view is taken. Here, the material properties are prescribed and the resulting electromagnetic fields are computed. In practice, the material properties are routinely approximated by a local continuum model or by the dielectric function of the system (such as the Debye, Lorentz, or Drude models) and then Maxwell's equations are solved for linear dielectric media arranged in appropriate geometries.²³

It is clear that both perspectives, a focus on matter dynamics alone or a focus on electromagnetic field dynamics alone, break down when the total induced currents become large and when electromagnetic near-field effects on the scale of the material system are not negligible anymore. Prime examples for such cases are nanoplasmonic systems, surface plasmon-polaritons, or tip-enhanced spectroscopies. In these cases, the back-action of the material response on the system itself has to be taken into account, leading to screening and retardation effects. The proper theoretical framework, which encompasses all these effects, is quantum electrodynamics. Starting with a generalized Pauli-Fierz field theory for the combined system of electrons, nuclei, and photons, we have recently derived different levels of self-consistent and microscopic couplings of light and matter. In the classical limit, this results in a coupled set of Ehrenfest–Maxwell–Pauli–Kohn–Sham equations.²¹ To implement these equations, we added a Maxwell solver to the Octopus code, which we couple self-consistently to the dynamics of the electrons and nuclei. In the following, we briefly summarize the basic ingredients for this implementation and show an example of self-consistent light–matter interactions for a nanoplasmonic system. Further details of the implementation and nano-optical applications can be found in Ref. 21.

Since over the years, Octopus has been optimized heavily to solve time-dependent Schrödinger and Kohn–Sham equations, we have exploited the fact that Maxwell's equations can be formulated in Schrödinger form²⁴ to benefit from the efficient time-evolution in the code. This reformulation is based on the Riemann–Silberstein vector,²⁵ which is a combination of the electric E(r, t) and magnetic field B(r, t),

$$F_{\pm}(r, t) = \sqrt{\frac{\epsilon_0}{2}} E(r, t) \pm i \sqrt{\frac{1}{2\mu_0}} B(r, t). \quad (1)$$

The sign of the imaginary part of the Riemann–Silberstein vector corresponds to different helicities. The reformulation of Maxwell's

equations in Schrödinger form is purely algebraic and starts out with the microscopic Maxwell's equations

$$\nabla \cdot \mathbf{E} = \frac{\rho}{\epsilon_0}, \qquad \nabla \cdot \mathbf{B} = 0, \quad (2)$$

$$\nabla \times \mathbf{B} = \frac{1}{c^2} \frac{\partial \mathbf{E}}{\partial t} + \mu_0 \mathbf{J}, \qquad \nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t}, \quad (3)$$

where E and B are the classical electric and magnetic fields, respectively, ρ and J are the charge and current densities, ε₀ and μ₀ are the vacuum permittivity and permeability, respectively, and c = (ε₀μ₀)⁻¹/² is the speed of light. Using the Riemann–Silberstein vector, the electric and magnetic Gauss laws may now be combined,

$$\nabla \cdot \mathbf{F} = \frac{1}{\sqrt{2\epsilon_0}} \rho, \quad (4)$$

and likewise, the Faraday and Ampere laws can be combined into one evolution equation for the Riemann–Silberstein vector,

$$i\hbar \frac{\partial \mathbf{F}}{\partial t} = c \left( \mathbf{S} \cdot \frac{\hbar}{i} \nabla \right) \mathbf{F} - \frac{i\hbar}{\sqrt{2\epsilon_0}} \mathbf{J}. \quad (5)$$

Here, S = (Sₓ, Sᵧ, S_z) denotes a vector of spin one matrices

$$S_x = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & -i \\ 0 & i & 0 \end{pmatrix}, \quad S_y = \begin{pmatrix} 0 & 0 & i \\ 0 & 0 & 0 \\ -i & 0 & 0 \end{pmatrix}, \quad S_z = \begin{pmatrix} 0 & -i & 0 \\ i & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad (6)$$

which are analogous to the Pauli matrices and show the spin-one character of the photon. Having cast Maxwell's equations as an inhomogeneous Schrödinger equation, it is now straightforward to use the time-evolution algorithms in Octopus to time-evolve the Riemann–Silberstein vector. The only difference to the matter propagation is that we are now dealing with the “Maxwell Hamiltonian”

$$H_{\text{EM}} = c \left( \mathbf{S} \cdot \frac{\hbar}{i} \nabla \right), \quad (7)$$

which acts upon six orbitals of the Riemann–Silberstein vector corresponding to the three components of the electric and magnetic field vectors. As in the matter case in Octopus, the discretization of the gradient in the Maxwell Hamiltonian is performed with finite-difference stencils, and the domain parallelization of Octopus can be used seamlessly for the Maxwell case as well. The difference to finite-difference time-domain (FDTD) codes based on the Yee algorithm is that we do not employ two shifted grids for the electric and magnetic fields,²³ but rather a single grid for the Riemann–Silberstein vector. This simplifies the coupling to matter and allows us to use higher-order finite-difference discretizations for the gradient. Since the spatial discretization is connected to the temporal discretization through the Courant condition, this, in turn, allows us to use larger time steps. Furthermore, from our experience, a unified grid also improves the stability compared to FDTD.

Instead of using the constitutive relations, we couple Maxwell's equations directly to the microscopic current density of the matter subsystem, consisting of the usual paramagnetic current term, the diamagnetic current term, and the magnetization current term. For the coupling of the electromagnetic fields to the matter subsystem,
we are relying on the Power–Zienau–Woolley transformation,²⁶,²⁷
which leads to a multipole expansion. We have implemented the
first two orders of this expansion: the dipole approximation in low-
est order and electrical-quadrupole and magnetic-dipole coupling in
the next order. In addition, we are currently working on implement-
ing the full minimal coupling with a full position dependent vector
potential.

The time evolution of the Kohn–Sham orbitals and the Maxwell
fields is performed side-by-side and the two subsystems are coupled
self-consistently in each time step, as illustrated in Fig. 1. To prop-
agate different subsystems with different Hamiltonians and differ-
ent sets of orbitals simultaneously, we implemented a multi-system
framework in Octopus (for more details, refer Subsection XVI C).

Similar to the matter propagation, also in the Maxwell case, out-
going waves that reach the boundary of the simulation box have to be
absorbed to avoid artificial reflections and backscattering. Our first
attempt was to use also the mask functions that are used for the mat-
ter propagation in Octopus. However, in the electromagnetic case,
the mask absorption of outgoing waves turned out to be not efficient
enough so that we implemented a perfectly matched layer (PML)²³
for the Maxwell propagation.

When considering incoming electromagnetic fields with opti-
cal wavelengths, the coupling to atomistic or nano-scale systems
leads to a multi-scale problem. The optical wavelength of the radi-
ation is in this case much larger than the de Broglie wavelength
of the matter. Likewise, the electromagnetic waves are traveling
with the speed of light, which requires sub-attosecond time steps.
We have, therefore, implemented different multi-scale couplings in
space and time. For example, the Maxwell simulation box can be
on the same scale as the matter box. In this case, the electromag-
netic waves are represented as incoming analytical time-dependent
boundary conditions and are propagated numerically inside the sim-
ulation box. Alternatively, the Maxwell simulation box can be much
larger than the matter box to fully encompass laser pulses with
optical wavelengths. In this case, prolongations and interpolations
have to be used similarly to multi-grid methods. Since the electronic
and nuclear motion is much slower than the time-evolution of the

electromagnetic waves, we also implemented a multi-scale approach
for the real-time propagation. The Riemann–Silberstein vector is
propagated with frozen electronic current from the last point of
interaction for many intermediate time steps before a coupling to
the matter subsystem takes place. The number of intermediate steps
is a convergence parameter and depends on the physical situation
at hand.

Since we now include the description of classical electromag-
netic fields explicitly in our real-time simulations, we have directly
access to the outgoing electromagnetic radiation. This allows us
to define electromagnetic detectors at the box boundaries, which
accumulate the outgoing electromagnetic waves. We implemented
such electromagnetic detectors in Octopus, and this allows us to
run simulations in close analogy with experiments and to directly
observe the outgoing radiation. For example, it is no longer needed
to Fourier transform the time signal of the matter dipole to get opti-
cal spectra, but we rather have access to the spectrum directly on the
Maxwell grid.

As an example of a coupled Ehrenfest–Maxwell–Pauli–Kohn–
Sham propagation with Octopus, we selected a nano-optical appli-
cation. We consider in this example two almost spherical sodium
nanoparticles with 297 sodium atoms each, which are arranged in
a dimer configuration. This system interacts with an incoming laser
pulse that excites either the internal dipole or quadrupole plasmon
motion of the dimer. In Fig. 2, we show the resulting electromag-
netic field enhancements for different levels of light–matter cou-
pling. In panels (a) and (b), we show the temporal profile of the
incoming laser pulse and the resulting current density at the cen-
ter point between the two nanoparticles. The field enhancement can
be seen in panels (c)–(e). Including a self-consistent back reaction
in the light–matter coupling, the field enhancement is reduced at
the center point between the two nanoparticles, as shown in panel
(c), while far away from the dimer, as shown in panels (d) and (e),
the field enhancement is larger than in the forward coupled case.
Furthermore, frequency shifts can also be observed, which are more
pronounced in the near-field. We found that the field enhancement
is also sensitive to the coupling terms of the multipole expansion,

<figure><img src="image_11.png" /><figcaption>FIG. 1. (Left) The standard forward-coupling approximation: the electromagnetic fields (in blue) propagate freely and only influence the propagation of the matter (in red). The back reaction of the matter currents on the electromagnetic fields is neglected. (Right) Illustration of a fully self-consistent predictor-corrector scheme for a coupled Maxwell–Pauli–Kohn–Sham time step. As before, the electromagnetic fields influence the propagation of the matter (forward coupling). However, in addition, here the currents from the matter propagation also influence the propagation of the electromagnetic fields (backward coupling). A given time step for the matter wavefunctions and the electromagnetic fields is repeated until self-consistency is found (self-consistent forward–backward coupling).</figcaption></figure>
we are relying on the Power–Zienau–Woolley transformation,²⁶,²⁷
which leads to a multipole expansion. We have implemented the
first two orders of this expansion: the dipole approximation in low-
est order and electrical-quadrupole and magnetic-dipole coupling in
the next order. In addition, we are currently working on implement-
ing the full minimal coupling with a full position dependent vector
potential.

The time evolution of the Kohn–Sham orbitals and the Maxwell
fields is performed side-by-side and the two subsystems are coupled
self-consistently in each time step, as illustrated in Fig. 1. To prop-
agate different subsystems with different Hamiltonians and differ-
ent sets of orbitals simultaneously, we implemented a multi-system
framework in Octopus (for more details, refer Subsection XVI C).

Similar to the matter propagation, also in the Maxwell case, out-
going waves that reach the boundary of the simulation box have to be
absorbed to avoid artificial reflections and backscattering. Our first
attempt was to use also the mask functions that are used for the mat-
ter propagation in Octopus. However, in the electromagnetic case,
the mask absorption of outgoing waves turned out to be not efficient
enough so that we implemented a perfectly matched layer (PML)²³
for the Maxwell propagation.

When considering incoming electromagnetic fields with opti-
cal wavelengths, the coupling to atomistic or nano-scale systems
leads to a multi-scale problem. The optical wavelength of the radi-
ation is in this case much larger than the de Broglie wavelength
of the matter. Likewise, the electromagnetic waves are traveling
with the speed of light, which requires sub-attosecond time steps.
We have, therefore, implemented different multi-scale couplings in
space and time. For example, the Maxwell simulation box can be
on the same scale as the matter box. In this case, the electromag-
netic waves are represented as incoming analytical time-dependent
boundary conditions and are propagated numerically inside the sim-
ulation box. Alternatively, the Maxwell simulation box can be much
larger than the matter box to fully encompass laser pulses with
optical wavelengths. In this case, prolongations and interpolations
have to be used similarly to multi-grid methods. Since the electronic
and nuclear motion is much slower than the time-evolution of the

electromagnetic waves, we also implemented a multi-scale approach
for the real-time propagation. The Riemann–Silberstein vector is
propagated with frozen electronic current from the last point of
interaction for many intermediate time steps before a coupling to
the matter subsystem takes place. The number of intermediate steps
is a convergence parameter and depends on the physical situation
at hand.

Since we now include the description of classical electromag-
netic fields explicitly in our real-time simulations, we have directly
access to the outgoing electromagnetic radiation. This allows us
to define electromagnetic detectors at the box boundaries, which
accumulate the outgoing electromagnetic waves. We implemented
such electromagnetic detectors in Octopus, and this allows us to
run simulations in close analogy with experiments and to directly
observe the outgoing radiation. For example, it is no longer needed
to Fourier transform the time signal of the matter dipole to get opti-
cal spectra, but we rather have access to the spectrum directly on the
Maxwell grid.

As an example of a coupled Ehrenfest–Maxwell–Pauli–Kohn–
Sham propagation with Octopus, we selected a nano-optical appli-
cation. We consider in this example two almost spherical sodium
nanoparticles with 297 sodium atoms each, which are arranged in
a dimer configuration. This system interacts with an incoming laser
pulse that excites either the internal dipole or quadrupole plasmon
motion of the dimer. In Fig. 2, we show the resulting electromag-
netic field enhancements for different levels of light–matter cou-
pling. In panels (a) and (b), we show the temporal profile of the
incoming laser pulse and the resulting current density at the cen-
ter point between the two nanoparticles. The field enhancement can
be seen in panels (c)–(e). Including a self-consistent back reaction
in the light–matter coupling, the field enhancement is reduced at
the center point between the two nanoparticles, as shown in panel
(c), while far away from the dimer, as shown in panels (d) and (e),
the field enhancement is larger than in the forward coupled case.
Furthermore, frequency shifts can also be observed, which are more
pronounced in the near-field. We found that the field enhancement
is also sensitive to the coupling terms of the multipole expansion,

<figure><img src="image_13.png" /><figcaption>FIG. 1. (Left) The standard forward-coupling approximation: the electromagnetic fields (in blue) propagate freely and only influence the propagation of the matter (in red). The back reaction of the matter currents on the electromagnetic fields is neglected. (Right) Illustration of a fully self-consistent predictor-corrector scheme for a coupled Maxwell–Pauli–Kohn–Sham time step. As before, the electromagnetic fields influence the propagation of the matter (forward coupling). However, in addition, here the currents from the matter propagation also influence the propagation of the electromagnetic fields (backward coupling). A given time step for the matter wavefunctions and the electromagnetic fields is repeated until self-consistency is found (self-consistent forward–backward coupling).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG. 2. Electric field and current density in the z-direction for a dimer of sodium nanoparticles. The centers of the two nanoparticles are located along the z-axis, and the distance between the two effective spheres of the nanoparticles is 0.5 nm. The electric field values are calculated at two different points: the center point $r_{cp}$ between the two nanoparticles, located at the origin, and an off-center point $r_{ocpx}$ located along the x-axis at a distance of 1.957 nm from $r_{cp}$. Panel (a) illustrates the incident cosinusoidal laser pulse with frequency $\omega_1 = 3.05$ eV (0.112 a.u.), $\lambda_1 = 406.5$ nm (7681.84 a.u.), and amplitude of $E_z^0 = 5.142 \times 10^7$ V/m ($10^{-4}$ a.u.), which drives the system. Panel (b) displays the current density at $r_{cp}$. In panels (c)–(e), we show the electric field enhancement at $r_{cp}$, the electric field enhancement at $r_{ocpx}$, and the average of the electric field over the detector surface close to the box boundary, respectively. The curve in bright gray in panel (c) was added to simplify the comparison and is identical to the laser pulse in panel (a). The period $T_1 = 1.36$ fs corresponding to the laser frequency $\omega_1$ is indicated with gray vertical lines.</figcaption></figure>

which are included, and that the quantitative difference of switching from local-density approximation (LDA) to the Perdew-Burke-Ernzerhof (PBE) functional for the exchange-correlation functional is in this case smaller than including the back reaction in the light-matter coupling (cf. Ref. 21).

To conclude, with our new efficient implementation for coupled Ehrenfest-Maxwell-Pauli-Kohn-Sham equations in Octopus, we have now a very versatile tool that allows us to compute fully self-consistent forward-backward light-matter coupling in real-time and real-space for a vast set of applications and in close analogy to experiments. As the coupling to Maxwell's equations of motion (EOM) represents the classical limit of the light-matter interaction, this development leads to the classical limit of QEDFT.

### III. STRONG ELECTRON-PHOTON INTERACTIONS IN REAL SPACE: QUANTUM-ELECTRODYNAMICAL DENSITY-FUNCTIONAL THEORY

The nascent field of strong light-matter interaction expanded over the past decades from small atomic structures to chemistry³ and materials science.²⁸ This development necessitates predictive first-principles methods capable to describe light and matter on the same footing. We introduced for the first time, a time-dependent density functional theory for quantum electrodynamics (QEDFT)⁴,¹⁰ to treat ab initio weak and strong light-matter interactions and its applications to chemistry and materials, which provides a unique framework to explore, predict, and control new states of matter out of equilibrium. This generalization of the time-dependent

density-functional method allows us, for the first time, to explore the effects of dressing electronic states with photons while retaining the electronic properties of real materials.

A general non-relativistic Hamiltonian for light-matter systems treating N interacting electrons coupled to $N_p$ photonic modes in the case of the so-called length-gauge, and after employing the long-wavelength (dipole) approximation,¹³ reads as follows:

$$ \hat{H} = \sum_{k=1}^{N} \left[ -\frac{1}{2} \partial_{\tau}^2 + v(\mathbf{r}_k, t) \right] + \frac{1}{2} \sum_{\alpha \neq l} w(\mathbf{r}_k, \mathbf{r}_l) + \frac{1}{2} \sum_{\alpha=1}^{N_p} \left[ -\frac{\partial^2}{\partial q_{\alpha}^2} + \left( \omega_{\alpha} q_{\alpha} - \lambda_{\alpha} \sum_{k=1}^{N} \mathbf{r}_k \right)^2 + 2 \frac{j_{\text{ext}}^{(\alpha)}(t)}{\omega_{\alpha}} q_{\alpha} \right]. \quad (8) $$

Here, the first two terms on the right-hand side correspond to the usual electronic many-body Hamiltonian, while the last term describes the photon modes, which is characterized for each photon mode $\alpha$ by its elongation $q_{\alpha}$, frequency $\omega_{\alpha}$, and electron-photon coupling strength vector $\lambda_{\alpha}$ that includes the polarization of the photon mode and introduces the coupling to the total dipole $\sum_{k=1}^{N} \mathbf{r}_k$ of the electronic system. For a more detailed discussion of the different regimes of weak/strong/ultra-strong light-matter coupling, we refer to Ref. 6. The external variable for the photon system is the time-dependent current $j_{\text{ext}}^{(\alpha)}(t)$.

QEDFT⁴,¹⁰ is structurally similar to time-dependent density-functional theory in that it is based on a one-to-one correspondence between internal and external variables. If now photons are also considered, the set of internal variables has to be expanded. In the
<figure><img src="image_14.png" /><figcaption>FIG. 2. Electric field and current density in the z-direction for a dimer of sodium nanoparticles. The centers of the two nanoparticles are located along the z-axis, and the distance between the two effective spheres of the nanoparticles is 0.5 nm. The electric field values are calculated at two different points: the center point $r_{cp}$ between the two nanoparticles, located at the origin, and an off-center point $r_{ocpx}$ located along the x-axis at a distance of 1.957 nm from $r_{cp}$. Panel (a) illustrates the incident cosinusoidal laser pulse with frequency $\omega_1 = 3.05$ eV (0.112 a.u.), $\lambda_1 = 406.5$ nm (7681.84 a.u.), and amplitude of $E_z^0 = 5.142 \times 10^7$ V/m ($10^{-4}$ a.u.), which drives the system. Panel (b) displays the current density at $r_{cp}$. In panels (c)–(e), we show the electric field enhancement at $r_{cp}$, the electric field enhancement at $r_{ocpx}$, and the average of the electric field over the detector surface close to the box boundary, respectively. The curve in bright gray in panel (c) was added to simplify the comparison and is identical to the laser pulse in panel (a). The period $T_1 = 1.36$ fs corresponding to the laser frequency $\omega_1$ is indicated with gray vertical lines.</figcaption></figure>

which are included, and that the quantitative difference of switching from local-density approximation (LDA) to the Perdew-Burke-Ernzerhof (PBE) functional for the exchange-correlation functional is in this case smaller than including the back reaction in the light-matter coupling (cf. Ref. 21).

To conclude, with our new efficient implementation for coupled Ehrenfest-Maxwell-Pauli-Kohn-Sham equations in Octopus, we have now a very versatile tool that allows us to compute fully self-consistent forward-backward light-matter coupling in real-time and real-space for a vast set of applications and in close analogy to experiments. As the coupling to Maxwell's equations of motion (EOM) represents the classical limit of the light-matter interaction, this development leads to the classical limit of QEDFT.

### III. STRONG ELECTRON-PHOTON INTERACTIONS IN REAL SPACE: QUANTUM-ELECTRODYNAMICAL DENSITY-FUNCTIONAL THEORY

The nascent field of strong light-matter interaction expanded over the past decades from small atomic structures to chemistry³ and materials science.²⁸ This development necessitates predictive first-principles methods capable to describe light and matter on the same footing. We introduced for the first time, a time-dependent density functional theory for quantum electrodynamics (QEDFT)⁴,¹⁰ to treat ab initio weak and strong light-matter interactions and its applications to chemistry and materials, which provides a unique framework to explore, predict, and control new states of matter out of equilibrium. This generalization of the time-dependent

density-functional method allows us, for the first time, to explore the effects of dressing electronic states with photons while retaining the electronic properties of real materials.

A general non-relativistic Hamiltonian for light-matter systems treating N interacting electrons coupled to $N_p$ photonic modes in the case of the so-called length-gauge, and after employing the long-wavelength (dipole) approximation,¹³ reads as follows:

$$ \hat{H} = \sum_{k=1}^{N} \left[ -\frac{1}{2} \partial_{\mathbf{r}_k}^2 + v(\mathbf{r}_k, t) \right] + \frac{1}{2} \sum_{\alpha \neq l} w(\mathbf{r}_k, \mathbf{r}_l) + \frac{1}{2} \sum_{\alpha=1}^{N_p} \left[ -\frac{\partial^2}{\partial q_{\alpha}^2} + \left( \omega_{\alpha} q_{\alpha} - \lambda_{\alpha} \sum_{k=1}^{N} \mathbf{r}_k \right)^2 + 2 \frac{j_{\text{ext}}^{(\alpha)}(t)}{\omega_{\alpha}} q_{\alpha} \right]. \quad (8) $$

Here, the first two terms on the right-hand side correspond to the usual electronic many-body Hamiltonian, while the last term describes the photon modes, which is characterized for each photon mode $\alpha$ by its elongation $q_{\alpha}$, frequency $\omega_{\alpha}$, and electron-photon coupling strength vector $\lambda_{\alpha}$ that includes the polarization of the photon mode and introduces the coupling to the total dipole $\sum_{k=1}^{N} \mathbf{r}_k$ of the electronic system. For a more detailed discussion of the different regimes of weak/strong/ultra-strong light-matter coupling, we refer to Ref. 6. The external variable for the photon system is the time-dependent current $j_{\text{ext}}^{(\alpha)}(t)$.

QEDFT⁴,¹⁰ is structurally similar to time-dependent density-functional theory in that it is based on a one-to-one correspondence between internal and external variables. If now photons are also considered, the set of internal variables has to be expanded. In the
framework of Eq. (8), the internal variables become the density $n(\mathbf{r}, t)$ and the mode-resolved contributions to the electric displacement field $q_\alpha(t)$. By introducing and exploiting the bijective mapping of these internal and the external variables [$v_{\text{ext}}(\mathbf{r}), t$] and $j_{\text{ext}}(t)$], the auxiliary Kohn-Sham system is described by the electronic Kohn-Sham equations as well as Maxwell's equations,⁵ leading to no exchange-correlation contribution in the photon subsystem [$j_{xc}^{(\alpha)}(t) = 0$]. This reformulation subsumes the “quantumness” of the light-matter interaction solely into the local exchange-correlation potential that now features a component due to the electron-photon interactions, in addition to the part due to electron-electron interaction, i.e., $v_{xc\sigma}(\mathbf{r}, t) = v_{xc\sigma}^{ee}(\mathbf{r}, t) + v_{xc\sigma}^{ep}(\mathbf{r}, t)$.⁴,¹⁰,¹² Extending the coupled Maxwell-Kohn-Sham equations to consider quantum photons is, thus, condensed into the calculation of the local exchange-correlation potential.

In practice, QEDFT requires the construction of an additional exchange-correlation potential that describes the electron-photon interaction. First attempts to construct this potential were based on many-body perturbation theory via the Sham-Schlüter equation²⁹ within the exact exchange (EXX) approximation¹³,³⁰ by utilizing the optimized-effective potential (OEP) method.³¹

### A. The optimized effective potential

The framework of the optimized effective potential (OEP) equation for electron-photon systems has been introduced in Ref. 30 for the static and time-dependent regimes. Solving the time-dependent OEP equation is unfortunately computationally challenging. Applications to electronic systems³²⁻³⁴ have been realized so far only for small and low-dimensional systems while their utilization for interacting electron-photon systems remained restricted to models.³⁰,³⁵,³⁶ This section provides a brief introduction to the implementation that allows us to solve the computationally more tractable static OEP equation.

As introduced in Ref. 30, the OEP photon energy depends on both the occupied and unoccupied electronic orbitals. Alternatively, the OEP photon energy can also be formulated using occupied orbitals and orbital shifts only.¹³ The full energy expression is given by

$$E_{xc} = E_{xc}^{(ee)} + \sum_{\alpha=1}^{N_p} E_{x}^{(\alpha)}, \qquad (9)$$

where $E_{xc}^{(ee)}$ describes the electronic exchange-correlation energy and $E_{x}^{(\alpha)}$ describes the exchange energy due to the interaction of the electrons with the photon mode $\alpha$. Avoiding unoccupied orbitals is computationally much more favorable for larger systems and the photonic induced exchange energy can be correspondingly expressed as a sum over the $N_\sigma$ occupied orbitals of spin channel $\sigma$,

$$E_{x}^{(\alpha)} = \sum_{\sigma=\uparrow,\downarrow} \sum_{i=1}^{N_{\sigma}} \sqrt{\frac{\omega_{\alpha}}{8}} \langle \Phi_{i\sigma,\alpha}^{(1)} | \hat{d}_{\alpha} | \varphi_{i\sigma} \rangle + \frac{1}{4} \langle \Phi_{i\sigma,\alpha}^{(2)} | \hat{d}_{\alpha} | \varphi_{i\sigma} \rangle + \text{c.c.}, \quad (10)$$

where $\omega_\alpha$ describes the $\alpha$th mode of the electromagnetic field and $\hat{d}_\alpha = \lambda_\alpha \cdot \mathbf{r}$ describes the dipole operator and the electron-photon coupling strength. We can now reformulate the problem in terms of two

electron-photon orbital shifts. The Kohn-Sham orbitals $\varphi_{i\sigma}$ contribute to both electron-photon orbital shifts $\Phi_{i\sigma,\alpha}^{(1)}$ and $\Phi_{i\sigma,\alpha}^{(2)}$ that can be calculated using the Sternheimer equations.¹³ The first electron-photon orbital shift can be obtained explicitly by the solution of a linear Sternheimer equation

$$[\hat{h}_{s\sigma} - (\epsilon_{i\sigma} - \omega_{\alpha})] \Phi_{i\sigma,\alpha}^{(1)}(\mathbf{r}) = -\sqrt{\frac{\omega_{\alpha}}{2}} \hat{d}_{\alpha} \varphi_{i\sigma}(\mathbf{r}) + \sqrt{\frac{\omega_{\alpha}}{2}} \sum_{k=1}^{N_\sigma} d_{ki\sigma}^{(\alpha)} \varphi_{k\sigma}(\mathbf{r}) \quad (11)$$

with the matrix element $d_{ij\sigma}^{(\alpha)} = \langle \varphi_{i\sigma} | \hat{d}_{\alpha} | \varphi_{j\sigma} \rangle$. In addition, the second electron-photon orbital shift $\Phi_{i\sigma,\alpha}^{(2)}(\mathbf{r})$ can be defined explicitly as follows:

$$\Phi_{i\sigma,\alpha}^{(2)}(\mathbf{r}) = \hat{d}_{\alpha} \varphi_{i\sigma}(\mathbf{r}) - \sum_{k=1}^{N_{\sigma}} d_{ki\sigma}^{(\alpha)} \varphi_{k\sigma}(\mathbf{r}). \quad (12)$$

From Eq. (10), we can now deduce the potential using

$$v_{xc\sigma}(\mathbf{r}) = \frac{\delta E_{xc}}{\delta n_{\sigma}(\mathbf{r})}. \quad (13)$$

After these reformulations, we find for the final OEP equation including electron-electron effects as well as electron-photon effects

$$\sum_{i=1}^{N_{\sigma}} \psi_{i\sigma}^{*}(\mathbf{r}) \varphi_{i\sigma}(\mathbf{r}) - \Lambda_{i\sigma}(\mathbf{r}) + \text{c.c.} = 0, \quad (14)$$

where the inhomogeneity $\Lambda_{i\sigma}(\mathbf{r})$ is given by

$$\Lambda_{i\sigma}(\mathbf{r}) = \frac{1}{2} \sum_{\alpha=1}^{N_p} \left[ \Phi_{i\sigma,\alpha}^{(1)*}(\mathbf{r}) \Phi_{i\sigma,\alpha}^{(1)}(\mathbf{r}) - (\Phi_{i\sigma,\alpha}^{(1)}|\Phi_{i\sigma,\alpha}^{(1)}) \varphi_{i\sigma}^{*}(\mathbf{r}) \varphi_{i\sigma}(\mathbf{r}) \right].$$

In Eq. (14), we defined a third orbital shift, the exchange-correlation orbital shift, which will be used to obtain the corresponding exchange-correlation potential and is defined along the lines of the orbital shift usually used in OEP calculations.³¹ We can also obtain $\psi_{i\sigma}^{*}(\mathbf{r})$ using a Sternheimer equation

$$(\hat{h}_{s\sigma} - \epsilon_{i\sigma}) \psi_{i\sigma}^{*}(\mathbf{r}) = M_{i\sigma}^{*}(\mathbf{r}) - \langle M_{i\sigma} | \varphi_{i\sigma} \rangle \varphi_{i\sigma}^{*}(\mathbf{r}), \quad (15)$$

where $M_{i\sigma}^{*}(\mathbf{r})$ now consist of the electron-photon orbital shifts and the Kohn-Sham orbitals, as described in Ref. 13. Accordingly, we define this quantity as

$$M_{i\sigma}^{*}(\mathbf{r}) = -[v_{x\sigma}(\mathbf{r}) - u_{x\sigma}(\mathbf{r})] \varphi_{i\sigma}^{*}(\mathbf{r}) + \sum_{\alpha=1}^{N_p} \left[ \hat{d}_{\alpha} \left( \sqrt{\frac{\omega_{\alpha}}{2}} \Phi_{i\sigma,\alpha}^{(1)*}(\mathbf{r}) + \frac{1}{2} \hat{d}_{\alpha} \varphi_{i\sigma}^{*}(\mathbf{r}) \right) - \sum_{k=1}^{N_{\sigma}} d_{ik\sigma}^{(\alpha)} \left( \sqrt{\frac{\omega_{\alpha}}{2}} \Phi_{k\sigma,\alpha}^{(1)*}(\mathbf{r}) + \hat{d}_{\alpha} \varphi_{k\sigma}^{*}(\mathbf{r}) \right) \right], \quad (16)$$

where the effects of the electron-electron interaction are included in the quantity $u_{x\sigma}(\mathbf{r})$. For instance, in exchange-only calculations, this quantity is defined as $u_{x\sigma}(\mathbf{r}) = \frac{1}{\varphi_{i\sigma}^{*}(\mathbf{r})} \frac{\delta E_{x}^{(ee)}}{\delta \varphi_{i\sigma}(\mathbf{r})}$, where $E_{x}^{(ee)}$ is the usual Fock exchange energy.

Equation (15) has to be solved self-consistently with Eq. (11). By this reformulation, we replaced the problem of calculating the
framework of Eq. (8), the internal variables become the density $n(\mathbf{r}, t)$ and the mode-resolved contributions to the electric displacement field $q_\alpha(t)$. By introducing and exploiting the bijective mapping of these internal and the external variables [$v_{\text{ext}}(\mathbf{r}), t)$ and $j_{\text{ext}}(t)$], the auxiliary Kohn-Sham system is described by the electronic Kohn-Sham equations as well as Maxwell's equations,⁵ leading to no exchange-correlation contribution in the photon subsystem [$j_{xc}^{(\alpha)}(t) = 0$]. This reformulation subsumes the “quantumness” of the light-matter interaction solely into the local exchange-correlation potential that now features a component due to the electron-photon interactions, in addition to the part due to electron-electron interaction, i.e., $v_{xc\sigma}(\mathbf{r}, t) = v_{xc\sigma}^{ee}(\mathbf{r}, t) + v_{xc\sigma}^{ep}(\mathbf{r}, t)$.⁴,¹⁰,¹² Extending the coupled Maxwell-Kohn-Sham equations to consider quantum photons is, thus, condensed into the calculation of the local exchange-correlation potential.

In practice, QEDFT requires the construction of an additional exchange-correlation potential that describes the electron-photon interaction. First attempts to construct this potential were based on many-body perturbation theory via the Sham-Schlüter equation²⁹ within the exact exchange (EXX) approximation¹³,³⁰ by utilizing the optimized-effective potential (OEP) method.³¹

## A. The optimized effective potential

The framework of the optimized effective potential (OEP) equation for electron-photon systems has been introduced in Ref. 30 for the static and time-dependent regimes. Solving the time-dependent OEP equation is unfortunately computationally challenging. Applications to electronic systems³²⁻³⁴ have been realized so far only for small and low-dimensional systems while their utilization for interacting electron-photon systems remained restricted to models.³⁰,³⁵,³⁶ This section provides a brief introduction to the implementation that allows us to solve the computationally more tractable static OEP equation.

As introduced in Ref. 30, the OEP photon energy depends on both the occupied and unoccupied electronic orbitals. Alternatively, the OEP photon energy can also be formulated using occupied orbitals and orbital shifts only.¹³ The full energy expression is given by

$$E_{xc} = E_{xc}^{(ee)} + \sum_{\alpha=1}^{N_p} E_{x}^{(\alpha)}, \quad (9)$$

where $E_{xc}^{(ee)}$ describes the electronic exchange-correlation energy and $E_{x}^{(\alpha)}$ describes the exchange energy due to the interaction of the electrons with the photon mode $\alpha$. Avoiding unoccupied orbitals is computationally much more favorable for larger systems and the photonic induced exchange energy can be correspondingly expressed as a sum over the $N_\sigma$ occupied orbitals of spin channel $\sigma$,

$$E_{x}^{(\alpha)} = \sum_{\sigma=\uparrow,\downarrow} \sum_{i=1}^{N_{\sigma}} \sqrt{\frac{\omega_{\alpha}}{8}} \langle \Phi_{i\sigma,\alpha}^{(1)} | \hat{d}_{\alpha} | \varphi_{i\sigma} \rangle + \frac{1}{4} \langle \Phi_{i\sigma,\alpha}^{(2)} | \hat{d}_{\alpha} | \varphi_{i\sigma} \rangle + \text{c.c.}, \quad (10)$$

where $\omega_\alpha$ describes the $\alpha$th mode of the electromagnetic field and $\hat{d}_\alpha = \lambda_\alpha \cdot \mathbf{r}$ describes the dipole operator and the electron-photon coupling strength. We can now reformulate the problem in terms of two

electron-photon orbital shifts. The Kohn-Sham orbitals $\varphi_{i\sigma}$ contribute to both electron-photon orbital shifts $\Phi_{i\sigma,\alpha}^{(1)}$ and $\Phi_{i\sigma,\alpha}^{(2)}$ that can be calculated using the Sternheimer equations.¹³ The first electron-photon orbital shift can be obtained explicitly by the solution of a linear Sternheimer equation

$$[\hat{h}_{s\sigma} - (\epsilon_{i\sigma} - \omega_{\alpha})] \Phi_{i\sigma,\alpha}^{(1)}(\mathbf{r}) = -\sqrt{\frac{\omega_{\alpha}}{2}} \hat{d}_{\alpha} \varphi_{i\sigma}(\mathbf{r}) + \sqrt{\frac{\omega_{\alpha}}{2}} \sum_{k=1}^{N_\sigma} d_{ki\sigma}^{(\alpha)} \varphi_{k\sigma}(\mathbf{r}) \quad (11)$$

with the matrix element $d_{ij\sigma}^{(\alpha)} = \langle \varphi_{i\sigma} | \hat{d}_{\alpha} | \varphi_{j\sigma} \rangle$. In addition, the second electron-photon orbital shift $\Phi_{i\sigma,\alpha}^{(2)}(\mathbf{r})$ can be defined explicitly as follows:

$$\Phi_{i\sigma,\alpha}^{(2)}(\mathbf{r}) = \hat{d}_{\alpha} \varphi_{i\sigma}(\mathbf{r}) - \sum_{k=1}^{N_\sigma} d_{ki\sigma}^{(\alpha)} \varphi_{k\sigma}(\mathbf{r}). \quad (12)$$

From Eq. (10), we can now deduce the potential using

$$v_{xc\sigma}(\mathbf{r}) = \frac{\delta E_{xc}}{\delta n_{\sigma}(\mathbf{r})}. \quad (13)$$

After these reformulations, we find for the final OEP equation including electron-electron effects as well as electron-photon effects

$$\sum_{i=1}^{N_\sigma} \psi_{i\sigma}^{*}(\mathbf{r}) \varphi_{i\sigma}(\mathbf{r}) - \Lambda_{i\sigma}(\mathbf{r}) + \text{c.c.} = 0, \quad (14)$$

where the inhomogeneity $\Lambda_{i\sigma}(\mathbf{r})$ is given by

$$\Lambda_{i\sigma}(\mathbf{r}) = \frac{1}{2} \sum_{\alpha=1}^{N_p} \left[ \Phi_{i\sigma,\alpha}^{(1)*}(\mathbf{r}) \Phi_{i\sigma,\alpha}^{(1)}(\mathbf{r}) - (\Phi_{i\sigma,\alpha}^{(1)}|\Phi_{i\sigma,\alpha}^{(1)}) \psi_{i\sigma}^{*}(\mathbf{r}) \varphi_{i\sigma}(\mathbf{r}) \right].$$

In Eq. (14), we defined a third orbital shift, the exchange-correlation orbital shift, which will be used to obtain the corresponding exchange-correlation potential and is defined along the lines of the orbital shift usually used in OEP calculations.³¹ We can also obtain $\psi_{i\sigma}^{*}(\mathbf{r})$ using a Sternheimer equation

$$(\hat{h}_{s\sigma} - \epsilon_{i\sigma}) \psi_{i\sigma}^{*}(\mathbf{r}) = M_{i\sigma}^{*}(\mathbf{r}) - \langle M_{i\sigma} | \varphi_{i\sigma} \rangle \psi_{i\sigma}^{*}(\mathbf{r}), \quad (15)$$

where $M_{i\sigma}^{*}(\mathbf{r})$ now consist of the electron-photon orbital shifts and the Kohn-Sham orbitals, as described in Ref. 13. Accordingly, we define this quantity as

$$M_{i\sigma}^{*}(\mathbf{r}) = -[v_{x\sigma}(\mathbf{r}) - u_{x\sigma}(\mathbf{r})] \varphi_{i\sigma}^{*}(\mathbf{r}) + \sum_{\alpha=1}^{N_p} \left[ \hat{d}_{\alpha} \left( \sqrt{\frac{\omega_{\alpha}}{2}} \Phi_{i\sigma,\alpha}^{(1)*}(\mathbf{r}) + \frac{1}{2} \hat{d}_{\alpha} \varphi_{i\sigma}^{*}(\mathbf{r}) \right) - \sum_{k=1}^{N_\sigma} d_{ik\sigma}^{(\alpha)} \left( \sqrt{\frac{\omega_{\alpha}}{2}} \Phi_{k\sigma,\alpha}^{(1)*}(\mathbf{r}) + \hat{d}_{\alpha} \varphi_{k\sigma}^{*}(\mathbf{r}) \right) \right], \quad (16)$$

where the effects of the electron-electron interaction are included in the quantity $u_{x\sigma}(\mathbf{r})$. For instance, in exchange-only calculations, this quantity is defined as $u_{x\sigma}(\mathbf{r}) = \frac{1}{\varphi_{i\sigma}^{*}(\mathbf{r})} \frac{\delta E_{x}^{(ee)}}{\delta \varphi_{i\sigma}(\mathbf{r})}$, where $E_{x}^{(ee)}$ is the usual Fock exchange energy.

Equation (15) has to be solved self-consistently with Eq. (11). By this reformulation, we replaced the problem of calculating the
OEP equation using all unoccupied states by a problem of solving $N_p+1$ Sternheimer equations that only involve the occupied orbitals. In this way, the formulation of the problem becomes similar to the one of Ref. 37 for electrons only.

For practical implementation, we reformulate the OEP equation in the following form, as it is commonly done to construct the electronic OEP:³¹

$$S_{\sigma}(\mathbf{r}) = \sum_{i=1}^{N_{\sigma}} \psi_{i\sigma}^{*}(\mathbf{r})\varphi_{i\sigma}(\mathbf{r}) - \Lambda_{i\sigma}(\mathbf{r}) + \text{c.c.}, \quad (17)$$

and update the potential with

$$v_{x\sigma}^{(\text{new})}(\mathbf{r}) = v_{x\sigma}^{(\text{old})}(\mathbf{r}) + c(\mathbf{r})S_{\sigma}(\mathbf{r}). \quad (18)$$

The quantity $S_σ(\mathbf{r})$ becomes a measure for convergence, since it is vanishing at the solution point [compare Eqs. (17) and (14)]. For the function $c(\mathbf{r})$, we have different possibilities, such as using a constant or using the inverse of the electron density, as used in Ref. 37. Other methods are also possible, such as the Barzilai-Borwein (BB) method.³⁸ We found stable algorithms when using a constant and the Barzilai-Borwein method.

While we will show the computational feasibility of this approach in the following, often a simplified solution is beneficial as a starting point for the self-consistency procedure. Such a simplified approximation can be deduced by reformulating Eq. (15) into the following equivalent form:

$$v_{x\sigma}(\mathbf{r}) = \frac{1}{2n_{\sigma}(\mathbf{r})} \sum_{i=1}^{N_{\sigma}} \left( \langle \varphi_{i\sigma} | v_{x\sigma} | \varphi_{i\sigma} \rangle |\varphi_{i\sigma}(\mathbf{r})|^2 \right. \\ \left. + \left[ \{M_{i\sigma}^*(\mathbf{r}) + v_{x\sigma}(\mathbf{r})\varphi_{i\sigma}^*(\mathbf{r})\} \right. \right. \\ \left. \left. - \{\langle M_{i\sigma} + v_{x\sigma}\varphi_{i\sigma} \rangle |\varphi_{i\sigma}\}\varphi_{i\sigma}^*(\mathbf{r}) \right. \right. \\ \left. \left. - (\hat{h}_{s\sigma}(\mathbf{r}) - \epsilon_{i\sigma})\psi_{i\sigma}^*(\mathbf{r})\right] \varphi_{i\sigma}(\mathbf{r}) \right) + \text{c.c.}, \quad (19)$$

and subsequently assuming $(\hat{h}_{s\sigma}(\mathbf{r}) - \epsilon_{i\sigma})\psi_{i\sigma}^*(\mathbf{r}) = 0$ to start the iterative process. In situations where $\Lambda_{i\sigma}(\mathbf{r}) = 0$, such as pure

electronic exact exchange, this approximation is exact for a single electron and is referred to as the Krieger-Li-Iafrate (KLI) approximation.¹³,³⁷,³⁹⁻⁴¹ By multiplying Eq. (19) by $|\varphi_{j\sigma}(\mathbf{r})|^2$ and integrating over the spatial coordinates, we arrive at a linear equation that, in turn, can be solved for the approximate $v_{x\sigma}$. This approximation scheme has proven to often deliver sufficiently accurate results for electronic structure calculations with significantly lower computational effort and reliable stability. As this approximation can be seen as a diagonal approximation to the response function, it unavoidably fails in accurately describing polarization features. In the context of light-matter correlated ground-states, this leads to a slight unbalance when approximating components including photonic excitations [introduced by $\Phi^{(1)}$] and self-polarization interaction [introduced by $\Phi^{(2)}$].¹³,³⁰,⁴²,⁴³ This results in a violation of translational invariance and introduces an artificial dependence on the permanent dipole. When performing KLI calculations including light-matter interaction, we thus suggest moving the set of coordinates into the electronic center-of-charge instead of the center-of-mass frame. To reduce the effect of this dependence during a self-consistent calculation, the optional input parameter KLIPt_coc has been introduced in the code. When activated, this option defines the dipole operator $\hat{d}_\alpha$ with respect to the electronic center-of-charge, thus improving the stability of the algorithm.

Finally, in Fig. 3, we show the capabilities of the new implementation, where we calculate two sodium dimers in the weak coupling regime under light-matter coupling. Figure 3(a) shows the electron density, and Fig. 3(b) shows the comparison of OEP and KLI results. As shown in Ref. 13, in the weak-coupling regime, the KLI is close to the OEP result. In Fig. 3(c), we show the convergence behavior when using a constant in Eq. (18) and when using the Barzilai-Borwein method.³⁸

For a more detailed analysis on how the OEP equation performs also in the strong and the ultra-strong coupling limit, we refer to Refs. 13 and 30, where we find accurate results for few-level systems and asymmetric external potentials up to high electron-photon coupling strength.

<figure><img src="image_13.png" /><figcaption>FIG. 3. Ground-state density of two sodium dimers affected by the vacuum-field of a cavity: In (a), we show the electron density, in (b), we show the difference of the electron density inside and outside the cavity for the OEP and the KLI approximations (see main text for details). We find for the exchange energy $E_x^{(\alpha)} = 6.52$ meV and the number of photons in the correlated light-matter ground-state $n_{pt} = \langle a_\alpha^\dagger a_\alpha \rangle = 2.03 \times 10^{-3}$ for the OEP case, and $E_x^{(\alpha)} = 6.67$ meV and $n_{pt} = 2.22 \times 10^{-3}$ for KLI. Panel (c) shows the convergence of $|S(\mathbf{r})|^2$ as defined by Eq. (17) using a constant with $c = 20$ and using the Barzilai-Borwein (BB) scheme. The simulation is set up as described in Refs. 12 and 13, with $\hbar\omega_\alpha = 2.19$ eV, $\lambda_\alpha = 2.95$ eV$^{1/2}$ nm$^{-1}$, and $\lambda_\alpha = \lambda_\alpha e_y$, and the real-space grid is sampled as 31.75 × 31.75 × 31.75 Å³ with a grid spacing of 0.265 Å.</figcaption></figure>
OEP equation using all unoccupied states by a problem of solving $N_p+1$ Sternheimer equations that only involve the occupied orbitals. In this way, the formulation of the problem becomes similar to the one of Ref. 37 for electrons only.

For practical implementation, we reformulate the OEP equation in the following form, as it is commonly done to construct the electronic OEP:³¹

$$S_{\sigma}(\mathbf{r}) = \sum_{i=1}^{N_{\sigma}} \psi_{i\sigma}^{*}(\mathbf{r})\varphi_{i\sigma}(\mathbf{r}) - \Lambda_{i\sigma}(\mathbf{r}) + \text{c.c.}, \quad (17)$$

and update the potential with

$$v_{x\sigma}^{(\text{new})}(\mathbf{r}) = v_{x\sigma}^{(\text{old})}(\mathbf{r}) + c(\mathbf{r})S_{\sigma}(\mathbf{r}). \quad (18)$$

The quantity $S_\sigma(\mathbf{r})$ becomes a measure for convergence, since it is vanishing at the solution point [compare Eqs. (17) and (14)]. For the function $c(\mathbf{r})$, we have different possibilities, such as using a constant or using the inverse of the electron density, as used in Ref. 37. Other methods are also possible, such as the Barzilai-Borwein (BB) method.³⁸ We found stable algorithms when using a constant and the Barzilai-Borwein method.

While we will show the computational feasibility of this approach in the following, often a simplified solution is beneficial as a starting point for the self-consistency procedure. Such a simplified approximation can be deduced by reformulating Eq. (15) into the following equivalent form:

$$v_{x\sigma}(\mathbf{r}) = \frac{1}{2n_{\sigma}(\mathbf{r})} \sum_{i=1}^{N_{\sigma}} \left( \langle \varphi_{i\sigma} | v_{x\sigma} | \varphi_{i\sigma} \rangle |\varphi_{i\sigma}(\mathbf{r})|^2 \right. \\ \left. + \left[ \{M_{i\sigma}^*(\mathbf{r}) + v_{x\sigma}(\mathbf{r})\varphi_{i\sigma}^*(\mathbf{r})\} \right. \right. \\ \left. \left. - \{\langle M_{i\sigma} + v_{x\sigma}\varphi_{i\sigma} \rangle |\varphi_{i\sigma}\}\varphi_{i\sigma}^*(\mathbf{r}) \right. \right. \\ \left. \left. - (\hat{h}_{s\sigma}(\mathbf{r}) - \epsilon_{i\sigma})\psi_{i\sigma}^*(\mathbf{r})\right] \varphi_{i\sigma}(\mathbf{r}) \right) + \text{c.c.}, \quad (19)$$

and subsequently assuming $(\hat{h}_{s\sigma}(\mathbf{r}) - \epsilon_{i\sigma})\psi_{i\sigma}^*(\mathbf{r}) = 0$ to start the iterative process. In situations where $\Lambda_{i\sigma}(\mathbf{r}) = 0$, such as pure

electronic exact exchange, this approximation is exact for a single electron and is referred to as the Krieger-Li-Iafrate (KLI) approximation.¹³,³⁷,³⁹⁻⁴¹ By multiplying Eq. (19) by $|\varphi_{j\sigma}(\mathbf{r})|^2$ and integrating over the spatial coordinates, we arrive at a linear equation that, in turn, can be solved for the approximate $v_{x\sigma}$. This approximation scheme has proven to often deliver sufficiently accurate results for electronic structure calculations with significantly lower computational effort and reliable stability. As this approximation can be seen as a diagonal approximation to the response function, it unavoidably fails in accurately describing polarization features. In the context of light-matter correlated ground-states, this leads to a slight unbalance when approximating components including photonic excitations [introduced by $\Phi^{(1)}$] and self-polarization interaction [introduced by $\Phi^{(2)}$].¹³,³⁰,⁴²,⁴³ This results in a violation of translational invariance and introduces an artificial dependence on the permanent dipole. When performing KLI calculations including light-matter interaction, we thus suggest moving the set of coordinates into the electronic center-of-charge instead of the center-of-mass frame. To reduce the effect of this dependence during a self-consistent calculation, the optional input parameter KLIPt_coc has been introduced in the code. When activated, this option defines the dipole operator $\hat{d}_\alpha$ with respect to the electronic center-of-charge, thus improving the stability of the algorithm.

Finally, in Fig. 3, we show the capabilities of the new implementation, where we calculate two sodium dimers in the weak coupling regime under light-matter coupling. Figure 3(a) shows the electron density, and Fig. 3(b) shows the comparison of OEP and KLI results. As shown in Ref. 13, in the weak-coupling regime, the KLI is close to the OEP result. In Fig. 3(c), we show the convergence behavior when using a constant in Eq. (18) and when using the Barzilai-Borwein method.³⁸

For a more detailed analysis on how the OEP equation performs also in the strong and the ultra-strong coupling limit, we refer to Refs. 13 and 30, where we find accurate results for few-level systems and asymmetric external potentials up to high electron-photon coupling strength.

<figure><img src="image_15.png" /><figcaption>FIG. 3. Ground-state density of two sodium dimers affected by the vacuum-field of a cavity: In (a), we show the electron density, in (b), we show the difference of the electron density inside and outside the cavity for the OEP and the KLI approximations (see main text for details). We find for the exchange energy $E_x^{(\alpha)} = 6.52$ meV and the number of photons in the correlated light-matter ground-state $n_{pt} = \langle a_\alpha^\dagger a_\alpha \rangle = 2.03 \times 10^{-3}$ for the OEP case, and $E_x^{(\alpha)} = 6.67$ meV and $n_{pt} = 2.22 \times 10^{-3}$ for KLI. Panel (c) shows the convergence of $|S(\mathbf{r})|^2$ as defined by Eq. (17) using a constant with $c = 20$ and using the Barzilai-Borwein (BB) scheme. The simulation is set up as described in Refs. 12 and 13, with $\hbar\omega_\alpha = 2.19$ eV, $\lambda_\alpha = 2.95$ eV$^{1/2}$ nm$^{-1}$, and $\lambda_\alpha = \lambda_\alpha e_y$, and the real-space grid is sampled as 31.75 × 31.75 × 31.75 Å³ with a grid spacing of 0.265 Å.</figcaption></figure>
Extensions of QEDFT to the regime of vibrational strong coupling,44 the linear-response regime,15 and multitrajectory methods that capture quantum fluctuations46 are currently work in progress and will further strengthen the capabilities of the Octopus code for the real-space description of strongly coupled light-matter systems. To describe the effects in the ultra-strong coupling regime, one can use an alternative method that is presented in Sec. IV.

## IV. DRESSED REDUCED DENSITY MATRIX FUNCTIONAL THEORY FOR ULTRA-STRONGLY COUPLED LIGHT-MATTER SYSTEMS

The accurate description of the (ultra-)strong coupling regime of light-matter systems is a formidable task. In many cases, the known functionals for QEDFT (see Sec. III) are inaccurate, and for complex electronic systems, typical few-level approximations become unreliable.42,47 In this section, we present the Octopus implementation of an alternative real-space ab initio method for coupled light-matter systems. Dressed Reduced Density Matrix Functional Theory (dressed RDMFT)48 extends standard electronic RDMFT to coupled light-matter systems similarly to how QEDFT extends DFT. First tests on simple model systems suggest that the dressed RDMFT remains accurate from the weak to the ultra-strong coupling regime. A proper introduction of the theory, examples, and convergence studies can be found in Ref. 49.

This approach allows for the description of an interacting N-electron system coupled to one photonic mode within the dipole approximation. The respective Hamiltonian is given in Eq. (8) of Sec. III. Note that we set $j_{ext} = 0$ throughout this section and that the ground state $\Psi$ of Eq. (8) depends on $4N + 1$ coordinates, i.e., $\Psi = \Psi(r_1, \sigma_1, \dots, r_N, \sigma_N, p)$, where $\{\sigma_i\}$ denotes the spin degrees of freedom and $p$ is the elongation of the photon mode.

Within the dressed RDMFT, the original Hamiltonian (8) of N electrons in d dimensions and one mode is replaced by an extended auxiliary Hamiltonian of N dressed fermions in d + 1 dimensions with coordinates $z = (r, q) \in \mathbb{R}^{d+1}$. This auxiliary Hamiltonian reads

$$ \hat{H}' = \sum_{k=1}^{N} \left[ -\frac{1}{2} \Delta'_{k} + v'(z_{k}) \right] + \frac{1}{2} \sum_{k \neq l} w'(z_{k}, z_{l}) \quad (20) $$

and gives access to the same physics (see Refs. 49, Sec. 4). For a d = 1 matter subsystem, the operators introduced in Eq. (20) read

as follows: the dressed Laplacian $\Delta' = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial q^2}$, the dressed local potential

$$ v'(z) = v(x) + \left[ \frac{1}{2} \omega^2 q^2 - \frac{\omega}{\sqrt{N}} \lambda q x + \frac{1}{2} (\lambda x)^2 \right], \quad (21) $$

and the dressed interaction kernel

$$ w'(z, z') = w(x, x') + \left[ -\frac{\omega}{\sqrt{N}} \lambda q x' - \frac{\omega}{\sqrt{N}} \lambda q' x + \lambda^2 x x' \right]. \quad (22) $$

The ground state $\Psi'(z_1, \sigma_1, \dots, z_N, \sigma_N) = \Psi(x_1, \sigma_1, \dots, x_N, \sigma_N) \otimes \chi(p_2, \dots, p_N)$ of $\hat{H}'$ is a product of the original physical ground state $\Psi$ and the ground state of $\chi$, which, in turn, is the product of $N-1$ harmonic oscillator ground states. The auxiliary Hamiltonian (20) contains only one-body and two-body terms in terms of the dressed coordinates, which allows us, in principle, to use any standard electronic structure method to solve it (see also Refs. 49, Secs. 4 and 5). We use this construction to develop the dressed RDMFT and dressed Hartree-Fock (HF). For that, we define the dressed (spin-summed) first order reduced density matrix (1RDM)

$$ y'(z, z') = N \sum_{\sigma_1, \dots, \sigma_N} \int d^{2(N-1)}z \Psi^{*}(z'_{1}, z_{2}, \dots, z_{N} \sigma_{N}) \times \Psi'(z \sigma_{1}, z_{2}, \dots, z_{N} \sigma_{N}). \quad (23) $$

To apply RDMFT theory on the auxiliary system, we have to replace the total energy functional of electronic RDMFT, given in Ref. 19, with the newly introduced quantities of the dressed system, i.e., the auxiliary Hamiltonian of Eq. (20) (approximately) evaluated by the dressed 1RDM $y'$ of Eq. (23). By that, common approximations for the two-body energy expression in terms of the 1RDM can be directly transferred from electronic theory to the dressed system.50 The minimization is performed like in the electronic case and is based on the RDMFT implementation of Octopus, though the convergence of the dressed system requires a more complicated protocol that can be found in the supplementary material of Ref. 49. The current implementation in Octopus approximates the conditions under which the dressed 1RDM corresponds to a wave function by ensuring the fermionic ensemble N-representability conditions.51 However, the auxiliary wave function also exhibits another exchange symmetry with respect to the auxiliary coordinates, which is currently neglected. For the practical validity of this approach, the reader is referred to Sec. 5 and the supplementary material of Ref. 49.

<figure><img src="image_14.png" /><figcaption>FIG. 4. Differences of the dressed HF (dHF) and the dressed RDMFT (dRDMFT) from the exact ground state energies (in hartree, left) and from the exact photon number (right) as a function of the coupling $g/\omega$ for the stretched H$_2$ molecule in the dressed orbital description. The dressed RDMFT improves considerably upon dressed HF for the energy, especially due to the better description of the electronic correlation. The photon number, an example of a photonic observable, is captured similarly with both methods.</figcaption></figure>
Extensions of QEDFT to the regime of vibrational strong coupling,44 the linear-response regime,15 and multitrajectory methods that capture quantum fluctuations46 are currently work in progress and will further strengthen the capabilities of the Octopus code for the real-space description of strongly coupled light-matter systems. To describe the effects in the ultra-strong coupling regime, one can use an alternative method that is presented in Sec. IV.

## IV. DRESSED REDUCED DENSITY MATRIX FUNCTIONAL THEORY FOR ULTRA-STRONGLY COUPLED LIGHT-MATTER SYSTEMS

The accurate description of the (ultra-)strong coupling regime of light-matter systems is a formidable task. In many cases, the known functionals for QEDFT (see Sec. III) are inaccurate, and for complex electronic systems, typical few-level approximations become unreliable.42,47 In this section, we present the Octopus implementation of an alternative real-space ab initio method for coupled light-matter systems. Dressed Reduced Density Matrix Functional Theory (dressed RDMFT)48 extends standard electronic RDMFT to coupled light-matter systems similarly to how QEDFT extends DFT. First tests on simple model systems suggest that the dressed RDMFT remains accurate from the weak to the ultra-strong coupling regime. A proper introduction of the theory, examples, and convergence studies can be found in Ref. 49.

This approach allows for the description of an interacting N-electron system coupled to one photonic mode within the dipole approximation. The respective Hamiltonian is given in Eq. (8) of Sec. III. Note that we set $j_{ext} = 0$ throughout this section and that the ground state $\Psi$ of Eq. (8) depends on $4N + 1$ coordinates, i.e., $\Psi = \Psi(r_1, \sigma_1, \dots, r_N, \sigma_N, p)$, where $\{\sigma_i\}$ denotes the spin degrees of freedom and $p$ is the elongation of the photon mode.

Within the dressed RDMFT, the original Hamiltonian (8) of N electrons in d dimensions and one mode is replaced by an extended auxiliary Hamiltonian of N dressed fermions in d + 1 dimensions with coordinates $z = (r, q) \in \mathbb{R}^{d+1}$. This auxiliary Hamiltonian reads

$$ \hat{H}' = \sum_{k=1}^{N} \left[ -\frac{1}{2} \Delta'_{k} + v'(z_{k}) \right] + \frac{1}{2} \sum_{k \neq l} w'(z_{k}, z_{l}) \quad (20) $$

and gives access to the same physics (see Refs. 49, Sec. 4). For a d = 1 matter subsystem, the operators introduced in Eq. (20) read

as follows: the dressed Laplacian $\Delta' = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial q^2}$, the dressed local potential

$$ v'(z) = v(x) + \left[ \frac{1}{2} \omega^2 q^2 - \frac{\omega}{\sqrt{N}} \lambda q x + \frac{1}{2} (\lambda x)^2 \right], \quad (21) $$

and the dressed interaction kernel

$$ w'(z, z') = w(x, x') + \left[ -\frac{\omega}{\sqrt{N}} \lambda q x' - \frac{\omega}{\sqrt{N}} \lambda q' x + \lambda^2 x x' \right]. \quad (22) $$

The ground state $\Psi'(z_1, \sigma_1, \dots, z_N, \sigma_N) = \Psi(x_1, \sigma_1, \dots, x_N, \sigma_N) \otimes \chi(p_2, \dots, p_N)$ of $\hat{H}'$ is a product of the original physical ground state $\Psi$ and the ground state of $\chi$, which, in turn, is the product of $N-1$ harmonic oscillator ground states. The auxiliary Hamiltonian (20) contains only one-body and two-body terms in terms of the dressed coordinates, which allows us, in principle, to use any standard electronic structure method to solve it (see also Refs. 49, Secs. 4 and 5). We use this construction to develop the dressed RDMFT and dressed Hartree-Fock (HF). For that, we define the dressed (spin-summed) first order reduced density matrix (1RDM)

$$ \begin{aligned} y'(z, z') = N \sum_{\sigma_1, \dots, \sigma_N} & \int d^{2(N-1)} z \Psi^{*}(z'_{1}, z_{2}, \dots, z_{N} \sigma_{N}) \\ & \times \Psi'(z \sigma_{1}, z_{2}, \dots, z_{N} \sigma_{N}). \end{aligned} \quad (23) $$

To apply RDMFT theory on the auxiliary system, we have to replace the total energy functional of electronic RDMFT, given in Ref. 19, with the newly introduced quantities of the dressed system, i.e., the auxiliary Hamiltonian of Eq. (20) (approximately) evaluated by the dressed 1RDM $y'$ of Eq. (23). By that, common approximations for the two-body energy expression in terms of the 1RDM can be directly transferred from electronic theory to the dressed system.50 The minimization is performed like in the electronic case and is based on the RDMFT implementation of Octopus, though the convergence of the dressed system requires a more complicated protocol that can be found in the supplementary material of Ref. 49. The current implementation in Octopus approximates the conditions under which the dressed 1RDM corresponds to a wave function by ensuring the fermionic ensemble N-representability conditions.51 However, the auxiliary wave function also exhibits another exchange symmetry with respect to the auxiliary coordinates, which is currently neglected. For the practical validity of this approach, the reader is referred to Sec. 5 and the supplementary material of Ref. 49.

<figure><img src="image_16.png" /><figcaption>FIG. 4. Differences of the dressed HF (dHF) and the dressed RDMFT (dRDMFT) from the exact ground state energies (in hartree, left) and from the exact photon number (right) as a function of the coupling $g/\omega$ for the stretched H$_2$ molecule in the dressed orbital description. The dressed RDMFT improves considerably upon dressed HF for the energy, especially due to the better description of the electronic correlation. The photon number, an example of a photonic observable, is captured similarly with both methods.</figcaption></figure>
As an example, we consider the one-dimensional (1D) H₂ molecule in a soft-Coulomb potential with a slightly stretched bond-length of b = 2.0 bohrs⁵² that is modeled by the local potential

$$v_{H_2}(x) = -\frac{1}{\sqrt{(x-b/2)^2 + 1}} - \frac{1}{\sqrt{(x+b/2)^2 + 1}} + \frac{1}{\sqrt{b^2 + 1}} \quad (24)$$

and the soft Coulomb interaction

$$w(x, x') = \frac{1}{\sqrt{|x - x'|^2 + 1}} \quad (25)$$

In Fig. 4, we show the total energy and the total photon number of the dressed RDMFT, the dressed HF, and the exact many-body calculations⁵³ for different coupling strengths. We see that for small couplings, both observables are captured well by the dressed RDMFT and dressed HF. With an increase in the coupling strength, both approximations fail to capture the strongly increasing photon number. For the total energy instead, the dressed RDMFT remains very close to the exact result, whereas the deviations to the dressed HF increase with increasing coupling strength. This shows the potential of the dressed RDMFT to describe correlated electron systems that are strongly coupled to a cavity mode. In the future, we plan to investigate better approximations to the polaritonic N-representability conditions that also account for the symmetry of the many-body wave function with respect to the exchange of photon coordinates.

## V. TOWARD DYNAMICS OF STRONGLY CORRELATED SYSTEMS: THE TDDFT+U FUNCTIONAL

It is well known that the standard local and semilocal functionals of DFT tend to over-delocalize the electrons, as usual approximations are based on the homogeneous electron gas. This leads to several failures of DFT for materials in which the localization of electrons plays a critical role in dictating the system's properties. This is, for instance, the case for transition metal oxides. The DFT+U method was originally proposed to compensate for some of the failures of the LDA for such materials.⁵⁴⁻⁵⁷ In essence, the DFT+U method aims at a better description of the local electron-electron interaction, which is achieved by adding the mean-field Hubbard model on a chosen localized subspace to the DFT total energy. The double counting of electron interaction in this localized subspace is then removed. The DFT+U total energy functional reads

$$E_{\text{DFT+U}}[n, \{n_{mm'}^{l,\sigma}\}] = E_{\text{DFT}}[n] + E_{(\text{ee})}[\{n_{mm'}^{l,\sigma}\}] - E_{\text{dc}}[\{n_{mm'}^{l,\sigma}\}], \quad (26)$$

where $E_{(\text{ee})}$ is the usual electron-electron interaction energy and $E_{\text{dc}}$ accounts for the double counting of the electron-electron interaction already present in $E_{\text{DFT}}$. The analytical expression of this double-counting term is not known in the general case, a general problem to all +U methods. Several approximated forms have been proposed along the years.⁵⁸,⁵⁹ In the Octopus code, we implemented the most commonly used double-counting terms: the fully localized limit (FLL)⁶⁰ and the around-mean field (AMF) double-counting terms.⁶¹ They read as⁵⁹

$$E_{\text{dc}}^{\text{FLL}}[\{n_{mm'}^{l,\sigma}\}] = \frac{U}{2}N(N-1) - \frac{J}{2}\sum_{\sigma} N^{\sigma}(N^{\sigma}-1) \quad (27)$$

and

$$E_{\text{dc}}^{\text{AMF}}[\{n_{mm'}^{l,\sigma}\}] = UN_{\uparrow}N_{\downarrow} - (U-J)\frac{2l}{2l+1}\sum_{\sigma} N_{\sigma}^{2}, \quad (28)$$

respectively, where $N_σ = \sum_m n_{mm}^σ$ and $N = N_↑ + N_↓$. The $E_{ee}$ and $E_{dc}$ energies depend on the density matrix of a localized orbital basis set {$φ_{I,m}$}, which are the localized orbitals attached to the atom I. In the following, we refer to the elements of the density matrix of the localized basis as occupation matrices and we denote them as {$n_{mm'}^{l,σ}$}. In the rotational-invariant form of DFT+U proposed by Dudarev et al.,⁶⁰ and for the FLL double-counting term, we obtain the $E_U$ energy to be added to the DFT total energy, which only depends on an effective Hubbard U parameter $U^{\text{eff}} = U - J$,

$$\begin{aligned} E_U[\{n_{mm'}^{l,\sigma}\}] &= E_{ee}[\{n_{mm'}^{l,\sigma}\}] - E_{dc}[\{n_{mm'}^{l,\sigma}\}] \\ &= \sum_{I,n,l} \frac{U_{I,n,l}^{\text{eff}}}{2} \sum_{m,\sigma} \left( n_{mm}^{l,n,l,\sigma} - \sum_{m'} n_{mm'}^{l,n,l,\sigma} n_{m'm}^{l,n,l,\sigma} \right), \quad (29) \end{aligned}$$

where I is an atom index, σ is the spin index, and n, l, and m refer to the principal, azimuthal, and angular quantum numbers, respectively. In the case of a periodic system, the occupation matrices $n_{mm'}^{l,n,l,σ}$ are given by

$$n_{mm'}^{l,n,l,\sigma} = \sum_{n} \sum_{k}^{\text{BZ}} w_{k} f_{nk}^{\sigma} \langle \psi_{n,k}^{\sigma} | \phi_{I,n,l,m} \rangle \langle \phi_{I,n,l,m'} | \psi_{n,k}^{\sigma} \rangle, \quad (30)$$

where $w_k$ is the k-point weight and $f_{nk}^σ$ is the occupation of the Bloch state $|\psi_{n,k}^σ\rangle$. Here, $|\phi_{I,n,l,m}\rangle$ are the localized orbitals that form the basis used to describe the electron localization. Details of the implementation can be found in Ref. 62. We recently extended our original implementation to be able to construct a localized subspace from localized states, such as Wannier states,⁶³ and to treat the intersite interaction.⁶⁴

In its usual formulation, the DFT+U method is an empirical method in which the effective U is a parameter of the calculation. However, it recently became possible to evaluate U and J fully ab initio and self-consistently using the ACBN0 functional.⁶⁵ We also implemented this method in Octopus and extended it to the time-dependent case⁶² to investigate strongly correlated materials out-of equilibrium. We showed that the absorption spectra of transition metal oxides, such as NiO or MnO, are well reproduced by our TDDFT+U simulations.⁶²

Figure 5 shows the calculated time profile of the effective Hubbard $U_{\text{eff}} = U - J$ for the 3d orbitals of Ni for light-driven NiO. The top panel exhibits the time profile of the driving vector potential. This shows that strongly driven correlated materials cannot be described by simply assuming that the effective electronic parameters (here, the effective Hubbard U) remain constant out of equilibrium. Moreover, the possibility to tune these effective electronic parameters offers opportunities for light-driven phase transitions, such as light-induced magnetic Weyl semimetals.⁶⁶

We finally note that TDDFT+U is not the only method that can be used to describe correlated materials such as transition metal oxides. An alternative approach to this problem would be to use hybrid functionals. However, such types of calculations in
As an example, we consider the one-dimensional (1D) H₂ molecule in a soft-Coulomb potential with a slightly stretched bond-length of b = 2.0 bohrs⁵² that is modeled by the local potential

$$v_{H_2}(x) = -\frac{1}{\sqrt{(x-b/2)^2+1}} - \frac{1}{\sqrt{(x+b/2)^2+1}} + \frac{1}{\sqrt{b^2+1}} \quad (24)$$

and the soft Coulomb interaction

$$w(x, x') = \frac{1}{\sqrt{|x - x'|^2 + 1}} \quad (25)$$

In Fig. 4, we show the total energy and the total photon number of the dressed RDMFT, the dressed HF, and the exact many-body calculations⁵³ for different coupling strengths. We see that for small couplings, both observables are captured well by the dressed RDMFT and dressed HF. With an increase in the coupling strength, both approximations fail to capture the strongly increasing photon number. For the total energy instead, the dressed RDMFT remains very close to the exact result, whereas the deviations to the dressed HF increase with increasing coupling strength. This shows the potential of the dressed RDMFT to describe correlated electron systems that are strongly coupled to a cavity mode. In the future, we plan to investigate better approximations to the polaritonic N-representability conditions that also account for the symmetry of the many-body wave function with respect to the exchange of photon coordinates.

## V. TOWARD DYNAMICS OF STRONGLY CORRELATED SYSTEMS: THE TDDFT+U FUNCTIONAL

It is well known that the standard local and semilocal functionals of DFT tend to over-delocalize the electrons, as usual approximations are based on the homogeneous electron gas. This leads to several failures of DFT for materials in which the localization of electrons plays a critical role in dictating the system's properties. This is, for instance, the case for transition metal oxides. The DFT+U method was originally proposed to compensate for some of the failures of the LDA for such materials.⁵⁴⁻⁵⁷ In essence, the DFT+U method aims at a better description of the local electron-electron interaction, which is achieved by adding the mean-field Hubbard model on a chosen localized subspace to the DFT total energy. The double counting of electron interaction in this localized subspace is then removed. The DFT+U total energy functional reads

$$E_{\text{DFT+U}}[n, \{n_{mm'}^{l,\sigma}\}] = E_{\text{DFT}}[n] + E_{(\text{ee})}[\{n_{mm'}^{l,\sigma}\}] - E_{\text{dc}}[\{n_{mm'}^{l,\sigma}\}], \quad (26)$$

where $E_{(\text{ee})}$ is the usual electron-electron interaction energy and $E_{\text{dc}}$ accounts for the double counting of the electron-electron interaction already present in $E_{\text{DFT}}$. The analytical expression of this double-counting term is not known in the general case, a general problem to all+U methods. Several approximated forms have been proposed along the years.⁵⁸,⁵⁹ In the Octopus code, we implemented the most commonly used double-counting terms: the fully localized limit (FLL)⁶⁰ and the around-mean field (AMF) double-counting terms.⁶¹ They read as⁵⁹

$$E_{\text{dc}}^{\text{FLL}}[\{n_{mm'}^{l,\sigma}\}] = \frac{U}{2}N(N-1) - \frac{J}{2}\sum_{\sigma} N^{\sigma}(N^{\sigma}-1) \quad (27)$$

and

$$E_{\text{dc}}^{\text{AMF}}[\{n_{mm'}^{l,\sigma}\}] = UN_{\uparrow}N_{\downarrow} - (U-J)\frac{2l}{2l+1}\sum_{\sigma} N_{\sigma}^2, \quad (28)$$

respectively, where $N_\sigma = \sum_m n_{mm}^\sigma$ and $N = N_\uparrow + N_\downarrow$. The $E_{ee}$ and $E_{dc}$ energies depend on the density matrix of a localized orbital basis set {$\phi_{l,m}$}, which are the localized orbitals attached to the atom I. In the following, we refer to the elements of the density matrix of the localized basis as occupation matrices and we denote them as {$n_{mm'}^{l,\sigma}$}. In the rotational-invariant form of DFT+U proposed by Dudarev et al.,⁶⁰ and for the FLL double-counting term, we obtain the $E_U$ energy to be added to the DFT total energy, which only depends on an effective Hubbard U parameter $U^{\text{eff}} = U - J$,

$$\begin{aligned} E_U[\{n_{mm'}^{l,\sigma}\}] &= E_{ee}[\{n_{mm'}^{l,\sigma}\}] - E_{dc}[\{n_{mm'}^{l,\sigma}\}] \\ &= \sum_{I,n,l} \frac{U_{I,n,l}^{\text{eff}}}{2} \sum_{m,\sigma} \left( n_{mm}^{l,n,l,\sigma} - \sum_{m'} n_{mm'}^{l,n,l,\sigma} n_{m'm}^{l,n,l,\sigma} \right), \end{aligned} \quad (29)$$

where I is an atom index, σ is the spin index, and n, l, and m refer to the principal, azimuthal, and angular quantum numbers, respectively. In the case of a periodic system, the occupation matrices $n_{mm'}^{l,n,l,\sigma}$ are given by

$$n_{mm'}^{l,n,l,\sigma} = \sum_{n} \sum_{k}^{\text{BZ}} w_{k} f_{nk}^{\sigma} \langle \psi_{n,k}^{\sigma} | \phi_{l,n,l,m} \rangle \langle \phi_{l,n,l,m'} | \psi_{n,k}^{\sigma} \rangle, \quad (30)$$

where $w_k$ is the k-point weight and $f_{nk}^\sigma$ is the occupation of the Bloch state $|\psi_{n,k}^\sigma\rangle$. Here, $|\phi_{l,n,l,m}\rangle$ are the localized orbitals that form the basis used to describe the electron localization. Details of the implementation can be found in Ref. 62. We recently extended our original implementation to be able to construct a localized subspace from localized states, such as Wannier states,⁶³ and to treat the intersite interaction.⁶⁴

In its usual formulation, the DFT+U method is an empirical method in which the effective U is a parameter of the calculation. However, it recently became possible to evaluate U and J fully ab initio and self-consistently using the ACBN0 functional.⁶⁵ We also implemented this method in Octopus and extended it to the time-dependent case⁶² to investigate strongly correlated materials out-of equilibrium. We showed that the absorption spectra of transition metal oxides, such as NiO or MnO, are well reproduced by our TDDFT+U simulations.⁶²

Figure 5 shows the calculated time profile of the effective Hubbard $U_{\text{eff}} = U - J$ for the 3d orbitals of Ni for light-driven NiO. The top panel exhibits the time profile of the driving vector potential. This shows that strongly driven correlated materials cannot be described by simply assuming that the effective electronic parameters (here, the effective Hubbard U) remain constant out of equilibrium. Moreover, the possibility to tune these effective electronic parameters offers opportunities for light-driven phase transitions, such as light-induced magnetic Weyl semimetals.⁶⁶

We finally note that TDDFT+U is not the only method that can be used to describe correlated materials such as transition metal oxides. An alternative approach to this problem would be to use hybrid functionals. However, such types of calculations in
<figure><img src="image_15.png" /><figcaption>FIG. 5. Self-consistent dynamics of Hubbard U for the Ni 3d orbitals (bottom panel) for pump intensities as indicated. The top panel represents the time-dependent vector potential, and the vertical dashed lines indicate the extrema of the vector potential, i.e., the minima of the driving electric field. Results are from Ref. 67.</figcaption></figure>

complex systems, potentially involving strong spin-orbit coupling interactions, are usually numerically very expensive in the context of real-time real-space TDDFT.

## VI. VAN DER WAALS INTERACTIONS

The van der Waals (vdW) interactions arise from correlations between electrons and are, in principle, described by the exact energy functional through the correlation energy functional $E_c[n(r)]$. However, the vdW interactions are inherently non-local and long-range and, by construction, cannot be described by usual local and semi-local functionals.⁶⁸ Therefore, much work has been devoted to finding consistent ways to enhance available functionals to correctly describe the necessary correlations yielding vdW forces.

Within a pure DFT approach, the inclusion of vdW interactions should be done through the correlation functional. To this end, a family of non-local density functionals has been proposed.⁶⁸ The so-called vdW density functionals (vdW-DF) are derived ab initio from the screened response of the homogeneous electron gas.⁶⁸ Another route to include vdW interactions in DFT calculations is by adding explicit corrections to the energy and forces based on atomic parameterizations. The most well-known method of this type is probably the vdW-D3 scheme of dispersion corrections from Grimme.⁶⁹ Although computationally more favorable than the vdW-DF method, the major disadvantage of this type of approach, particularly from the perspective of time-dependent simulations, is that it fails to correctly describe the effects that depend on the electron dynamics. A trivial example would be a lone electron traveling through space. Another scheme based on explicit corrections to the energy and forces was proposed by Tkatchenko and Scheffler (vdW-TS).⁷⁰ This scheme has the advantage of retaining much of the low computational cost of the vdW-D3 scheme while making the atomic parameterization dependent on the electronic density. Recently, these three schemes (vdW-DF, vdW-TS, and vdW-D3) were implemented in the Octopus code to deal with vdW interactions in isolated and periodic systems.

In the context of real-time TDDFT, it is important to note that these various levels of the description of vdW interactions have different implications. Indeed, in the vdW-D3 method, the energy and forces do not depend on the electronic density and will not change in time, unless atoms are moved, as already mentioned above. This should be a good approximation only when the time-dependent density remains close enough to the ground-state density. On the contrary, the expressions for the vdW-TS and the vdW-DF methods have a direct dependence on the electronic density and should thus yield a better description of vdW energies and forces out of equilibrium. However, we note that the use of the vdW-TS and vdW-DF methods within real-time TDDFT must be done using the adiabatic approximation and will, therefore, behave poorly when this approximation breaks down.

### A. vdW-DF

Octopus supports vdW-DF functionals⁶⁸ through the libvdwxc library.⁷¹ The vdW-DF functionals are expressed as a sum,

$$E_c^{\text{vdW-DF}}[n] = E_c^{\text{LDA}}[n] + E_c^{\text{nl}}[n], \quad (31)$$

of the LDA correlation energy⁷² and a fully non-local correlation term. The latter is the integral over a kernel function $\phi(q_0, q_0', r)$,

$$E_c^{\text{nl}} = \frac{1}{2} \iint n(\mathbf{r}) \phi(q_0(\mathbf{r}), q_0(\mathbf{r}'), |\mathbf{r} - \mathbf{r}'|) n(\mathbf{r}') \, d\mathbf{r} \, d\mathbf{r}', \quad (32)$$

where $q_0(r)$ depends on the local density and its gradient. An explicit evaluation of this six-dimensional integral is very expensive and scales as the volume squared, $O(N^2)$. Roman-Pérez and Soler proposed an efficient method, which approximates it as a sum of three-dimensional integrals.⁷³ The method works by expressing the integrand as a convolution using a limited set of helper functions and then applying the Fourier convolution theorem for a scaling of $O(N \log N)$ and a much lower prefactor. This has since become the standard method for evaluating vdW-DF functionals.

As implemented in Octopus, the density is redistributed from its normal uniform grid of arbitrary shape onto a uniform 3D grid forming a cube or parallelepiped, which is suitable for 3D Fourier transforms. The library libvdwxc then evaluates the non-local energy and contributions to the potential, relying on the FFTW library⁴ for efficient parallel Fourier transforms. After calculating the energy and potential, the potential is redistributed back to the original form.

Octopus supports the standard functionals vdW-DF1,⁷⁵,⁷⁶ vdW-DF2,⁷⁷ and vdW-DF-cx,⁷⁸ as well as other common forms that differ by combining the correlation with a different exchange functional. Some common supported variations are vdW-DF-optPBE, vdW-DF-optB88,⁷⁹ and vdW-DF-C09.⁸⁰

### B. vdW-TS

Since the vdW-TS approach depends on the density, the effect of the van der Waals interaction can be observed in properties other than the forces. In particular, we expect to observe an effect in the excited states of systems that interact through vdW forces. We use the hydrogen fluoride dimer as a simple model system for a proof-of-concept application of the modular implementation of the vdW-TS functional correction on TDDFT calculations. The dimer geometry is setup as shown in Fig. 6. The hydrogen fluoride monomers are
<figure><img src="image_17.png" /><figcaption>FIG. 5. Self-consistent dynamics of Hubbard U for the Ni 3d orbitals (bottom panel) for pump intensities as indicated. The top panel represents the time-dependent vector potential, and the vertical dashed lines indicate the extrema of the vector potential, i.e., the minima of the driving electric field. Results are from Ref. 67.</figcaption></figure>

complex systems, potentially involving strong spin-orbit coupling interactions, are usually numerically very expensive in the context of real-time real-space TDDFT.

## VI. VAN DER WAALS INTERACTIONS

The van der Waals (vdW) interactions arise from correlations between electrons and are, in principle, described by the exact energy functional through the correlation energy functional $E_c[n(r)]$. However, the vdW interactions are inherently non-local and long-range and, by construction, cannot be described by usual local and semi-local functionals.⁶⁸ Therefore, much work has been devoted to finding consistent ways to enhance available functionals to correctly describe the necessary correlations yielding vdW forces.

Within a pure DFT approach, the inclusion of vdW interactions should be done through the correlation functional. To this end, a family of non-local density functionals has been proposed.⁶⁸ The so-called vdW density functionals (vdW-DF) are derived ab initio from the screened response of the homogeneous electron gas.⁶⁸ Another route to include vdW interactions in DFT calculations is by adding explicit corrections to the energy and forces based on atomic parameterizations. The most well-known method of this type is probably the vdW-D3 scheme of dispersion corrections from Grimme.⁶⁹ Although computationally more favorable than the vdW-DF method, the major disadvantage of this type of approach, particularly from the perspective of time-dependent simulations, is that it fails to correctly describe the effects that depend on the electron dynamics. A trivial example would be a lone electron traveling through space. Another scheme based on explicit corrections to the energy and forces was proposed by Tkatchenko and Scheffler (vdW-TS).⁷⁰ This scheme has the advantage of retaining much of the low computational cost of the vdW-D3 scheme while making the atomic parameterization dependent on the electronic density. Recently, these three schemes (vdW-DF, vdW-TS, and vdW-D3) were implemented in the Octopus code to deal with vdW interactions in isolated and periodic systems.

In the context of real-time TDDFT, it is important to note that these various levels of the description of vdW interactions have different implications. Indeed, in the vdW-D3 method, the energy and forces do not depend on the electronic density and will not change in time, unless atoms are moved, as already mentioned above. This should be a good approximation only when the time-dependent density remains close enough to the ground-state density. On the contrary, the expressions for the vdW-TS and the vdW-DF methods have a direct dependence on the electronic density and should thus yield a better description of vdW energies and forces out of equilibrium. However, we note that the use of the vdW-TS and vdW-DF methods within real-time TDDFT must be done using the adiabatic approximation and will, therefore, behave poorly when this approximation breaks down.

### A. vdW-DF

Octopus supports vdW-DF functionals⁶⁸ through the libvdwxc library.⁷¹ The vdW-DF functionals are expressed as a sum,

$$E_c^{\text{vdW-DF}}[n] = E_c^{\text{LDA}}[n] + E_c^{\text{nl}}[n], \quad (31)$$

of the LDA correlation energy⁷² and a fully non-local correlation term. The latter is the integral over a kernel function $\phi(q_0, q_0', r)$,

$$E_c^{\text{nl}} = \frac{1}{2} \iint n(\mathbf{r}) \phi(q_0(\mathbf{r}), q_0(\mathbf{r}'), |\mathbf{r} - \mathbf{r}'|) n(\mathbf{r}') \, d\mathbf{r} \, d\mathbf{r}', \quad (32)$$

where $q_0(r)$ depends on the local density and its gradient. An explicit evaluation of this six-dimensional integral is very expensive and scales as the volume squared, $O(N^2)$. Roman-Pérez and Soler proposed an efficient method, which approximates it as a sum of three-dimensional integrals.⁷³ The method works by expressing the integrand as a convolution using a limited set of helper functions and then applying the Fourier convolution theorem for a scaling of $O(N \log N)$ and a much lower prefactor. This has since become the standard method for evaluating vdW-DF functionals.

As implemented in Octopus, the density is redistributed from its normal uniform grid of arbitrary shape onto a uniform 3D grid forming a cube or parallelepiped, which is suitable for 3D Fourier transforms. The library libvdwxc then evaluates the non-local energy and contributions to the potential, relying on the FFTW library⁴ for efficient parallel Fourier transforms. After calculating the energy and potential, the potential is redistributed back to the original form.

Octopus supports the standard functionals vdW-DF1,⁷⁵,⁷⁶ vdW-DF2,⁷⁷ and vdW-DF-cx,⁷⁸ as well as other common forms that differ by combining the correlation with a different exchange functional. Some common supported variations are vdW-DF-optPBE, vdW-DF-optB88,⁷⁹ and vdW-DF-C09.⁸⁰

### B. vdW-TS

Since the vdW-TS approach depends on the density, the effect of the van der Waals interaction can be observed in properties other than the forces. In particular, we expect to observe an effect in the excited states of systems that interact through vdW forces. We use the hydrogen fluoride dimer as a simple model system for a proof-of-concept application of the modular implementation of the vdW-TS functional correction on TDDFT calculations. The dimer geometry is setup as shown in Fig. 6. The hydrogen fluoride monomers are
<figure><img src="image_16.png" /><figcaption>FIG. 6. Absorption cross section spectrum of the hydrogen fluoride dimer with and without van der Waals effects. A zoomed-in view of the absorption feature at around 9.3 eV illustrates a small red shift when van der Waals effects are considered. Both hydrogen fluoride molecules are placed on the same plane at a separation of 2.8 Å.</figcaption></figure>

placed in an anti-parallel fashion, each one with its main symmetry axis oriented along the y-axis. The hydrogen fluoride bond in each monomer is 0.92 Å long, and the molecules are separated by 2.8 Å along the z-axis. At this distance, the van der Waals interaction between the monomers is the strongest according to the vdW-TS model. We choose this dimer model because it is a conveniently small system in which the effects of including the dispersion correction can be shown at an affordable computational cost.

To calculate absorption, we excite the system with an infinitesimal electric-field pulse and then propagate the time-dependent Kohn-Sham (TDKS) equations for 30.385 35 h/eV. The singlet dipole spectrum is evaluated from the time-dependent dipole moment. The strength of the perturbation is set to 0.01 Å⁻¹, and the polarization is in the z-axis. The time evolution is carried out using the enforced time-reversal symmetry propagator with (default) time steps of 0.033 52 h/eV.

The results (Fig. 6) show a small van der Waals-induced bathochromic-like (red) shift in the optical spectrum of the hydrogen fluoride dimer calculated with the LDA. This example opens the door for a new series of applications in supra-molecular chemistry, structural biology, and polymer science, which incorporate van der Waals effects on real-time electron dynamics.

### C. vdW-D3

Octopus also supports the DFT-D3 van der Waals correction.⁶⁹ This correction depends only on the atomic positions and does not depend on the atomic density. It is implemented in Octopus by linking to the DFT-D3 library provided by the authors. As we have validated Octopus results with the reference data provided with the library, the results are guaranteed to be consistent with other codes that implement this correction.

The only modification we did to the library was to move the very large set of coefficients from a source file to a stand alone text file that is parsed only when necessary. This speeds up compilation and reduces the size of the binaries. The modified library is distributed with Octopus, so users need not compile it separately.

## VII. POLARIZABLE CONTINUUM MODEL

The Polarizable Continuum Model (PCM)⁸¹ comprises a family of implicit-solvent approaches to tackle quantum-mechanical calculations of molecules in solution. The PCM assumes that (i) the solvent is a continuum and infinite dielectric medium characterized by a frequency-dependent dielectric function and (ii) a void cavity of appropriate shape and size encapsulates the solute molecule, separating it from the solvent by a sharp interface. The numerical implementation of the PCM relies on the Apparent Surface Charge (ASC) approach and the Boundary Element Method (BEM).⁸² In this framework, the solvent polarization response induced by the molecule's charge density is modeled by a reaction potential defined by a set of point charges $q = \{q_1, q_2, ..., q_T\}$ that spread over a tessellated cavity surface consisting of T finite surface elements or tesserae.⁸³ The reaction potential is the object through which the complex dielectric environment is accounted for.⁸⁴ The implementation of the PCM in Octopus rests on the Integral Equation Formalism (IEF).⁸⁵ The key IEF-PCM equation for computing the polarization charges is

$$q = QV, \qquad (33)$$

where $q$ and $V$ are column vectors of size T storing the induced polarization charges and the molecule's electrostatic potential at the tesserae representative points, respectively. $Q$ is the $T \times T$ PCM response matrix, which depends on the geometry of the cavity and the dielectric function of the solvent.⁸¹

A realistic description of the molecular cavity is key to capture accurate electronic and optical properties of molecules in solution. In principle, the cavity should (i) exclude the solvent, (ii) comprise most of the solute electronic density, and (iii) conform with the molecular shape. Octopus uses the GEPOL algorithm,⁸⁶ which builds up the van der Waals cavity from the union of interlocking spheres with element-specific radii centered at each atom position (by default, no spheres are built around hydrogen atoms). Within GEPOL, the BEM tessellation of the solute-solvent interface is done starting from a 60- or 240-face circumscribed polyhedron per sphere, selecting only exposed tesserae and properly reshaping those that are partially exposed. The BEM tessellation is a surface grid constructed independently from the real-space three-dimensional grid used in Octopus to represent both the Kohn-Sham electronic Hamiltonian and molecular orbitals. The mismatch between the two discrete representations might cause numerical problems arising from the Coulomb singularities whenever tesserae and grid points are close to each other. The implementation of the PCM in Octopus regularizes such singularities by using normalized spherical Gaussian functions to smooth the discretized polarization charges $q$.⁸⁴

Having briefly described the implementation, we now look at specific and relevant cases for chemical reactions, namely, solvation energies and ground-state stabilization. The most evident effect of the presence of a solvent is to change the total energy of a system compared to its value in vacuum. In the framework of DFT, the Kohn-Sham Hamiltonian of the solvated molecule contains the solvent reaction potential, which is a functional of the electronic and nuclear densities of the solute molecule.⁸⁴ The latter implies that the Kohn-Sham and IEF-PCM equations become coupled and the polarization charges $q_c = QV_{Hartree}[\rho^e]$, induced by the molecule's
<figure><img src="image_18.png" /><figcaption>FIG. 6. Absorption cross section spectrum of the hydrogen fluoride dimer with and without van der Waals effects. A zoomed-in view of the absorption feature at around 9.3 eV illustrates a small red shift when van der Waals effects are considered. Both hydrogen fluoride molecules are placed on the same plane at a separation of 2.8 Å.</figcaption></figure>

placed in an anti-parallel fashion, each one with its main symmetry axis oriented along the y-axis. The hydrogen fluoride bond in each monomer is 0.92 Å long, and the molecules are separated by 2.8 Å along the z-axis. At this distance, the van der Waals interaction between the monomers is the strongest according to the vdW-TS model. We choose this dimer model because it is a conveniently small system in which the effects of including the dispersion correction can be shown at an affordable computational cost.

To calculate absorption, we excite the system with an infinitesimal electric-field pulse and then propagate the time-dependent Kohn-Sham (TDKS) equations for 30.385 35 h/eV. The singlet dipole spectrum is evaluated from the time-dependent dipole moment. The strength of the perturbation is set to 0.01 Å⁻¹, and the polarization is in the z-axis. The time evolution is carried out using the enforced time-reversal symmetry propagator with (default) time steps of 0.033 52 h/eV.

The results (Fig. 6) show a small van der Waals-induced bathochromic-like (red) shift in the optical spectrum of the hydrogen fluoride dimer calculated with the LDA. This example opens the door for a new series of applications in supra-molecular chemistry, structural biology, and polymer science, which incorporate van der Waals effects on real-time electron dynamics.

### C. vdW-D3

Octopus also supports the DFT-D3 van der Waals correction.⁶⁹ This correction depends only on the atomic positions and does not depend on the atomic density. It is implemented in Octopus by linking to the DFT-D3 library provided by the authors. As we have validated Octopus results with the reference data provided with the library, the results are guaranteed to be consistent with other codes that implement this correction.

The only modification we did to the library was to move the very large set of coefficients from a source file to a stand alone text file that is parsed only when necessary. This speeds up compilation and reduces the size of the binaries. The modified library is distributed with Octopus, so users need not compile it separately.

## VII. POLARIZABLE CONTINUUM MODEL

The Polarizable Continuum Model (PCM)⁸¹ comprises a family of implicit-solvent approaches to tackle quantum-mechanical calculations of molecules in solution. The PCM assumes that (i) the solvent is a continuum and infinite dielectric medium characterized by a frequency-dependent dielectric function and (ii) a void cavity of appropriate shape and size encapsulates the solute molecule, separating it from the solvent by a sharp interface. The numerical implementation of the PCM relies on the Apparent Surface Charge (ASC) approach and the Boundary Element Method (BEM).⁸² In this framework, the solvent polarization response induced by the molecule's charge density is modeled by a reaction potential defined by a set of point charges $q = \{q_1, q_2, ..., q_T\}$ that spread over a tessellated cavity surface consisting of T finite surface elements or tesserae.⁸³ The reaction potential is the object through which the complex dielectric environment is accounted for.⁸⁴ The implementation of the PCM in Octopus rests on the Integral Equation Formalism (IEF).⁸⁵ The key IEF-PCM equation for computing the polarization charges is

$$q = QV, \qquad (33)$$

where $q$ and $V$ are column vectors of size $T$ storing the induced polarization charges and the molecule's electrostatic potential at the tesserae representative points, respectively. $Q$ is the $T \times T$ PCM response matrix, which depends on the geometry of the cavity and the dielectric function of the solvent.⁸¹

A realistic description of the molecular cavity is key to capture accurate electronic and optical properties of molecules in solution. In principle, the cavity should (i) exclude the solvent, (ii) comprise most of the solute electronic density, and (iii) conform with the molecular shape. Octopus uses the GEPOL algorithm,⁸⁶ which builds up the van der Waals cavity from the union of interlocking spheres with element-specific radii centered at each atom position (by default, no spheres are built around hydrogen atoms). Within GEPOL, the BEM tessellation of the solute-solvent interface is done starting from a 60- or 240-face circumscribed polyhedron per sphere, selecting only exposed tesserae and properly reshaping those that are partially exposed. The BEM tessellation is a surface grid constructed independently from the real-space three-dimensional grid used in Octopus to represent both the Kohn-Sham electronic Hamiltonian and molecular orbitals. The mismatch between the two discrete representations might cause numerical problems arising from the Coulomb singularities whenever tesserae and grid points are close to each other. The implementation of the PCM in Octopus regularizes such singularities by using normalized spherical Gaussian functions to smooth the discretized polarization charges $q$.⁸⁴

Having briefly described the implementation, we now look at specific and relevant cases for chemical reactions, namely, solvation energies and ground-state stabilization. The most evident effect of the presence of a solvent is to change the total energy of a system compared to its value in vacuum. In the framework of DFT, the Kohn-Sham Hamiltonian of the solvated molecule contains the solvent reaction potential, which is a functional of the electronic and nuclear densities of the solute molecule.⁸⁴ The latter implies that the Kohn-Sham and IEF-PCM equations become coupled and the polarization charges $q_c = QV_{Hartree}[\rho^e]$, induced by the molecule's
electronic density, have to be computed at each Kohn–Sham iteration until convergence is reached. The resulting electronic density can be used to compute the electrostatic contribution to the solvation free energy $\Delta G_{el} = G[\rho_{solv}^e] - E[\rho_{vac}^e]$, where G and E denote the free and total energy functionals of the molecule in solvent and in vacuum, respectively. The current PCM implementation in Octopus⁸⁴ has been thoroughly tested on a benchmark set of organic molecules and compared with analogous calculations performed with the quantum chemistry package GAMESS.⁸⁷ As an example of the solvent stabilization effect in the ground state, we studied the nitrobenzene molecule in water (static dielectric constant set to $\epsilon_s = 80$).

In Fig. 7, we show the convergence of the solvation free energy as a function of the Kohn–Sham iteration. The solvation free energy is stabilized already after 10 Kohn–Sham iterations out of the 19 required to optimize the molecular orbitals of the solvated molecule. We have also verified that the numerical error inherent to the discretization of the solute cavity surface is very small. For example, the total polarization charge $Q_{PCM} = \sum_{i=1}^{T} q_i^e$ at each Kohn–Sham iteration deviates by only 0.03% from the actual number of valence electrons in the nitrobenzene molecule (the relation between q and the solute charge is determined by Gauss’s theorem⁸¹). However, the magnitude of this error will depend, in general, on the size and the geometry of the molecule.

Now we move forward and describe the extension of the previous PCM method to the time domain using real-time PCM and non-equilibrium solvation. The ground state PCM is unable to capture the complex dynamical interactions between the solvent and the solute when the latter is in an excited state. Fortunately, the PCM admits a generalization to account for solvation dynamics.⁸⁹ The time-dependent PCM (TD-PCM) implementation in Octopus comes in three flavors of increasing complexity, all coupled with real-time TDDFT calculations of the solute molecules. The first

one, called equilibrium TD-PCM,⁸⁴ assumes that the solvent is fast enough to instantaneously equilibrate the solute charge density fluctuations and to polarize accordingly. This approach is physically sound for weakly polar solvents, having similar values for the static and dynamic dielectric constants. The second is a nonequilibrium approach called inertial TD-PCM and amounts to partitioning the solvent response in a fast (dynamic) and a slow (inertial) part.⁹⁰ Faster degrees of freedom respond instantaneously to the changes of the applied potential (either of molecular or external in origin), whereas slower degrees of freedom remain “frozen” and in equilibrium with the initial value of the field. This approximation works well when the solvent relaxation times are large enough with respect to the electronic excitations in the solute molecule (e.g., within the picosecond scale). The third TD-PCM approach, called equation of motion (EOM) TD-PCM, considers the full history-dependent evolution of the solvent polarization through a set of equations of motion for the polarization charges.⁹¹ Nonequilibrium polarization effects of this sort originate from the frequency-dependent dielectric response of the solvent, encoding the fact that it takes a different non-negligible time to adjust to fast or slow electrostatic perturbations. Solvation dynamics affect strongly and non-trivially the absorption spectrum of molecules, especially for fast and polar solvents, by inducing solvatochromic shifts of the peaks in the UV-visible absorption spectrum and modifying their relative amplitudes. The details about the implementation of all of these schemes and a detailed discussion of their effects can be found in Refs. 84 and 92.

Here, we show the differences among the TD-PCM flavors described above for the case of the nitrobenzene molecule in aqueous solution, excited with an electrical dipole perturbation in the x direction. We take the dynamic dielectric constant of water, entering in the inertial and EOM TD-PCMs, as $\epsilon_d = 1.786$. Figure 8 shows the photo-absorption spectrum of the gas-phase and solvated

<figure><img src="image_17.png" /><figcaption>FIG. 7. (Left panel) Nitrobenzene molecule in water ($\epsilon_s = 80$) surrounded by the ground state PCM charges distributed on the solvation cavity. The actual distribution of the charges in 3D is given in the supplementary material, SM1. (Right panel) Electrostatic contribution to the solvation free energy computed at each Kohn–Sham iteration. Kohn–Sham equations were solved in real-space as implemented in Octopus using the GGA-PBE approximation⁸⁸ to the exchange–correlation energy. The simulation box was built using spheres of radius 5 Å and an uniform spacing of 0.19 Å between grid points. The radii of the spheres used to build the cavity surface of the solute molecule are 2.4 Å for carbon, 1.8 Å for oxygen, and 1.9 Å for nitrogen.</figcaption></figure>
electronic density, have to be computed at each Kohn–Sham iteration until convergence is reached. The resulting electronic density can be used to compute the electrostatic contribution to the solvation free energy $ΔG_{el} = G[ρ_{solv}^e] - E[ρ_{vac}^e]$, where G and E denote the free and total energy functionals of the molecule in solvent and in vacuum, respectively. The current PCM implementation in Octopus⁸⁴ has been thoroughly tested on a benchmark set of organic molecules and compared with analogous calculations performed with the quantum chemistry package GAMESS.⁸⁷ As an example of the solvent stabilization effect in the ground state, we studied the nitrobenzene molecule in water (static dielectric constant set to $ε_s = 80$).

In Fig. 7, we show the convergence of the solvation free energy as a function of the Kohn–Sham iteration. The solvation free energy is stabilized already after 10 Kohn–Sham iterations out of the 19 required to optimize the molecular orbitals of the solvated molecule. We have also verified that the numerical error inherent to the discretization of the solute cavity surface is very small. For example, the total polarization charge $Q_{PCM} = \sum_{i=1}^{T} q_{i}^{e}$ at each Kohn–Sham iteration deviates by only 0.03% from the actual number of valence electrons in the nitrobenzene molecule (the relation between q and the solute charge is determined by Gauss’s theorem⁸¹). However, the magnitude of this error will depend, in general, on the size and the geometry of the molecule.

Now we move forward and describe the extension of the previous PCM method to the time domain using real-time PCM and non-equilibrium solvation. The ground state PCM is unable to capture the complex dynamical interactions between the solvent and the solute when the latter is in an excited state. Fortunately, the PCM admits a generalization to account for solvation dynamics.⁸⁹ The time-dependent PCM (TD-PCM) implementation in Octopus comes in three flavors of increasing complexity, all coupled with real-time TDDFT calculations of the solute molecules. The first

one, called equilibrium TD-PCM,⁸⁴ assumes that the solvent is fast enough to instantaneously equilibrate the solute charge density fluctuations and to polarize accordingly. This approach is physically sound for weakly polar solvents, having similar values for the static and dynamic dielectric constants. The second is a nonequilibrium approach called inertial TD-PCM and amounts to partitioning the solvent response in a fast (dynamic) and a slow (inertial) part.⁹⁰ Faster degrees of freedom respond instantaneously to the changes of the applied potential (either of molecular or external in origin), whereas slower degrees of freedom remain “frozen” and in equilibrium with the initial value of the field. This approximation works well when the solvent relaxation times are large enough with respect to the electronic excitations in the solute molecule (e.g., within the picosecond scale). The third TD-PCM approach, called equation of motion (EOM) TD-PCM, considers the full history-dependent evolution of the solvent polarization through a set of equations of motion for the polarization charges.⁹¹ Nonequilibrium polarization effects of this sort originate from the frequency-dependent dielectric response of the solvent, encoding the fact that it takes a different non-negligible time to adjust to fast or slow electrostatic perturbations. Solvation dynamics affect strongly and non-trivially the absorption spectrum of molecules, especially for fast and polar solvents, by inducing solvatochromic shifts of the peaks in the UV-visible absorption spectrum and modifying their relative amplitudes. The details about the implementation of all of these schemes and a detailed discussion of their effects can be found in Refs. 84 and 92.

Here, we show the differences among the TD-PCM flavors described above for the case of the nitrobenzene molecule in aqueous solution, excited with an electrical dipole perturbation in the x direction. We take the dynamic dielectric constant of water, entering in the inertial and EOM TD-PCMs, as $ε_d = 1.786$. Figure 8 shows the photo-absorption spectrum of the gas-phase and solvated

<figure><img src="image_19.png" /><figcaption>FIG. 7. (Left panel) Nitrobenzene molecule in water ($ε_s = 80$) surrounded by the ground state PCM charges distributed on the solvation cavity. The actual distribution of the charges in 3D is given in the supplementary material, SM1. (Right panel) Electrostatic contribution to the solvation free energy computed at each Kohn–Sham iteration. Kohn–Sham equations were solved in real-space as implemented in Octopus using the GGA-PBE approximation⁸⁸ to the exchange–correlation energy. The simulation box was built using spheres of radius 5 Å and an uniform spacing of 0.19 Å between grid points. The radii of the spheres used to build the cavity surface of the solute molecule are 2.4 Å for carbon, 1.8 Å for oxygen, and 1.9 Å for nitrogen.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**FIG. 8.** Absorption spectrum of nitrobenzene *in vacuo* and in water ($\epsilon_s$ = 80, $\epsilon_d$ = 1.786) for all the TD-PCM schemes. The inset shows the normalized cavity field factor (CFF) vs the excitation energy for the two EOM TD-PCMs with solvent relaxation times of 3.37 ps (slow) and 1 fs (fast) corresponding to CFFmax = 1.27 and 2.55, respectively. Real-time TDDFT simulations of 20 fs with a time step of 1.7 × 10⁻³ fs were performed to obtain the absorption spectra. An electric dipole perturbation within the linear regime and oriented along x was applied as an initial perturbation. The rest of the computational details are shared with the ground-state DFT calculations leading to Fig. 7.</figcaption></figure>

molecule. We can see that, with respect to the absorption *in vacuo*, all of the TD-PCM methods produce shifts of the features (in this case, toward lower excitation energies). The shifts are not rigid overall but depend on the excited state and on the specific solvation scheme (equilibrium, inertial, and EOM). Within Debye's model, the equilibrium and inertial TD-PCMs are limiting cases for the dynamics when the relaxation time is zero and infinite, respectively. The EOM TD-PCM, with a finite relaxation time, interpolates between these limits. Water has a large relaxation time of 3.37 ps, and therefore, the absorption spectrum is almost coincident for the EOM and the inertial TD-PCMs, as shown in Fig. 8. Whenever the solvent is as slow as water, there is not much gain in selecting the EOM over inertial TD-PCM and the latter is the method of choice in terms of computational performance. Instead, for faster solvents, the EOM results depart from those of the inertial algorithm, approaching those of the equilibrium TD-PCM. In Fig. 8, we considered an effective solvent having the same static and dynamic dielectric constants as water, but with a 1 fs relaxation time. The EOM TD-PCM for such a model solvent produces almost the same excitation energy shift equilibrium TD-PCM and peak intensities in between the inertial and the equilibrium TD-PCMs, as expected.

A real-time representation of the molecular dipole coupled with the PCM surface charges is found in the *supplementary material*, SM2, which also highlights how the different TD-PCM approaches affect the time-evolution of the dipole.

When studying the evolution of a solvated molecule under the effect of an external time-dependent electromagnetic field, a separate treatment is required to take into account that there are two solvent polarization contributions interacting with the solute molecule, namely, the reaction and the cavity fields.⁹³ The reaction field comes from the polarization induced by the solute molecule itself, while the cavity field arises as a polarization induced by the external electric field. Both reaction and cavity field effects can be

accounted for in static and real-time quantum-mechanical calculations of molecules within the PCM and TD-PCM simulations in Octopus.⁹²

All TD-PCM calculations shown here were performed using both reaction- and cavity-field effects, although the redshift of the peaks in each TD-PCM scheme is mainly a feature of the reaction-field, as ascertained in test calculations (not shown) and as expected.⁹⁴ Cavity field effects impact directly on the peak intensities by making the absorption more favorable depending on how large the effective local field acting on the molecule is allowed to be by solvent dielectric properties and the geometry of the cavity (normally, reassembling the molecular shape). Still, cavity field effects can have a non-trivial influence in the absorption spectrum shape when considering non-equilibrium solvation dynamics (EOM TD-PCM) by changing the relative peak intensities. This effect can be seen in the inset of Fig. 8, where the normalized cavity field factor (CFF)—the ratio between absorption cross section with and without cavity field effects—is plotted against the excitation energy for the EOM TD-PCM simulations with water and the aforementioned faster water-like solvent. The large difference in absorption peak intensities between the equilibrium and the rest of the TD-PCMs in Fig. 8 is also related to the modification of the probing electromagnetic field inside the dielectric. The photo-absorption cross section is, by definition, the ratio between the absorbed and the incoming power of light. In our case, both increase with the dielectric constant: the larger the dielectric constant of a medium, the larger the CFF, but also a larger dielectric constant implies a smaller phase velocity of light, therefore increasing the power of the traveling electromagnetic wave. The impact of the latter effect ($\propto$ refractive index = $\sqrt{(|\epsilon(\omega)| + \Re\epsilon(\omega))/2}$) is stronger than that of the former [$\propto$ CFF, e.g., for a spherical cavity⁹³ $\propto \left| \frac{3\epsilon(\omega)}{2\epsilon(\omega)+1} \right|^2$] for a large enough dielectric constant such as the one corresponding to the equilibrium TD-PCM for water ($\epsilon_s$ = 80).

In conclusion, Octopus is now capable of including an implicit dielectric continuum model as an environment for quantum mechanical calculations of excited states in the time-domain. The different versions of the TD-PCM scheme implemented allow us to select the most suitable to capture the relevant physics, accounting for a full range of different relaxation and response times.

VIII. MAGNONS FROM REAL-TIME TDDFT

In the last couple of years, the first studies investigating magnetization dynamics from first principles in real time have emerged.⁹⁵⁻¹⁰⁰ Here, we are interested in transverse magnetic excitations, specifically magnons, which are long wavelength collective excitations with a typical energy of tens to hundreds of millielectron volts. We recently developed an alternative to the linear-response TDDFT formulation⁹⁸,⁹⁹,¹⁰¹ to compute the spin susceptibilities of magnetic systems based on real-time TDDFT.¹⁰² Our approach follows the work of Bertsch *et al.* for optical excitations. In the original work of Bertsch *et al.*, the system is perturbed by a sudden change in the vector potential, which induces a charge- and current-density response. To investigate magnons, we employ a “transverse magnetic kick,” which induces magnetic fluctuations in the system. To be more precise, our perturbation corresponds to an infinitely short
<figure><img src="image_20.png" /><figcaption>**FIG. 8.** Absorption spectrum of nitrobenzene *in vacuo* and in water ($\epsilon_s$ = 80, $\epsilon_d$ = 1.786) for all the TD-PCM schemes. The inset shows the normalized cavity field factor (CFF) vs the excitation energy for the two EOM TD-PCMs with solvent relaxation times of 3.37 ps (slow) and 1 fs (fast) corresponding to CFFmax = 1.27 and 2.55, respectively. Real-time TDDFT simulations of 20 fs with a time step of 1.7 × 10⁻³ fs were performed to obtain the absorption spectra. An electric dipole perturbation within the linear regime and oriented along *x* was applied as an initial perturbation. The rest of the computational details are shared with the ground-state DFT calculations leading to Fig. 7.</figcaption></figure>

molecule. We can see that, with respect to the absorption *in vacuo*, all of the TD-PCM methods produce shifts of the features (in this case, toward lower excitation energies). The shifts are not rigid overall but depend on the excited state and on the specific solvation scheme (equilibrium, inertial, and EOM). Within Debye's model, the equilibrium and inertial TD-PCMs are limiting cases for the dynamics when the relaxation time is zero and infinite, respectively. The EOM TD-PCM, with a finite relaxation time, interpolates between these limits. Water has a large relaxation time of 3.37 ps, and therefore, the absorption spectrum is almost coincident for the EOM and the inertial TD-PCMs, as shown in Fig. 8. Whenever the solvent is as slow as water, there is not much gain in selecting the EOM over inertial TD-PCM and the latter is the method of choice in terms of computational performance. Instead, for faster solvents, the EOM results depart from those of the inertial algorithm, approaching those of the equilibrium TD-PCM. In Fig. 8, we considered an effective solvent having the same static and dynamic dielectric constants as water, but with a 1 fs relaxation time. The EOM TD-PCM for such a model solvent produces almost the same excitation energy shift equilibrium TD-PCM and peak intensities in between the inertial and the equilibrium TD-PCMs, as expected.

A real-time representation of the molecular dipole coupled with the PCM surface charges is found in the *supplementary material*, SM2, which also highlights how the different TD-PCM approaches affect the time-evolution of the dipole.

When studying the evolution of a solvated molecule under the effect of an external time-dependent electromagnetic field, a separate treatment is required to take into account that there are two solvent polarization contributions interacting with the solute molecule, namely, the reaction and the cavity fields.⁹³ The reaction field comes from the polarization induced by the solute molecule itself, while the cavity field arises as a polarization induced by the external electric field. Both reaction and cavity field effects can be

accounted for in static and real-time quantum-mechanical calculations of molecules within the PCM and TD-PCM simulations in Octopus.⁹²

All TD-PCM calculations shown here were performed using both reaction- and cavity-field effects, although the redshift of the peaks in each TD-PCM scheme is mainly a feature of the reaction-field, as ascertained in test calculations (not shown) and as expected.⁹⁴ Cavity field effects impact directly on the peak intensities by making the absorption more favorable depending on how large the effective local field acting on the molecule is allowed to be by solvent dielectric properties and the geometry of the cavity (normally, reassembling the molecular shape). Still, cavity field effects can have a non-trivial influence in the absorption spectrum shape when considering non-equilibrium solvation dynamics (EOM TD-PCM) by changing the relative peak intensities. This effect can be seen in the inset of Fig. 8, where the normalized cavity field factor (CFF)—the ratio between absorption cross section with and without cavity field effects—is plotted against the excitation energy for the EOM TD-PCM simulations with water and the aforementioned faster water-like solvent. The large difference in absorption peak intensities between the equilibrium and the rest of the TD-PCMs in Fig. 8 is also related to the modification of the probing electromagnetic field inside the dielectric. The photo-absorption cross section is, by definition, the ratio between the absorbed and the incoming power of light. In our case, both increase with the dielectric constant: the larger the dielectric constant of a medium, the larger the CFF, but also a larger dielectric constant implies a smaller phase velocity of light, therefore increasing the power of the traveling electromagnetic wave. The impact of the latter effect ($\propto$ refractive index = $\sqrt{(|\epsilon(\omega)| + \Re\epsilon(\omega))}/2$) is stronger than that of the former [$\propto$ CFF, e.g., for a spherical cavity⁹³ $\propto \left| \frac{3\epsilon(\omega)}{2\epsilon(\omega)+1} \right|^2$] for a large enough dielectric constant such as the one corresponding to the equilibrium TD-PCM for water ($\epsilon_s$ = 80).

In conclusion, Octopus is now capable of including an implicit dielectric continuum model as an environment for quantum mechanical calculations of excited states in the time-domain. The different versions of the TD-PCM scheme implemented allow us to select the most suitable to capture the relevant physics, accounting for a full range of different relaxation and response times.

VIII. MAGNONS FROM REAL-TIME TDDFT

In the last couple of years, the first studies investigating magnetization dynamics from first principles in real time have emerged.⁹⁵⁻¹⁰⁰ Here, we are interested in transverse magnetic excitations, specifically magnons, which are long wavelength collective excitations with a typical energy of tens to hundreds of millielectron volts. We recently developed an alternative to the linear-response TDDFT formulation⁹⁸,⁹⁹,¹⁰¹ to compute the spin susceptibilities of magnetic systems based on real-time TDDFT.¹⁰² Our approach follows the work of Bertsch *et al.* for optical excitations. In the original work of Bertsch *et al.*, the system is perturbed by a sudden change in the vector potential, which induces a charge- and current-density response. To investigate magnons, we employ a “transverse magnetic kick,” which induces magnetic fluctuations in the system. To be more precise, our perturbation corresponds to an infinitely short
application of a Zeeman term,

$$ \delta\hat{H}_q(t) = \frac{\mathcal{B}}{2}\delta(t) \int d^3\mathbf{r} [e^{-i\mathbf{q}\cdot\mathbf{r}}\hat{\sigma}^+(\mathbf{r}) + e^{i\mathbf{q}\cdot\mathbf{r}}\hat{\sigma}^-(\mathbf{r})], \quad (34) $$

where $\mathbf{q}$ is the momentum of the spin wave we are exciting, $\mathcal{B}$ is the strength of the perturbation, and $\sigma^\pm = \sigma_x \pm i\sigma_y$ are Pauli matrices. In this expression, the $z$ axis is taken to be along the direction of the magnetization of the system before excitation. In case the system has a preferred magnetization direction due to the presence of spin-orbit coupling, usually referred as an “easy axis,” we perform a kick in the transverse direction with respect to this “easy axis” of the system.

The subsequent time evolution of the spin magnetization $\mathbf{m}(\mathbf{r}, t)$, governed by the time-dependent Schrödinger equation, is then computed and analyzed in Fourier space. If we perturb our system from its ground state and we assume linear response, we have directly that

$$ m_+(\mathbf{q}; \omega) = \chi_+(\mathbf{q}; \omega) \frac{\mathcal{B}}{2}, \quad (35) $$

where $\chi_+(\mathbf{q}; \omega)$ is the spin susceptibility we want to extract. A similar expression is obtained for $\chi_-(\mathbf{q}; \omega)$.

In order to access finite momenta, one typically has to employ large supercells to perform the dynamics, which can be computationally very expensive when a few meV energy resolution is needed. Fortunately, there is a way to circumvent the construction of supercells by using the so-called generalized Bloch theorem (GBT). The GBT has been introduced by Sandradskii¹⁰³ for the calculation of ground-state spin waves and requires the implementation of specific boundary conditions. We implemented the GBT for the real-time calculation of magnons, taking advantage that Octopus is a real-space finite difference code, for which we can easily specify any type of boundary conditions. However, it is important to note that this applies only in the absence of spin-orbit coupling. The boundary condition, described below, depends on the momentum $\mathbf{q}$ of the perturbation and acts differently depending on whether a state was originally “up” or “down” with respect to the unperturbed magnetization. This is determined by the sign of $\langle\Phi_{nk}|S_z|\Phi_{nk}\rangle$ just before the perturbation for each spinor state $|\Phi_{nk}\rangle$. If we label these states $\alpha$ and $\beta$, the boundary condition reads

$$ \begin{aligned} \Phi_{\alpha,kn}(\mathbf{r}, t) &= e^{i\mathbf{k}\cdot\mathbf{r}} \left( \begin{array}{c} u_{\alpha,kn}^{\uparrow}(\mathbf{r}, t) \\ e^{i\mathbf{q}\cdot\mathbf{r}} u_{\alpha,kn}^{\downarrow}(\mathbf{r}, t) \end{array} \right), \\ \Phi_{\beta,kn}(\mathbf{r}, t) &= e^{i\mathbf{k}\cdot\mathbf{r}} \left( \begin{array}{c} e^{-i\mathbf{q}\cdot\mathbf{r}} u_{\beta,kn}^{\uparrow}(\mathbf{r}, t) \\ u_{\beta,kn}^{\downarrow}(\mathbf{r}, t) \end{array} \right). \end{aligned} \quad (36) $$

We checked that performing the simulation using the GBT or the supercell approach leads to the same results up to numerical precision. We tested our approach for cubic Ni, Fe, and Co, which have been widely studied using linear response TDDFT, and we found that our results are in very good agreement with previous work,⁹⁶ thus validating our implementation. Figure 9 shows the results obtained for bulk nickel using the adiabatic local-density approximation. We used here a real-space grid spacing of 0.27 bohr, norm-conserving pseudopotentials, and we employed a 16 × 16 × 16 k-point grid shifted 4 times in order to resolve momenta $\mathbf{q}$ that are multiples of $2\pi/(16a)$, where $a$ is the lattice parameter of Ni,

<figure><img src="image_19.png" /><figcaption>FIG. 9. Calculated magnon dispersion: spin susceptibility of bulk Ni along the $\overline{\Gamma}X$ direction, displayed in logarithmic scale. Numerical details are given in the main text. The experimental data are taken from Ref. 104.</figcaption></figure>

which we took to be 3.436 Å. In order to obtain the linear response, we took $\mathcal{B} = 0.02$ and we propagated for 435.4 fs using a time step of 1.81 as. To reduce the numerical burden, symmetries were employed.

Let us now comment on the interest of the proposed approach. So far, we only used this new method for investigating weak magnetic kicks, where we recover the results from linear response theory. However, our approach does not rely on the assumptions of small perturbations and can be used to investigate nonlinear phenomena induced by a strong magnetic field as well as out-of-equilibrium situations where the system is kicked from an excited state. This is the strength of a real-time method.¹⁰⁵⁻¹⁰⁸ Moreover, we can investigate directly the coupling with other degrees of freedom, such as phonons or photons, without any new theory or code development. Using the real-space method, we benefit from the favorable scaling of the time propagation, which is linear in the number of states, whereas sum-over-state approaches usually scale quadratically with the number of states. Finally, our approach offers the great advantage of not requiring the use of an exchange-correlation kernel, as only the exchange-correlation potential is needed to perform a time propagation. This is very interesting in order to test new functionals and theory levels for which deriving the expression of the exchange-correlation kernel can become complicated. These different aspects will be investigated in future works.

## IX. ORBITAL MAGNETO-OPTICAL RESPONSE OF SOLIDS AND MOLECULES FROM A STERNHEIMER APPROACH

In this section, we review the implemented routines for magneto-optical phenomena, which arise from the loss of symmetry between left and right circularly polarized light in the presence of a magnetic field.¹⁰⁹,¹¹⁰ While magneto-optical spectra can be straightforwardly computed for molecules,¹¹¹⁻¹¹⁵ the theory for solids has been developed only recently.¹¹⁶ The reason is that external electromagnetic fields break the translational symmetry of periodic systems, which is formally expressed through the
application of a Zeeman term,

$$ \delta\hat{H}_q(t) = \frac{\mathcal{B}}{2}\delta(t) \int d^3\mathbf{r} [e^{-i\mathbf{q}\cdot\mathbf{r}}\hat{\sigma}^+(\mathbf{r}) + e^{i\mathbf{q}\cdot\mathbf{r}}\hat{\sigma}^-(\mathbf{r})], \quad (34) $$

where $\mathbf{q}$ is the momentum of the spin wave we are exciting, $\mathcal{B}$ is the strength of the perturbation, and $\sigma^\pm = \sigma_x \pm i\sigma_y$ are Pauli matrices. In this expression, the $z$ axis is taken to be along the direction of the magnetization of the system before excitation. In case the system has a preferred magnetization direction due to the presence of spin-orbit coupling, usually referred as an “easy axis,” we perform a kick in the transverse direction with respect to this “easy axis” of the system.

The subsequent time evolution of the spin magnetization $\mathbf{m}(\mathbf{r}, t)$, governed by the time-dependent Schrödinger equation, is then computed and analyzed in Fourier space. If we perturb our system from its ground state and we assume linear response, we have directly that

$$ m_+(\mathbf{q}; \omega) = \chi_+(\mathbf{q}; \omega) \frac{\mathcal{B}}{2}, \quad (35) $$

where $\chi_+(\mathbf{q}; \omega)$ is the spin susceptibility we want to extract. A similar expression is obtained for $\chi_-(\mathbf{q}; \omega)$.

In order to access finite momenta, one typically has to employ large supercells to perform the dynamics, which can be computationally very expensive when a few meV energy resolution is needed. Fortunately, there is a way to circumvent the construction of supercells by using the so-called generalized Bloch theorem (GBT). The GBT has been introduced by Sandradskii¹⁰³ for the calculation of ground-state spin waves and requires the implementation of specific boundary conditions. We implemented the GBT for the real-time calculation of magnons, taking advantage that Octopus is a real-space finite difference code, for which we can easily specify any type of boundary conditions. However, it is important to note that this applies only in the absence of spin-orbit coupling. The boundary condition, described below, depends on the momentum $\mathbf{q}$ of the perturbation and acts differently depending on whether a state was originally “up” or “down” with respect to the unperturbed magnetization. This is determined by the sign of $\langle\Phi_{nk}|S_z|\Phi_{nk}\rangle$ just before the perturbation for each spinor state $|\Phi_{nk}\rangle$. If we label these states $\alpha$ and $\beta$, the boundary condition reads

$$ \Phi_{\alpha,kn}(\mathbf{r}, t) = e^{i\mathbf{k}\cdot\mathbf{r}} \begin{pmatrix} u_{\alpha,kn}^{\uparrow}(\mathbf{r}, t) \\ e^{i\mathbf{q}\cdot\mathbf{r}} u_{\alpha,kn}^{\downarrow}(\mathbf{r}, t) \end{pmatrix}, \quad (36) \\ \Phi_{\beta,kn}(\mathbf{r}, t) = e^{i\mathbf{k}\cdot\mathbf{r}} \begin{pmatrix} e^{-i\mathbf{q}\cdot\mathbf{r}} u_{\beta,kn}^{\uparrow}(\mathbf{r}, t) \\ u_{\beta,kn}^{\downarrow}(\mathbf{r}, t) \end{pmatrix}. $$

We checked that performing the simulation using the GBT or the supercell approach leads to the same results up to numerical precision. We tested our approach for cubic Ni, Fe, and Co, which have been widely studied using linear response TDDFT, and we found that our results are in very good agreement with previous work,⁹⁶ thus validating our implementation. Figure 9 shows the results obtained for bulk nickel using the adiabatic local-density approximation. We used here a real-space grid spacing of 0.27 bohr, norm-conserving pseudopotentials, and we employed a 16 × 16 × 16 k-point grid shifted 4 times in order to resolve momenta $\mathbf{q}$ that are multiples of $2\pi/(16a)$, where $a$ is the lattice parameter of Ni,

<figure><img src="image_21.png" /><figcaption>FIG. 9. Calculated magnon dispersion: spin susceptibility of bulk Ni along the $\overline{\Gamma}X$ direction, displayed in logarithmic scale. Numerical details are given in the main text. The experimental data are taken from Ref. 104.</figcaption></figure>

which we took to be 3.436 Å. In order to obtain the linear response, we took $\mathcal{B} = 0.02$ and we propagated for 435.4 fs using a time step of 1.81 as. To reduce the numerical burden, symmetries were employed.

Let us now comment on the interest of the proposed approach. So far, we only used this new method for investigating weak magnetic kicks, where we recover the results from linear response theory. However, our approach does not rely on the assumptions of small perturbations and can be used to investigate nonlinear phenomena induced by a strong magnetic field as well as out-of-equilibrium situations where the system is kicked from an excited state. This is the strength of a real-time method.¹⁰⁵⁻¹⁰⁸ Moreover, we can investigate directly the coupling with other degrees of freedom, such as phonons or photons, without any new theory or code development. Using the real-space method, we benefit from the favorable scaling of the time propagation, which is linear in the number of states, whereas sum-over-state approaches usually scale quadratically with the number of states. Finally, our approach offers the great advantage of not requiring the use of an exchange-correlation kernel, as only the exchange-correlation potential is needed to perform a time propagation. This is very interesting in order to test new functionals and theory levels for which deriving the expression of the exchange-correlation kernel can become complicated. These different aspects will be investigated in future works.

## IX. ORBITAL MAGNETO-OPTICAL RESPONSE OF SOLIDS AND MOLECULES FROM A STERNHEIMER APPROACH

In this section, we review the implemented routines for magneto-optical phenomena, which arise from the loss of symmetry between left and right circularly polarized light in the presence of a magnetic field.¹⁰⁹,¹¹⁰ While magneto-optical spectra can be straightforwardly computed for molecules,¹¹¹⁻¹¹⁵ the theory for solids has been developed only recently.¹¹⁶ The reason is that external electromagnetic fields break the translational symmetry of periodic systems, which is formally expressed through the
unboundness of the position operator.¹¹⁷⁻¹¹⁹ The orbital response to magnetic fields is especially complicated to describe, as such fields lead to non-perturbative changes in wavefunctions and introduce vector coupling to electron dynamics.¹²⁰⁻¹²² Here, we will focus on calculations of changes in the linear optical response in the presence of the magnetic field within the Sternheimer approach.¹⁹,¹²³⁻¹²⁵

To treat uniform magnetic fields in periodic systems, we use the approach based on perturbation theory for the one-particle density matrix.¹¹⁶,¹²⁰,¹²² Such an approach allows us to work under purely periodic boundary conditions and to take into account automatically gauge invariance. A periodic and gauge-invariant counterpart $\tilde{O}$ is distinguished for any operator $O = O_{r_1 r_2}$ defined for two points $\mathbf{r}_1$ and $\mathbf{r}_2$ in real space

$$O_{\mathbf{r}_1 \mathbf{r}_2} = \tilde{O}_{\mathbf{r}_1 \mathbf{r}_2} \exp\left(-\frac{i}{c} \int_{\mathbf{r}_2}^{\mathbf{r}_1} A(\mathbf{r}) d\mathbf{r}\right), \quad (37)$$

where $A$ is the vector potential associated with external electric ($\mathbf{E} = -c^{-1} \partial/\partial t \mathbf{A}$, where $c$ is the speed of light) and magnetic ($\mathbf{B} = \nabla \times \mathbf{A}$) fields and the integral is taken along the straight line between points $\mathbf{r}_2$ and $\mathbf{r}_1$.

The time-dependent Liouville equation for the gauge-invariant counterpart $\tilde{\rho}$ of the one-particle density matrix to the first order in $\mathbf{E}$ and $\mathbf{B}$ takes the form¹¹⁶

$$-i \partial_t \tilde{\rho} + [H_0, \tilde{\rho}] = -\frac{1}{2} \left\{ \mathbf{E} + \frac{1}{c} \mathbf{V} \times \mathbf{B}, [\mathbf{r}, \tilde{\rho}] \right\} - [\delta \tilde{H}, \tilde{\rho}]. \quad (38)$$

Here, the commutator and anticommutator of operators $O^{(1)}$ and $O^{(2)}$ are denoted by $[O^{(1)}, O^{(2)}]$ and $\{O^{(1)}, O^{(2)}\}$, respectively, the velocity operator $\mathbf{V} = -i[\mathbf{r}, \tilde{\mathbf{H}}]$ is computed taking into account all non-local contributions to the Hamiltonian (e.g., non-local pseudopotentials), and the Hamiltonian is represented as $\tilde{\mathbf{H}} = H_0 + \delta \tilde{\mathbf{H}}$, where the difference $\delta \tilde{\mathbf{H}}$ between the gauge-invariant counterpart $\tilde{\mathbf{H}}$ of the Hamiltonian and unperturbed Hamiltonian $H_0$ is related to the local-field effects. Unlike the singular position operator $\mathbf{r}$, the commutator $[\mathbf{r}, \tilde{\rho}]$ of the position operator with the periodic function $\tilde{\rho}$ is well defined in Eq. (38) and corresponds to the derivative with respect to the wave vector, $i \partial_k \tilde{\rho}_k$, in reciprocal space. Differentiating Eq. (38), one finds the derivatives of the density matrix $\tilde{\rho}^{(P)} = \partial \tilde{\rho} / \partial P$ with respect to external perturbations $P(\mathbf{E}, \mathbf{B}, \text{etc.})$.

For TDDFT calculations in Octopus, $\rho$ is the Kohn–Sham density matrix. The nth order derivative $\tilde{\rho}^{(P)}$ describing the joint response to the perturbations $P = P_1 P_2 \dots P_n$ is divided into four blocks within and between the occupied (V) and unoccupied subspaces (C): $\tilde{\rho}_{VV}^{(P)} = P_v \tilde{\rho}^{(P)} P_v$, $\tilde{\rho}_{CC}^{(P)} = P_c \tilde{\rho}^{(P)} P_c$, $\tilde{\rho}_{VC}^{(P)} = P_v \tilde{\rho}^{(P)} P_c$, and $\tilde{\rho}_{CV}^{(P)} = P_c \tilde{\rho}^{(P)} P_v$, where $P_v = \rho^{(0)}$ and $P_c = 1 - P_v$ are the projectors onto the occupied and unoccupied bands. In accordance with the density matrix perturbation theory,¹²⁶ to get the elements $\tilde{\rho}_{CV}^{(P)}$, Eq. (38) is projected onto unperturbed Kohn–Sham wavefunctions $|\psi_{v\mathbf{k}}^{(0)}\rangle$ of occupied bands $v$ to give an equation for the response function $|\eta_{v\mathbf{k}}^{(P)}\rangle = \tilde{\rho}_{CV}^{(P)}(\Omega)|\psi_{v\mathbf{k}}^{(0)}\rangle$,

$$L_{v\mathbf{k}}(\Omega)|\eta_{v\mathbf{k}}^{(P)}\rangle = P_c R^{(P)}[\tilde{\rho}^{(n-1)}, \dots, \rho^{(0)}, n^{(P)}]|\psi_{v\mathbf{k}}^{(0)}\rangle. \quad (39)$$

Here, the operator on the left-hand side is given by $L_{v\mathbf{k}}(\Omega) = \Omega + H_0 - \epsilon_{v\mathbf{k}}$, where $\Omega$ is the frequency considered and $\epsilon_{v\mathbf{k}}$ is the energy of the unperturbed state $|\psi_{v\mathbf{k}}^{(0)}\rangle$. The operator $R$ comes from the right-hand side of Eq. (38) and is determined by the derivatives of the density matrix of the previous orders. If the local-field effects are taken into account, the right-hand side $R$ also depends on the derivative of the electron density $n^{(P)}(\mathbf{r}_1) = \rho^{(P)}(\mathbf{r}_1, \mathbf{r}_2)\delta(\mathbf{r}_1 - \mathbf{r}_2)$. In this case, Eq. (39) needs to be solved self-consistently. The calculations, in practice, work with the periodic parts of the wavefunctions, $|\mu_{v\mathbf{k}}^{(0)}\rangle$. The commutator $[\mathbf{r}, \tilde{\rho}]$ corresponding to $i \partial_k \tilde{\rho}_k$ in the reciprocal space is computed within the $\mathbf{k} \cdot \mathbf{p}$ theory.¹⁹,¹²⁴,¹²⁵ Equation (39) is solved using the efficient Sternheimer approach,¹⁹,¹²³⁻¹²⁵ where the function $|\eta_{v\mathbf{k}}^{(P)}(\Omega)\rangle$ that fits into Eq. (39) is found iteratively at each frequency $\Omega$. To avoid divergences at resonances, a small but finite imaginary frequency $i\delta$ is added to the frequency $\Omega_0$ of the external perturbation so that $\Omega = \Omega_0 + i\delta$.

Once the solution of Eq. (39) is known, the elements $\tilde{\rho}_{CV}^{(P)}$ are obtained as

$$\tilde{\rho}_{CV}^{(P)}(\Omega) = \int_{\text{BZ}} \frac{d\mathbf{k}}{(2\pi)^3} \sum_v |\eta_{v\mathbf{k}}^{(P)}(\Omega)\rangle\langle\psi_{v\mathbf{k}}^{(0)}| \quad (40)$$

The elements of $\tilde{\rho}^{(P)}$ between the occupied and unoccupied subspaces are found using the relation $\tilde{\rho}_{VC}^{(P)}(\Omega) = (\tilde{\rho}_{CV}^{(P)}(-\Omega^*))^*$, and for that, Eq. (39) is also solved for the frequency $-\Omega^*$.

To find the elements within the occupied $\tilde{\rho}_{VV}^{(P)}$ and unoccupied $\tilde{\rho}_{CC}^{(P)}$ subspaces, the idempotency condition $\rho = \rho\rho$ is used. In terms of the periodic counterpart $\tilde{\rho}$ of the density matrix and to the first order in the magnetic field, it is written as¹²⁶,¹²⁷

$$\tilde{\rho} = \tilde{\rho}\tilde{\rho} + \frac{i}{2c}\mathbf{B} \cdot [\mathbf{r}, \tilde{\rho}] \times [\mathbf{r}, \tilde{\rho}]. \quad (41)$$

The contribution $\alpha_{v\mu,\gamma}$ to the polarizability in the presence of the magnetic field ($\alpha_{v\mu} = \alpha_{0v\mu} + \alpha_{v\mu,\gamma}B_\gamma$) is finally obtained from the current response as

$$\alpha_{v\mu,\gamma}(\Omega) = \frac{i}{\Omega} \mathrm{Tr} [V_v \tilde{\rho}^{(E_\mu B_\gamma)}(\Omega)], \quad (42)$$

where indices $v$, $\mu$, and $\gamma$ are used to denote components of the vectors $\mathbf{V}$, $\mathbf{E}$, and $\mathbf{B}$. The dielectric tensor in the presence of the magnetic field is computed as $\epsilon_{v\mu} = \delta_{v\mu} + 4\pi\alpha_{v\mu}/V$, where $V$ is the volume of the unit cell. Note that according to the “$2n+1$” theorem,¹²⁷,¹²⁸ there is no need to calculate explicitly the second-order derivative $\rho^{(E_\mu B_\gamma)}$. Instead, $\alpha_{v\mu,\gamma}$ is expressed through the first-order derivatives to the perturbations $P = E_\mu$, $B_\gamma$ and a supplementary perturbation corresponding to a vector potential $P = A_\nu$ at frequency $-\Omega$.¹¹⁶

To test the formalism for solids, we applied it to bulk silicon and the corresponding results are shown in Fig. 10. The full details of the calculations can be found in Ref. 116. It is seen in Fig. 10(b) that even without accounting for excitonic effects, the calculated spectra Re/Im $\epsilon_{xy}$ for the transverse component of the dielectric tensor are already qualitatively similar to the experimental curves¹¹⁰ at the direct absorption edge. To model excitonic effects, we used the model from Ref. 129. The spectra computed by
unboundness of the position operator.¹¹⁷⁻¹¹⁹ The orbital response to magnetic fields is especially complicated to describe, as such fields lead to non-perturbative changes in wavefunctions and introduce vector coupling to electron dynamics.¹²⁰⁻¹²² Here, we will focus on calculations of changes in the linear optical response in the presence of the magnetic field within the Sternheimer approach.¹⁹,¹²³⁻¹²⁵

To treat uniform magnetic fields in periodic systems, we use the approach based on perturbation theory for the one-particle density matrix.¹¹⁶,¹²⁰,¹²² Such an approach allows us to work under purely periodic boundary conditions and to take into account automatically gauge invariance. A periodic and gauge-invariant counterpart $\tilde{O}$ is distinguished for any operator $O = O_{r_1 r_2}$ defined for two points $\mathbf{r}_1$ and $\mathbf{r}_2$ in real space

$$O_{\mathbf{r}_1 \mathbf{r}_2} = \tilde{O}_{\mathbf{r}_1 \mathbf{r}_2} \exp\left(-\frac{i}{c} \int_{\mathbf{r}_2}^{\mathbf{r}_1} A(\mathbf{r}) d\mathbf{r}\right), \quad (37)$$

where $A$ is the vector potential associated with external electric ($\mathbf{E} = -c^{-1} \partial/\partial t \mathbf{A}$, where $c$ is the speed of light) and magnetic ($\mathbf{B} = \nabla \times \mathbf{A}$) fields and the integral is taken along the straight line between points $\mathbf{r}_2$ and $\mathbf{r}_1$.

The time-dependent Liouville equation for the gauge-invariant counterpart $\tilde{\rho}$ of the one-particle density matrix to the first order in $\mathbf{E}$ and $\mathbf{B}$ takes the form¹¹⁶

$$-i \partial_t \tilde{\rho} + [H_0, \tilde{\rho}] = -\frac{1}{2} \left\{ \mathbf{E} + \frac{1}{c} \mathbf{V} \times \mathbf{B}, [\mathbf{r}, \tilde{\rho}] \right\} - [\delta \tilde{H}, \tilde{\rho}]. \quad (38)$$

Here, the commutator and anticommutator of operators $O^{(1)}$ and $O^{(2)}$ are denoted by $[O^{(1)}, O^{(2)}]$ and $\{O^{(1)}, O^{(2)}\}$, respectively, the velocity operator $\mathbf{V} = -i[\mathbf{r}, \tilde{\mathbf{H}}]$ is computed taking into account all non-local contributions to the Hamiltonian (e.g., non-local pseudopotentials), and the Hamiltonian is represented as $\tilde{\mathbf{H}} = H_0 + \delta \tilde{\mathbf{H}}$, where the difference $\delta \tilde{\mathbf{H}}$ between the gauge-invariant counterpart $\tilde{\mathbf{H}}$ of the Hamiltonian and unperturbed Hamiltonian $H_0$ is related to the local-field effects. Unlike the singular position operator $\mathbf{r}$, the commutator $[\mathbf{r}, \tilde{\rho}]$ of the position operator with the periodic function $\tilde{\rho}$ is well defined in Eq. (38) and corresponds to the derivative with respect to the wave vector, $i \partial_k \tilde{\rho}_k$, in reciprocal space. Differentiating Eq. (38), one finds the derivatives of the density matrix $\tilde{\rho}^{(P)} = \partial \tilde{\rho} / \partial P$ with respect to external perturbations $P(\mathbf{E}, \mathbf{B}, \text{etc.})$.

For TDDFT calculations in Octopus, $\rho$ is the Kohn-Sham density matrix. The nth order derivative $\tilde{\rho}^{(P)}$ describing the joint response to the perturbations $P = P_1 P_2 \dots P_n$ is divided into four blocks within and between the occupied (V) and unoccupied subspaces (C): $\tilde{\rho}_{VV}^{(P)} = P_v \tilde{\rho}^{(P)} P_v$, $\tilde{\rho}_{CC}^{(P)} = P_c \tilde{\rho}^{(P)} P_c$, $\tilde{\rho}_{VC}^{(P)} = P_v \tilde{\rho}^{(P)} P_c$, and $\tilde{\rho}_{CV}^{(P)} = P_c \tilde{\rho}^{(P)} P_v$, where $P_v = \rho^{(0)}$ and $P_c = 1 - P_v$ are the projectors onto the occupied and unoccupied bands. In accordance with the density matrix perturbation theory,¹²⁶ to get the elements $\tilde{\rho}_{CV}^{(P)}$, Eq. (38) is projected onto unperturbed Kohn-Sham wavefunctions $|\psi_{v\mathbf{k}}^{(0)}\rangle$ of occupied bands $v$ to give an equation for the response function $|\eta_{v\mathbf{k}}^{(P)}\rangle = \tilde{\rho}_{CV}^{(P)}(\Omega)|\psi_{v\mathbf{k}}^{(0)}\rangle$,

$$L_{v\mathbf{k}}(\Omega)|\eta_{v\mathbf{k}}^{(P)}\rangle = P_c R^{(P)}[\tilde{\rho}^{(n-1)}, \dots, \rho^{(0)}, n^{(P)}]|\psi_{v\mathbf{k}}^{(0)}\rangle. \quad (39)$$

Here, the operator on the left-hand side is given by $L_{v\mathbf{k}}(\Omega) = \Omega + H_0 - \epsilon_{v\mathbf{k}}$, where $\Omega$ is the frequency considered and $\epsilon_{v\mathbf{k}}$ is the energy of the unperturbed state $|\psi_{v\mathbf{k}}^{(0)}\rangle$. The operator $R$ comes from the right-hand side of Eq. (38) and is determined by the derivatives of the density matrix of the previous orders. If the local-field effects are taken into account, the right-hand side $R$ also depends on the derivative of the electron density $n^{(P)}(\mathbf{r}_1) = \rho^{(P)}(\mathbf{r}_1, \mathbf{r}_2)\delta(\mathbf{r}_1 - \mathbf{r}_2)$. In this case, Eq. (39) needs to be solved self-consistently. The calculations, in practice, work with the periodic parts of the wavefunctions, $|\mu_{v\mathbf{k}}^{(0)}\rangle$. The commutator $[\mathbf{r}, \tilde{\rho}]$ corresponding to $i \partial_k \tilde{\rho}_k$ in the reciprocal space is computed within the $\mathbf{k} \cdot \mathbf{p}$ theory.¹⁹,¹²⁴,¹²⁵ Equation (39) is solved using the efficient Sternheimer approach,¹⁹,¹²³⁻¹²⁵ where the function $|\eta_{v\mathbf{k}}^{(P)}(\Omega)\rangle$ that fits into Eq. (39) is found iteratively at each frequency $\Omega$. To avoid divergences at resonances, a small but finite imaginary frequency $i\delta$ is added to the frequency $\Omega_0$ of the external perturbation so that $\Omega = \Omega_0 + i\delta$.

Once the solution of Eq. (39) is known, the elements $\tilde{\rho}_{CV}^{(P)}$ are obtained as

$$\tilde{\rho}_{CV}^{(P)}(\Omega) = \int_{\text{BZ}} \frac{d\mathbf{k}}{(2\pi)^3} \sum_v |\eta_{v\mathbf{k}}^{(P)}(\Omega)\rangle\langle\psi_{v\mathbf{k}}^{(0)}| \quad (40)$$

The elements of $\tilde{\rho}^{(P)}$ between the occupied and unoccupied subspaces are found using the relation $\tilde{\rho}_{VC}^{(P)}(\Omega) = (\tilde{\rho}_{CV}^{(P)}(-\Omega^*))^*$, and for that, Eq. (39) is also solved for the frequency $-\Omega^*$.

To find the elements within the occupied $\tilde{\rho}_{VV}^{(P)}$ and unoccupied $\tilde{\rho}_{CC}^{(P)}$ subspaces, the idempotency condition $\rho = \rho\rho$ is used. In terms of the periodic counterpart $\tilde{\rho}$ of the density matrix and to the first order in the magnetic field, it is written as¹²⁶,¹²⁷

$$\tilde{\rho} = \tilde{\rho}\tilde{\rho} + \frac{i}{2c}\mathbf{B} \cdot [\mathbf{r}, \tilde{\rho}] \times [\mathbf{r}, \tilde{\rho}]. \quad (41)$$

The contribution $\alpha_{v\mu,\gamma}$ to the polarizability in the presence of the magnetic field ($\alpha_{v\mu} = \alpha_{0v\mu} + \alpha_{v\mu,\gamma}B_\gamma$) is finally obtained from the current response as

$$\alpha_{v\mu,\gamma}(\Omega) = \frac{i}{\Omega} \text{Tr} [V_v \tilde{\rho}^{(E_\mu B_\gamma)}(\Omega)], \quad (42)$$

where indices $v$, $\mu$, and $\gamma$ are used to denote components of the vectors $\mathbf{V}$, $\mathbf{E}$, and $\mathbf{B}$. The dielectric tensor in the presence of the magnetic field is computed as $\epsilon_{v\mu} = \delta_{v\mu} + 4\pi\alpha_{v\mu}/V$, where $V$ is the volume of the unit cell. Note that according to the “$2n+1$” theorem,¹²⁷,¹²⁸ there is no need to calculate explicitly the second-order derivative $\rho^{(E_\mu B_\gamma)}$. Instead, $\alpha_{v\mu,\gamma}$ is expressed through the first-order derivatives to the perturbations $P = E_\mu, B_\gamma$ and a supplementary perturbation corresponding to a vector potential $P = A_\nu$ at frequency $-\Omega$.¹¹⁶

To test the formalism for solids, we applied it to bulk silicon and the corresponding results are shown in Fig. 10. The full details of the calculations can be found in Ref. 116. It is seen in Fig. 10(b) that even without accounting for excitonic effects, the calculated spectra Re/Im $\epsilon_{xy}$ for the transverse component of the dielectric tensor are already qualitatively similar to the experimental curves¹¹⁰ at the direct absorption edge. To model excitonic effects, we used the model from Ref. 129. The spectra computed by
<figure><img src="image_20.png" /><figcaption>FIG. 10. Calculated components (a) $\epsilon_{xx}$ and (b) $\epsilon_{xy}$ of the dielectric tensor of silicon for different frequencies of light $\Omega_0$ (in eV) and the magnetic field of 1 T along the z axis (with and without account of excitonic effects). The linewidth $\delta = 0.1$ eV is used. The blue shift by 0.7 eV is applied to take into account the GW correction to the bandgap.<sup>130,131</sup> The experimental data for $\epsilon_{xx}$<sup>132</sup> and $\epsilon_{xy}$<sup>110</sup> (scaled by 1/2) are indicated by symbols.</figcaption></figure>

accounting for excitonic effects show a better agreement with the experimental results [Figs. 10(a) and 10(b)]. Although the magnitudes of the peaks in the magneto-optical spectra are about a factor of two smaller than in the magneto-optical measurements,<sup>110</sup> they can be corrected by reducing the linewidth $\delta$ assumed in the calculations.

In the limit of a large supercell, this formalism becomes equivalent to the simpler standard formulation for finite systems, which we have also implemented for reference. In this case, the Liouville equation for the density matrix is written in the Coulomb gauge as

$$ \Omega\rho + [H_0, \rho] = [\mathbf{d} \cdot \mathbf{E} + \mathbf{m} \cdot \mathbf{B} - \delta H, \rho], \quad (43) $$

where $\mathbf{d} = -\mathbf{r}$ is the electric dipole moment, $\mathbf{m} = -\mathbf{r} \times \mathbf{V}/2c$ is the orbital magnetic dipole moment, and $\delta H$ describes local-field effects. The change in the polarizability $\alpha_{v\mu,\gamma}B_\gamma$ in the presence of the magnetic field is calculated from the dipole response

$$ \alpha_{v\mu,\gamma}(\Omega) = \mathrm{Tr}[d_v\rho^{(E_\mu B_\gamma)}(\Omega)]. \quad (44) $$

As in the periodic case, due to the “$2n + 1$” theorem,<sup>127,128</sup> the explicit calculation of the second-order derivative $\rho^{(E_\mu B_\gamma)}$ is avoided. Instead, a supplementary electric field at frequency $-\Omega$, $E_v^\prime$, is introduced.<sup>116</sup> In the case of real wavefunctions, $|\eta_v^{(E_\mu)}(-\Omega)\rangle = (|\eta_v^{(E_\mu)}(-\Omega^*)\rangle)^*$, and there is no need to solve additionally the Liouville equation for this supplementary electric field.

The efficiency of the present scheme to compute magneto-optics is comparable to standard linear-response calculations of simple optical polarizability in the absence of the magnetic field.<sup>116</sup> When local-field effects are included self-consistently, the calculations of magneto-optical spectra for solids take only twice as long as those of polarizability. For finite systems, the computational effort is the same as for the simple optical polarizability.

## X. TIME-DEPENDENT ANGULAR RESOLVED PHOTOELECTRON SPECTROSCOPY

The real-space representation of the dynamics of the electronic structure allows for a seamless and straightforward description of dynamical processes outside the material, i.e., processes where electrons are excited into vacuum. Beyond simply describing the ionization process, Octopus has routines implemented that compute photoelectron spectroscopy in different flavors. Photoelectron spectroscopy is particularly ubiquitous for the characterization of the electronic structure in solids because it provides a direct observable of the energy and momentum distribution of electronic states, known as the bandstructure. In contrast to other electronic structure methods that compute quantities linked to the photoelectron spectrum, most commonly the single-quasiparticle spectral function obtained through the GW approximation to the many-body self-energy, the approach described here does not consider unit cells of bulk materials but instead computes directly the energy and momentum resolved ionization probability.<sup>107,133</sup> Besides the study of the electronic structure of solids, photoelectron spectra of atoms and molecules are of equal interest and the implementation in Octopus is capable to describe accurately all such systems on an equal footing.

The most detailed quantity available in the experiments is the momentum-resolved photoelectron probability $P(\mathbf{p})$, i.e., the probability to detect an electron with a given momentum $\mathbf{p}$. In some cases, the experimental setup offers the possibility to measure the spin polarization along a given axis. The formalism we are going to outline can easily accommodate a non-collinear spin structure and, therefore, calculate spin-resolved quantities, but for the sake of simplicity, in the following, we are going to restrict ourselves to closed shell systems with collinear spins, i.e., where all the orbitals are doubly occupied with electrons having opposite spins. The reader interested in the most general case can find a detailed description in Ref. 134. From $P(\mathbf{p})$, one can obtain other derived quantities by simple manipulation. For instance, the energy-resolved spectrum $P(E)$, used to identify the occupied energy levels, can be obtained with a change of variable using the free electron energy dispersion relation $E = \rho^2/2$. The angle-resolved photoelectron spectrum (ARPES), normally employed to measure the quasiparticle bandstructure, is simply obtained by taking the energy resolved spectrum as a function of the electron momentum parallel to the surface $P(\mathbf{p}_{\parallel}, E)$.

The t-SURFF method was first proposed by Tao and Scrinzi<sup>135</sup> for one-electron systems and later extended to many electrons with TDDFT for periodic<sup>134</sup> and non-periodic systems.<sup>136</sup> It is based on the assumption that the Kohn-Sham Hamiltonian describing the full experimental process, i.e., including the ionization and detection, can be decomposed into the sum of two Hamiltonians acting into complementary spatial regions, *inner* and *outer*, and that can be approximated in different ways. In the *inner* region surrounding the system, the electron dynamics is governed by the interacting Kohn-Sham Hamiltonian $\hat{H}_{\mathrm{KS}}(\mathbf{r}, t)$, while in the *outer* region, electrons are free from the Coulomb tails of the parent system and behave as independent particles driven by an external field and, therefore, are described by the Volkov Hamiltonian $\hat{H}_{\mathrm{V}}(\mathbf{r}, t) = 1/2(-i\nabla - \mathbf{A}(t)/c)^2$. We express the field with a time dependent vector potential $\mathbf{A}(t)$ in the dipole approximation, i.e., by discarding the spatial dependence of the field $\mathbf{A}(\mathbf{r}, t) \approx \mathbf{A}(t)$. The