Ab initio modeling of excitons: from perfect crystals to biomaterials

ABSTRACT Excitons, or coupled electron-hole excitations, are important both for fundamental optical properties of materials as well as and for the functionality of materials in opto-electronic devices. Depending on the material they are created in, excitons can come in many forms, from Wannier–Mott excitons in inorganic semiconductors, to molecular Frenkel or bi-molecular charge-transfer excitons in disordered organic or biological heterostructures. This multitude of materials and exciton types poses tremendous challenges for ab initio modeling. Following a brief overview of typical ab initio techniques, we summarize our recent work based on Many-Body Green’s Functions Theory in the GW approximation and Bethe–Salpeter Equation (BSE) as a method applicable to a wide range of material classes from perfect crystals to disordered materials. In particular, we emphasize the current challenges of embedding this GW-BSE method into multi-method, mixed quantum-classical (QM/MM) models for organic materials and illustrate them with examples from organic photovoltaics and fluorescence spectroscopy. Our perspectives on future studies are also presented. Graphical Abstract


Introduction
When an incident photon is absorbed by a material, an electron is excited from the valence to the conduction band and leaves behind a positively charged vacancy, a hole. The attractive Coulomb interaction between the excited electron and the hole binds them together to form a bound neutral compound system of the two charge carriers, reminiscent of a hydrogen atom. This coupled excited electron-hole pair is called an exciton. An exciton can also be formed by different processes, e.g. from the injection of free electrons and holes into the material.
Excitons play a pivotal role in many different phenomena in materials. They determine their optical properties, such as absorption but also emission via fluorescence or phosphorescence [1]. Strongly or weakly coupled electron-hole pairs are not only static objects: they can diffuse through a material, thereby transporting energy with zero net charge [2][3][4]. Excitons can also play a role in chemical reactions [5], be intermediates to the generation of free charges in devices [6], and drive biological processes [7], such as photosynthesis. The initial step in these processes of providing chemical energy for plants, algae, and bacteria involves the capture of energy from sunlight. Specialized pigmentprotein complexes, called light-harvesting antenna complexes, consist of lightabsorbing chromophores, typically attached to a protein structure that holds them in place. Before the molecule can relax, the electronic excitation must be 'harvested'. That is, the excitation is transferred through space among the chromophores until it eventually reaches a reaction center where it initiates charge separation. Technological applications often exploit the static and dynamic properties of excitons in a similar fashion. Considerable effort is directed at designing either pure or composite materials with target properties to enhance device characteristics, such as optimizing the color of emitted light and increasing stability in organic light-emitting diodes (OLEDs) [8][9][10], or boosting the power conversion efficiency in solar cells, to name only a few.
A rational design of materials fulfilling these goals relies on the understanding of how the excitonic properties are related to material attributes, e.g. composition, chemical bonding, structural order, etc. On the nanoscale such an understanding is often sought with the help of ab initio modeling approaches, which can predict -if employed appropriately -the energetics of excitons, the binding characteristics of the electron-hole pair, and reveal how they emerge from the underlying material structure. There is a plethora of possible structures ranging form perfectly ordered, crystalline atomic lattice structures to complex disordered morphologies in biomaterials.
From the perspective of ab initio modeling, this multitude of material types constitutes a tremendous challenge, and no one-size-fits-all solutions can be expected. In this review, we will in the following briefly first revisit the relation between material properties and exciton characteristics using a rough distinction between inorganic and organic materials. We will give a few examples of the role exciton processes play in organic light-emitting diodes and photovoltaic cells. We will also provide a concise overview of typical ab initio techniques employed in the literature for the study of excitons and their general applicability from a methodological and computational perspective to the various material types. Special attention is given to Many-Body Green's Functions Theory in the GW approximation and Bethe-Salpeter Equation (BSE) as a method applicable to a wide range of material classes from perfect crystals to biomaterials. We showcase its use and the related challenges with examples from the literature for inorganic and organic crystals and summarize our recent work on complex, disordered molecular materials. In particular, we emphasize the current challenges of embedding this GW-BSE method into multi-method, mixed quantumclassical (QM/MM) models and illustrate them with examples from organic photovoltaics and fluorescence spectroscopy.

Relation between materials properties and exciton character
Details of the static and dynamic properties of excitons depend on the material in which they form. Materials are in broad terms classified as either inorganic or organic as they differ in their chemical bonding and electronic characteristics. We illustrate the different materials and types of excitons schematically in Figure 1.
Inorganic materials are usually well-ordered systems with atoms arranged in regular crystal structures, with either covalent or ionic bonds. As a result, they are structurally rigid and allow very little room for modifications. Traditional inorganic semiconductors such as silicon, germanium, and GaAs have low band gaps ([0.67] eV, [1.11] eV, and [1.43] eV respectively [11]). The dielectric constant in these materials is large (� r > 10) so that Coulomb effects between electrons and holes are small due to dielectric screening, with exciton binding energies of only several meV. The excitonic wave function typically extends over several lattice spacings, as indicated in Figure 1(a). This is known as a Wannier-Mott exciton [12]. For larger bandgap and/or lower dielectric constant, the screening effects become smaller, the exciton binding energy increases to the order of [0.1-1] eV as the excitonic wave function becomes more localized in space, e.g. in alkali halide crystals. Typically, such small-radius excitons are referred to as Frenkel excitons [13]. It should be noted, however, that in many situations, the distinction between Wannier-Mott and Frenkel excitons, or Frenkel excitons of different localization lengths, is not sharp, in particular in disordered materials without any clear lattice spacing.
Organic materials are formed from molecular building blocks and can host Frenkel excitons of vastly varying localization character [14]. While the intramolecular bonding is determined by covalent and ionic interactions, weak cohesive electrostatic and van-der-Waals interactions are responsible for the inter-molecular structure formation. As a result, there is in general a plethora of different material structures, ranging from well-ordered organic crystals (mostly at low temperature) to statically and/or thermally disordered structures with varying dimensionality. Often, there are even multiple phases present in the same organic material or material composite as technologically relevant materials are often created from solution processing. Organic materials combine chemical and mechanical benefits of organic compounds, such as tailoring of electronic properties by a modification by chemical synthesis, their light-weight, and flexibility with peculiar properties of semiconductor materials, e.g. the absorption/emission of light in the visible spectral range and conductivity that is sufficient for the operation of devices such as light-emitting diodes (LEDs), solar cells, and field-effect-transistors. As a consequence of the intrinsically weak inter-molecular interactions and resulting disorder, electronic states are typically Highly ordered molecular crystal structure where each oval represents a molecular unit. The exciton is delocalized over several molecules due to strong intermolecular excitonic coupling. (c) Disordered (amorphous) molecular material, with Frenkel excitons strongly localized on single molecules. (d) Interface of a donor-acceptor heterostructure of two disordered molecular materials with bimolecular charge-transfer excitons. localized on one or several molecular building blocks, leading to low values of the dielectric constant usually in the region of � r ¼ 3-4. Absorption and emission take place mostly in the range of 2-3eV, hindering any significant chargecarrier concentration by thermal excitation at room temperature. As mentioned above, the low dielectric constant implies large exciton binding energies and small exciton radii. How small depends on the structural details of the material. Figure 1(b) schematically depicts a highly ordered molecular crystal. Tight and regular packing motifs can lead to high π -orbital overlap and strong excitonic coupling between the molecules and as a result, the exciton can extend over several units. With increasing disorder in the system as depicted in Figure 1(c), the inter-molecular coupling is weak, and typically single-molecule Frenkel excitons can be observed.
In multicomponent materials, as shown in Figure 1(d), one can often find another type of exciton that extends over two molecules. Unlike an extended Frenkel exciton in ordered organic morphologies as in Figure 1(b), where both electron and hole part delocalize similarly over the molecular building blocks, the so-called charge-transfer (CT) exciton [14] is characterized by separation of the two charges on donor and acceptor parts. Due to the increased distance between electron and hole, the exciton binding energy is reduced compared to the Frenkel excitons in the bulk phase. As we will discuss in Section 3, the details of the conversion process between localized Frenkel and these bi-molecular CT excitons is a significant step, e.g. in the generation of free charges in organic solar cells.
We emphasize that the above distinction between inorganic and organic materials, as well as between Wannier-Mott, Frenkel, and charge-transfer excitons is simplified for the sake of a compact presentation and therefore far from exhaustive. In fact, there are many examples of materials with mixed characteristics. These include perovskites or perovskite-like structures with embedded molecules, such as CH 3 NH 3 PbX 3 with X = I, Br, Cl, in which a CH 3 NH þ 3 is surrounded by PbX 6 octahedra [15]. Other examples are organic-inorganic hybrids like metal or semiconductor nanoparticles functionalized by organic ligands, or dye-sensitized solar cell materials (e.g. TiO 2 with perylene-based dyes [16]). Furthermore, in soft conjugated polymers, electronic states and thereby also excitons can be subject to changes in localizaton due to dynamical or static variations in conformations of the πconjugated backbone, such as torsion angles between repeat units [17][18][19].

Role of exciton processes in devices
The nature of the different excitons as determined from material properties also has immediate consequences for dynamical processes involving the electron-hole pairs, and their exploitation in device applications [20].
While the binding energy of Wannier-Mott excitons in, e.g. silicon is comparable to the thermal energy and it is hence possible to generate free charge carriers after photon absorption in silicon solar cells, more complex processes need to be considered in disordered, organic materials. Figure 2 llustrates two examples for the vital role excitons and excitonic processes play in opto-electronic device applications based on organic heterostructures: organic light-emitting diodes (OLEDs, Figure 2(a)) and organic solar cells (Figure 2(b)).
An OLED device typically consists of several layers, each of which contributes a specific task to the overall device functionality. The schematic in Figure 2(a) shows a simple case, in which an emission layer is sandwiched between hole and electron transporting layers, respectively, both contacted to electrodes. The electronic processes needed for actual light emission to occur comprise injection of charge carriers from the electrodes into the transport layers, their drift-diffusion through the transport into the emission layer. Careful tuning of electron and hole transport level, as well as the energies of singlet and triplet excitons, is required to ensure that the excitons are formed and they efficiently emit light via fluorescence or phosphorescence. Significant effort is currently directed at optimizing materials for thermally activated delayed fluorescence, a process in which a molecule is initially in a non-emitting excited state before thermal energy of the surrounding allows it to change to an emissive state, or exciplex emission, i.e. emission from a bi-molecular charge-transfer exciton.
The active layer of an organic solar cell is usually a mix of donor and acceptor materials, and the interface between them, as depicted in Figure 2(b), plays a prominent role in the charge generation process. As mentioned before, the excitons created upon light absorption in either of the two layers are strongly bound with an exciton binding energy on the order of several tens ofeV. Thermal energy alone is therefore not sufficient to separate the excited electron and hole. Instead, the Frenkel excitons diffuse towards the interface at which they can, ideally, undergo a conversion process to a bi-molecular CT exciton with reduced binding energy. This conversion combined with electrostatic energy profiles near the interface, eventually allow the charges to separate and to transfer through the bulk material to the electrodes.
In both application examples, the device functionality is directly linked to the static (absorption/emission, recombination) or dynamic (diffusion, conversion, separation) processes involving excitons of various types. Also in more general cases, optimizing the device performance therefore targets to a significant extent the design of materials and material combinations with tailored and well-controlled exciton properties.
Ab initio modeling can be a very powerful tool in this rational design of materials. Instead of being descriptive, like hydrogenic models which omit explicit atomistic details of materials, it holds the potential of not just being predictive (energies and wave functions) but also allows to analyze and understand how the predicted properties emerge from atomistic detail. To fulfill this potential is, however, an extremely challenging task in general. Due to the two-particle nature of the excitation, popular one-electron methods such as Hartree-Fock or Density-Functional Theory (see Section 4.1 and Section 4.2.1) are not sufficient to properly describe the excitation process. Advanced -and therefore usually more computationally demanding -techniques are required instead. Each of these methods have their own advantages and disadvantages, arising either from the approximations made, their computational complexity, or the need to apply them to the vastly different material types as discussed above.

Overview of ab initio techniques
Ab initio, a Latin expression for 'from the beginning', refers to theories and computational methods that treat the many-electron problem using as inputs physical quantities only. The starting point is the (non-relativistic) many-body Hamiltonian of an atomic, molecular or solid system made out of M atoms and N electrons. If electronic and the nuclear degrees of freedom are decoupled and the nuclear motion is negligible (Born-Oppenheimer approximation [21]), the exact many-electron Hamiltonian is (in atomic units) Here, the first term denotes the kinetic energy operator of the electrons, the second term refers to the attractive electron-nuclei interactions (Z α is the nuclear charge of atom α), the third term describes the electron-electron repulsion, and the final term describes the repulsion among the nuclei of the system. Unless explicitly stated, throughout this review, we refer to electronic coordinates with lowercase (fr i g) while we use capital letters for nuclear coordinates (fR α g). We will interchangeably use the Dirac (bra-ket) and wavefunction notation for all the equations in the text. Once the Hamiltonian is known, the solution of the stationary Schrödinger equation yields the exact anti-symmetric wave function of the many-electron system. Analytic solutions are not available except for minimal cases.
Seeking for a solution numerically is also impractical. A naive attempt is to discretize the wave function on a spatial grid. Even when choosing a very coarse grid of 20 grid points in each direction, 20 3N values are needed. A water molecule has 10 electrons. Therefore, more than 10 39 values would be required. This is not feasible even on large computers. This opened the quest for methods that can lead to accurate results within a reasonable computational cost. The goal of ab initio methods is finding ingenious ways to numerically determine an approximate wave function. In the following we give a brief overview of different methods, distinguishing between approaches that work with an N-electron wave function and those that are based on the electron density. For a more exhaustive perspective of these methods and numerical techniques, we refer the reader to Refs. [22,23].

Many-electron wave-function approaches
One of the first methods employed to tackle the many-body problem of Eq. (1) is the Hartree-Fock (HF) method. It assumes that the exact N-body wave function of the system can be approximated by a single Slater determinant of single particle wave functions to enforce the antisymmetry property of a fermionic system: Invoking the variational principle, one can derive a set of N coupled equations for the N spin orbitals from a minimization of the expectation value of under the constraint of normalization of the effective single-particle functions. The resulting Hartree-Fock equations for the orbitals fϕ HF is the exact exchange operator defining the electron exchange energy due to the antisymmetry of the total N-electron wave function. Using the ϕ HF i in Eq. (3) yields the ground-state Slater determinant Ψ HF 0 . It has been observed that the HF method yields total energies (or derived properties) that can deviate significantly from experiments, as it includes no further quantummechanical effects apart from the exchange interaction from the antisymmetry requirement. The missing ingredient going beyond exchange is also referred to as correlation. It should be noted that the above is a treatment of the ground state, whereas excited states need special consideration, e.g. in a time-dependent HF (see related discussion in Section 4.2.2).
A number of approaches, called post-Hartree-Fock methods, have been devised to include electron correlation to the multi-electron wave function and also allow for the construction of excited state wave functions, as we summarize in the following.

Coupled cluster theory
Coupled cluster (CC) theory exploits the basic Hartree-Fock molecular orbital method and constructs multi-electron wave functions using the exponential cluster operator to account for electron correlation: with T ¼ P i¼1 T i the cluster operator and T i the i-th excitation operator (i ¼ 1 single excitation, i ¼ 2 double excitation and so on). This formulation provides in principle the exact solution to the time-independent Schrödinger equation Eq. (2) but makes the problem non-Hermitian and the obtained energy non-variational. Variationally optimized coupled cluster [24] overcomes this problem but can only be applied to small systems [25][26][27]. Another drawback is the computational cost of most coupledcluster implementations that makes these methods suitable only for small molecules, in general. CC methods scale at best as OðN p Þ, where p is a relatively high power (e.g. p ¼ 7 for coupled cluster including single and double excitations, with triples treated perturbatively) and development of a efficient techniques is an active field of research. CC has been successfully used for excited state calculations [28][29][30]. Wang and Berkelback [31] have recently shown that Equation-of-Motion Coupled-Cluster Theory can yield promising results for optical excitation energies, exciton binding energies, exciton dispersion relations, and exciton-phonon interaction energies of simple crystalline solids like Si, C, SiC, MgO, or LiF. Other attempts are made in this direction for ground-state and excited-state methods that combine CC and perturbation theory based on a partitioning of excitations that are internal or external to an active space [32][33][34][35].

Configuration interaction
Configuration interaction (CI) is another post-Hartree-Fock linear variational method for quantum chemical multi-electron systems. Mathematically, configuration describes the linear combination of Slater determinants used for the wave function. In order to account for electron correlation, CI uses a variational wave function that is a linear combination of configuration state functions (CSFs) built from spin orbital Slater determinants of the kind of Eq. (3) If the expansion includes all possible CSFs of the appropriate symmetry, then this is a full configuration interaction procedure which exactly solves the electronic Schrödinger equation within the space spanned by the oneparticle basis set. If only one spin orbital differs, we describe this as a single excitation determinant. If two spin orbitals differ it is a double excitation determinant and so on. Due to the long CPU time and large memory required for CI calculations, the method is limited to relatively small systems.

Quantum Monte Carlo
Quantum Monte Carlo (QMC) has been very successful in performing large-scale calculations for extended systems in quantum chemistry. The method relies on the stochastic estimation of the energy of a trial wave function according to For the calculation of excited-states, many refined QMC methods have been developed. Variational Monte Carlo (VMC) as presented by Zhao et al. [36] allows for the application of a rigorous variational principle to both ground and excited state wave functions. Diffusion Quantum Monte Carlo (DQMC) is very promising for applications to condensed-matter because it explicitly includes electron-electron correlation effects, and it scales reasonably well with system size as shown in the case of a silicon crystal in [37]. These QMC methods offer an accurate means of probing both the groundand excited-state properties of atoms, molecules, and solids from first principles and with system-size scaling as OðN 3 Þ, albeit with a large prefactor. The methods have been thoroughly reviewed by Hunt et al. [38] discussing isolated molecules (anthracene, tetracyanoethylene, benzothiazole and boron trifluoride) and three-dimensional systems (diamond, silicon, cubic boron nitride) and free-standing monolayer phosphorene. While describing excited states, a major reason of failure is the intrinsic difficulty of maintaining orthogonality with the lower-lying states when the targeted many-body excited state is being represented stochastically in an imaginary-time projection method. Auxiliary-field Quantum Monte Carlo (AFQMC) as discussed in [39] offers a new framework for addressing this difficulty and for performing excited state calculations in solid. Ma et al. studied the fundamental gap of prototypical semiconductors, Si and diamond, and of the more challenging wurtzite ZnO crystal being in good agreement with GW calculations (see Section 4.2.5) and experiment, offering a non-perturbative and free of empirical parameters methods for correlated materials.

Mean-field theories
Completely different approaches for the study of excitons are the so-called mean-field theories. In a nutshell, an effective Hamiltonian capable of including many-body effects in response to an external excitation (such as incident light) can be used to describe charged (quasiparticles) and neutral excitations (excitons). Mean-field theories can come in different flavors depending on which quantity under perturbation is taken as a reference, either the time-dependent electron density nðr; tÞ (see Section 4.2.2) or the one-particle Green's function Gðr; t; r 0 ; t 0 Þ (see Section 4.2.5). The first choice leads to a time-dependent extension of Density Functional Theory (TD-DFT), while the second choice stems from the many-body perturbation theory. This theory formally describes the physics of charged excitations (i.e. electron addition and removal energies), using, e.g. Hedin's GW approximation for the electron self-energy and neutral excitations (e.g. optical and energy-loss spectra) via the solution of the Bethe-Salpeter equation (BSE). Both theories proved to be quite successful at the calculation of the absorption spectra of a large variety of systems: insulators, semiconductors, atoms, clusters, surfaces, or polymers. Despite being different theories, TD-DFT and GW-BSE share some common ingredients. In the following, a small summary of the main equations for both theories is presented. For the sake of simplicity and brevity, we also restrict the discussion to the non-periodic formulation of DFT, TD-DFT, and GW-BSE. For more details about these theories, also in the case of periodic potentials, we direct to the work of Onida et al. [40] and Rohlfing and Louie [41].

Density functional theory
Density Functional Theory (DFT) [42] is the base ingredient for both TD-DFT and GW-BSE approaches since it describes the electronic ground state. The ground-state energy of a system of interacting electrons in an external potential can be written as a functional of the ground-state electronic density E½n 0 ðrÞ�. Compared to the methods in [sec:wfmethods] Section sec:wfmethods this approach is particularly appealing since it does not rely on a complete knowledge of the explicit N -electron wave function but only on the electronic density. For a closed-shell system of N electrons, following the formulation of Kohn and Sham [43], ground-state properties are derived from the solutions of the Kohn- Here b V xc ½n 0 � is the exchange-correlation (xc) potential, formally containing all many-body effects. Although the theory is exact, the xc contribution to the energy functional, and hence also the potential, is not known and has to be approximated in practical implementations, leading to a zoo of approximate exchange-correlation functionals. A thorough discussion of specific functional choices is beyond the scope of this work and refer the reader to Refs. [44,45]. Instead, we briefly mention three main types of approximations commonly used among the DFT community. The simplest of these approximation is the local density approximation (LDA). The assumption behind this approximation is that the charge density of the system, not homogeneous overall, is locally similar to the one of the homogeneous electron gas, whose exchange-correlation energy is known [43]. An improvement upon the LDA can be obtained by semi-local Generalized Gradient Approximation (GGA) functionals [46,47]. These depend not just on the value of the density at a point (as in the LDA case) but also on its gradient. The last popular approximation are called hybrid functionals [48,49]. Hybrid functionals are based on the ansatz that the exact exchange energy is situated between the GGA exchange energy functional and the Hartree-Fock exchange integral. In these, the Hartree-Fock exchange integral is mixed with GGA exchange functionals at a constant ratio.
DFT proved to be the game-changer in all those calculations that requires the knowledge of the ground-state density and total-energy (atomic structure, lattice parameters, total energy, phase stability, electronic density, elastic constants, and phonon frequencies). On the other hand, KS eigenvalues cannot be considered as true electronic energy levels, that is, as experimentally measurable electron-addition or electron-removal energies in direct or inverse photoemission experiments. The most noticeable drawback is that the HOMO-LUMO gaps in molecular systems, or the band gap in extended semiconductors or insulators, are dramatically underestimated [50,51] (see Figure 3). On top of that, DFT cannot take into account the Figure 3. Schematics of the theoretical steps for a description of exciton formation within MBGFT. When incident light of energy � hω is absorbed by the material, an electron is promoted to the unoccupied manifold (blue curve) leaving an hole behind in the occupied one (red curve). These two charged particles, under some circumstances depending on the dielectric properties of the material, can be independent, hence non-interacting. Despite its success in describing ground-state properties, KS-DFT (a) gives wrong quantitative results for charged excited states (e.g. underestimate of the band gap). A correct description can be obtained using a quasiparticle picture (b) in which real excitation is formally treated including many-body effects embodying the electron/hole states with the perturbation of its own surrounding. Band gaps are thus improved. These two oppositely charged quasiparticles can interact forming a bound state known as an exciton (c). They are attracted via a screened Coulomb potential with v c the bare Coulomb interaction and � � the dielectric function of the material. This neutral excited state has a lower energy (the missing energy is the so-called exciton binding energy). Exciton formation is highly dependent on the material. In fact this can be hindered (thus leaving only free charged carriers) by different dielectric responses in different materials.
coupling of electrons and holes being an independent particle theory. These intrinsically excited state quantities need special treatment.

Time-dependent density functional theory
Neutral excitations of the DFT ground-state can be described via timedependent density functional theory (TD-DFT) [52]. In analogy to the traditional time-independent Kohn-Sham scheme, all exchange and correlation effects in TD-DFT are collected in the time-dependent V xc ðr; tÞ potential. In the linear response formulation of TD-DFT [53], the susceptibility χðx; x 0 Þ (with x ¼ ðr; tÞ and x 0 ¼ ðr 0 ; t 0 Þ space-time variables) is defined as the linear response kernel between the variations of the electron density δnðxÞ ¼ nðxÞ À n 0 ðxÞ with respect to a local external perturbation Uðx 0 Þ: In its spectral representation (Laplace transform) this is usually connected to the independent-particle susceptibility χ 0 (a response function whose poles are the KS energies) via a Dyson-like equation: χðr; r 0 ; ωÞ À 1 ¼ χ 0 ðr; r 0 ; ωÞ À 1 þ K DFT ðr; r 0 ; ωÞ ¼ χ 0 ðr; r 0 ; ωÞ À 1 þ f Hxc ðr; r 0 ; ωÞ; with the excited-state wave function being A S vc stand for the resonant (occupied to unoccupied) components whereas B S vc the provide the anti-resonant (unoccupied to occupied) contributions. Then, the matrix elements of H res and K are calculated as where κ ¼ 2 ð0Þ for spin singlet (triplet) excitations, and The exact time-dependent exchange-correlation action functional (from which V xc and f xc are derived) is not known, and approximations have to be introduced to perform numerical calculations on real systems. The most used approximation is the so-called Adiabatic Local Density Approximation. This approximation implies frequency-independent and real xc kernels in linear response calculations. The main advantage of TD-DFT is its efficiency (all interactions are local) and scales better than other computational methods. In fact OðN 2 Þ and OðNÞ implementations can be found [55]. Excited-state forces are available [56,57] and TD-DFT is easily combined with classical Molecular Mechanics models for mixed quantumclassical simulations (QM/MM) [58]. The major drawback of TD-DFT is its sensitivity on the choice of the exchange-correlation functional. Using the wrong functional may lead to wrong results. For instance, ALDA-TD-DFT can yield poor results because of the absence of an attractive electron-hole term in the present semilocal exchange-correlation functionals [59], prompting the research of a completely new class of functionals. Remedies for some aspects of these problems have been suggested, among them orbital-dependent density functionals, such as the exact-exchange optimized effective potential method. An approach that has recently attracted considerable interest is hybrid functionals mixing the nonlocal Hartree-Fock exchange term to the otherwise semilocal functional. In the letter by Sharma et al. [60] instead, a new parameter-free approximation for f xc is discussed. This kernel gives accurate results with the computational cost of ALDA. Another possibility is the nanoquanta kernel by Sottile et al. [61], derived from the four-point Bethe-Salpeter kernel (see Section 4.2.6). This is very accurate but has the drawback of being nearly as computationally demanding as solving the BSE itself. The long-range correction (LRC) kernel [62] instead has a particularly simple form in reciprocal space which limits its computational cost. This kernel produces the desired excitonic peak but depends on the choice of the parameter, which turns out to be strongly material-dependent, thereby limiting the predictiveness of this approximation. A well known problem in TD-DFT is the correct detection of charge-transfer states. The reason for this failure is the wrong decay of the exchange-correlation potential with increasing electron-electron distance. Hybrid [63] and double hybrid Functional [64] have been proposed as a solution to this problem.

Many-body Green's functions theory
As stated before, KS-DFT energies and wave functions are not suitable for describing excited electrons and holes. A formal description of electrons and holes resulting from an external excitation can be given within the framework of Many-body Green's Function Theory (MBGFT), specifically with the GW-Bethe-Salpeter (GW-BSE) scheme. We can anticipate that the GW spectrum provides accurate occupied/virtual energy levels, that can be directly compared to experimental electron addition or removal energies, while the BSE scheme adds the electron-hole interaction leading to a reduction of the optical gap as compared to the (photoemission) HOMO-LUMO gap. This is sketched in Figure 3.

Charged excitations: quasiparticle picture
Excitations that add ðN ! N þ 1Þ or remove ðN ! N À 1Þ an electron from the system, referred to as quasiparticles (QPs), are determined by the one-particle Green's function [65] G 1 ðr 1 ; t 1 ; r 2 ; t 2 Þ ¼ À ihN; 0j b T ðψðr 1 ; t 1 Þψðr 2 ; t 2 ÞÞ N; 0 j i; where b T is the time-ordering operator, and b ψ and ψ are the annihilation and the creation electron field operator, respectively. It describes the propagation of electrons and holes (lack of electrons) measuring the amplitude of probability of finding, on top of the ground-state, an electron in (r 1 ; t 1 ) that was previously introduced in (r 2 ; t 2 ), while for t 1 < t 2 the propagation of an hole. The frequency-dependent spectral representation of the Green's function has poles at the charged excitation energies: where the index n runs over all the charged excitations. In other words its poles are proper addition/removal energies of the N-electron system. This Green's function obeys a Dyson-type equation of motion, which reads in spectral representation as where the electron self-energy operator b �ðEÞ contains the exchangecorrelation effects. This equation is part of a closed set of coupled equations, known as Hedin's equations [65,66].

GW approximation
An approximate solution to this system is provided by the GW approximation, in which the self-energy takes the form �ðr; r 0 ; ωÞ ¼ i 2π i.e. it is a convolution of G 1 with the screened Coulomb interaction where v c ðr; r 0 Þ ¼ jr À r 0 j À 1 is the bare Coulomb interaction and � � À 1 ðr; r 0 ; ωÞ is the inverse dielectric function calculated in the Random-Phase Approximation (RPA) [67].
Using this approximation, Eq. (20) is converted into a Dyson-type equation of motion for the quasiparticles (i.e. the QP electron and hole states) [68,69]: where ε i QP are the one-particle excitation energies of the system, and jϕ QP i � are the QP wave functions. As the self-energy is energy-dependent, and thus depends on ε QP i , the solution of Eq. (22) must be found self-consistently. Constructing both G and W based on the KS eigenvalues and eigenfunctions is known as a 'one-shot' G 0 W 0 calculation. Alternatively, it is possible to use updated QP energies in G and W until eigenvlaue self-consistency is reached (evGW) [70][71][72].
The frequency integration in Eq. (21) is one of the major difficulties for GW calculations since both functions that are integrated have poles infinitesimally above and below the real frequency axis. A numerical integration is not feasible because of the instability of the integrand to be evaluated in regions in which it is ill-behaved. Approximations and exact alternatives are available.
The simplest way to calculate the frequency integral is to approximate the frequency dependence of the dielectric function � � À 1 and thus the screened Coulomb interaction by a plasmon pole model (PPM) [67,[73][74][75]. The PPM approximation takes advantage of the fact that � � À 1 is usually dominated by a pole at the plasma frequency ω p . This pole corresponds to a collective charge-neutral excitation (a plasmon) in the material. Assuming that only one plasmon branch is excited, the shape of � � À 1 can be modeled by a singlepole function. A comparison of different PPM methods has been done by Larson et al. [76] and Stankovski et al. [77].
A full-frequency integration technique for the calculation of �ðωÞ is the contour deformation (CD) approach [78][79][80][81][82][83]. In the CD approach, the realfrequency integration is carried out by using the contour integral extending the integrand to the complex plane, the numerically unstable integration along the real-frequency axis, where the poles of G and W are located, is avoided. Another method that enables an integration over the fullfrequency range is the Analytic Continuation (AC). The AC technique exploits the fact that the integral of the self-energy along the imaginary frequency axis is smooth and easy to evaluate, unlike the integral along the real-frequency axis. The self-energy is first calculated for a set of imaginary frequencies iω and then continued to the real-frequency axis by fitting the matrix elements �ðiωÞ to a multipole model (either the so-called 2-polemodel [84,85] or the Padé approximant [86][87][88]). A fully-analytical approach (FAA) [89][90][91][92] can be used to perform the integration in Eq. (21) . In fact the screened Coulomb interaction can be written as a Dysonlike equation WðωÞ ¼ v þ v � χðωÞ � v with χðωÞ is the frequency-dependent susceptibility. As in the TD-DFT case, poles of this function can be found solving a Casida-like eigenvalue problem similar to the one of Eq. (11) . In this case, the exchange-correlation kernel f xc is omitted (otherwise it would be TD-DFT). The self-energy integral can then be solved analytically and a closed expression for �ðωÞ is obtained. The FAA is the computationally most expensive method. Solving the eigenvalue problem to obtain the poles of χ is an OðN 6 Þ step. The scaling of the CD and AC approach is generally lower, but depends on the details of the implementation. The PPM is computationally the most efficient method, because the dielectric function � � À 1 used to compute W has to be calculated only at a few frequency points to determine the parameters of the PPM. By design, the fully analytic approach is in principle the most exact one since it is parameter-free. The same accuracy can already be achieved with the CD using a moderately-sized numerical integration grid for the imaginary frequency term [83]. In the AC approach the self-energy is calculated on the imaginary frequency axis, which is fairly featureless. The accuracy of the AC approach depends on the features of the self-energy on the real axis and on the flexibility of the model function, which continues the self-energy to real frequencies. For this reason the AC is likely to fail for deeper states. In the PPM approximation the accuracy is determined by the chosen model function and generally difficult to estimate. Further information about the details of the GW approximation can be found in the work of Golze et al. [72].

Bethe-Salpeter equation
Neutral excitations with a conserved number of electrons and a change in their configuration S ( N; 0 j i ! N; S j i) can be described within MBGFT with some analogy to the TD-DFT case. Here, the approach relies on the four-point response function (also known as two-body Green's function) [93], which is formally defined as the variation of G with respect to an external perturbation U Similarly to TD-DFT this response-function can be written in a Dyson-like equation Þ . This kernel is different to the TD-DFT one and it is more cumbersome to compute. Again, this can be turned into an effective Hamiltonian problem in the same form as Eq. (11) . In this case however KS energies and wave functions are substituted with the QP ones. If the GW approximation is used for the selfenergy we have in lieu of Eq. (15), Eq. (16) and Eq. (17) of TD-DFT the following diagonal, exchange, and direct term for BSE The direct interaction K d contains the attractive, but screened, interaction W between electron and hole. This interaction is responsible for the binding of the electron-hole pair (compare to Eq. (15) in TD-DFT). In systems for which the elements of the off-diagonal blocks of Eq. (11) are negligible (K ' 0), it is possible to use the Tamm-Dancoff Approximation (TDA) [94]. On the other hand, the coupling between resonant and anti-resonant parts is significant, and its neglect can cause deviations of several [0.1] eV from results obtained with the full approach [95,96]. More computational cost can be saved employing the so-called static approximation for Eq. (27) in which the static limit of the screened Coulomb interaction is taken (WðωÞ ! Wð0Þ). This approach is sufficient most of the time to obtain meaningful results. If an improvement is needed a perturbative approach can be used [97,98]. GW-BSE is expensive (OðN 4 Þ) but less than QMC, CC, etc. A very remarkable success of the BSE formalism lies in the description of chargetransfer (CT) excitations, a notoriously difficult problem for TD-DFT adopting standard (semi)local functionals as discussed before. The screened Coulomb potential in the BSE kernel includes such long-range (nonlocal) electron-hole interactions, including in environments (solvents, molecular solid, etc.) where the screening reduces the long-range electron-hole interactions [99]. Overall, GW-BSE has proven to be very successful at treating different types of excitons equally well [100][101][102][103][104][105][106][107][108][109][110][111][112][113].
A significant limitation of the BSE formalism, as compared to TD-DFT, lies in the lack of analytic nuclear gradients for both the ground and excited states, preventing efficient studies of many key excited-state processes. At the moment only one study introduced an approximate analytic form of the excited-state BSE gradients for small molecules (CO and NH 3 ) has been published [114]. A more recent evaluation of the geometry optimization of small molecules with numerical gradients has shown that in principle reliable excited state geometries can be obtained with GW-BSE [115].

Perfect crystals
Initial applications of GW-BSE have been reported for perfect inorganic crystals, such as the elementary semiconductors silicon (Si) and germanium (Ge), or insulators such as magnesium oxide (MgO) or lithium fluoride (LiF). The atomic structure of these materials can be represented by small unit cells with only a few atoms, having made these early studies tractable. Current implementations [116][117][118] with periodic boundary conditions often employ a plane-wave basis, a natural choice since electrons experience a periodic potential and thus they can be described as Bloch states [119]. Hybrid schemes combining localized orbitals for the real space wave function representation and plane waves for the evaluation of Coulomb interaction in Fourier space also exist, as implemented in the CP2K code [120]. All-electron linearized augmented plane wave (LAPW) approaches [113] enable the computation of core excitations on the same footing as optical ones [121]. Figure 4 shows a typical result of a periodic GW-BSE calculation for small crystals. Rohlfing and Louie [41] calculated the optical absorption spectra of bulk (a) GaAs and (b) LiF with particular focus on the effect of the inclusion of electron-hole interaction via the BSE. Whereas the absorption spectra resulting from the free-quasiparticle transitions show systematic deviations from experiments in both systems, the inclusion of the electron-hole interaction leads to a good agreement with the measured data. In particular, for the more ionic LiF, a pronounced exciton peak can be observed in Figure 4 (b) that is missing in the independent particle picture. Figure 4(c) shows visualizations of the electronic part of the two-particle wave function 1 jψ S ðr e ; r h Þj 2 with the hole coordinate r h fixed on top of the central fluorine atom. One can see the Wannier-Mott character of the exciton as its wave functions extend over several of the primitive unit cells.
Recent developments [122,123] and increased computing power have allowed to treat larger systems, such as reconstructed surfaces, defect structures, or in general materials with more complex structures. Yet, studies of the exciton properties in some of the currently most intriguing systems are still extremely difficult to perform. An example of a such a class of materials are perovskites of the type OMX 3 where O ¼ Organic, M ¼ Metal and X ¼ Halide [124]. They are attractive materials due to their simple production procedure and high efficiency of above 20% in photovoltaic cells. Optical gaps in the range of 1:6-3:1eV have been reported. Actual details of the atomic structure of perovskites do however show a significant temperature dependence, which needs to be accounted for in ab initio modeling approaches. Of particular interest is the fact that high efficiencies have been reported, even though the exciton binding energy is much higher than the thermal energy at room temperature (as in organic materials). Naturally, using ab initio GW-BSE modeling to resolve this issue is appealing but theoretically challenging. Structure prediction can be complicated and must be carefully performed due to temperature dependence. Relevant structures are then described in a large unit or supercells, translating into significant computational costs. Furthermore, the presence of heavy elements (e.g. Pb) requires the inclusion of spin-orbit coupling.
From the point of view of modeling, organic crystals share many complications arising from large unit cells, temperature-dependent structures, and the mixed character of the excitons with perovskites. Typical molecules forming these crystals are flat, aromatic molecules such as the polyacenes, in particular naphthalene, anthracene, tetracene, and pentacene, as well as pyrene, perylene, and similar compounds. As we mentioned before, the attractive force between the neutral and nonpolar molecules is provided by comparatively weak van-der-Waals interactions. The crystalline phases can be described in the same way as the inorganic crystals by well-defined lattice structures and periodic boundary conditions. However, the primitive cells typically contain a few molecules in distinct packing motifs and are therefore relatively large. While in principle the same types of periodic GW-BSE calculations as mentioned before for inorganics can be employed, they face the challenge of having to describe strongly localized (with respect to the primitive cell) excitations. High energy cutoffs are required in plane wave implementations and large calculation costs are the consequence.
For example, Cudazzo et al. [106] investigated the excitonic properties of picene and pentacene crystals, which show remarkable differences in their optical spectra despite otherwise similar structural and electronic properties. They employed a periodic, plane-wave implementation of GW-BSE in the TDA based on an LDA ground state. They could show that a difference in localization of the excitons gives rise to the different absorption spectra and related this to the competition between the direct and exchange interactions in the BSE. Figure  5 shows visualizations of the wave functions of the lowest-energy singlet (a,b) and triplet (c,d) excitons in picene and pentacene, respectively. For a fixed hole position, the isosurfaces represent the electronic part of jψ S ðr e ; r h Þj 2 of the singlet in Figure 5(a,b) and of the triplet in Figure 5(c,d), respectively. The singlet exciton in picene is localized on a single molecule, while it extends over several molecules in pentacene. In contrast, the triplet exciton appears to be localized in both cases. Due to the structure of the BSE Hamiltonian, it is clear that the delocalization in pentacene is driven by the repulsive exchange interaction.

Disordered molecular materials
Disordered systems cannot be described efficiently by repeated small unit cells. Structure simulations of, e.g. amorphous materials are typically based on either ab initio or (more likely) classical Molecular Dynamics and can contain several hundred to a few thousand of molecules. The investigations of the exciton properties for systems of this kind cannot use the same techniques as those in Section 5. As a consequence, the study of electronic excitations in complex disordered systems remains a challenge for computational methods. Typical solution strategies on the other hand exploit that excitons strongly localize on single molecules due to the disorder of the large-scale morphology, and thus the interaction of this active exciton with the surrounding can be described with lower-level (i.e, less expensive) methods. This leads to various strategies to split the system into an 'active' portion of space, where the highest desired level of Quantum Mechanical (QM) treatment is performed, and an 'environmental' surrounding medium (a solvent, a metallic nanoparticle, a disordered (polymeric phase), etc.) that is considered at a lower level of theory. This lower level can be either a simplified QM model (QM/QM') (Dvorak et al. [125] explored embedding of wave function theories with Green's functions as Manby et al. [126] proposed a generalized and flexible QM-in-QM embedding scheme used to explore excited states [127,128]), a molecular mechanics (MM) approach (QM/MM) [97,[129][130][131], or a continuum model [132] representing the environment as a structure-less material having realistic macroscopic dielectric properties. The effects of the embedding on the quasiparticle and electron-hole energies is schematically shown in Figure 6.
Specifically using GW-BSE as a high-level QM method, Duchemin et al. [132] have incorporated the polarizable continuum model (PCM) into the calculation of the screened Coulomb interaction W ¼ � � À 1 v c . In this way W accounts for both the polarizability of the solute and that of the solvent when combined with the PCM. This has been coupled to the BSE [133] to analyze the contributions to the BSE dynamic shift originating from the solvent-induced renormalization of the ionization potential and electron affinity, namely the variation of the GW HOMO-LUMO gap in the water solvent, and the reduction of the electron-hole interaction for excitation in acrolein, the CT excitation in the benzene-TCNE complex, and the (planar) PNA mixed excitation.
A prominent class of atomistic multi-method embedding approaches renders the environment by classical multipole potentials, as in classical Molecular Dynamics. We illustrate how such a choice of the embedding method affects GW-BSE excitation energies using an example from our work on a bulk mixture of a rubrene guest donor in a host matrix of fullerene C 60 acceptor [131]. Here, the QM region contains two molecules forming a donor-acceptor charge-transfer exciplex. GW-BSE calculations on such isolated molecular structures are efficiently performed with an implementation based on atom-centered Gaussian-type orbitals, and is available in the VOTCA-XTP package [131,134], among others [91]. Figure 7 shows two selected configurations of the exciplex as found in the bulk MD morphology: one in which the fullerene is close to one of the rubrene's phenyl rings (CT P , (a)) and one in which it is close to the anthracene core (CT A , (b)). Such differences of mutual position and orientation can have a massive influence on the nature of the CT excitons in particular, as can be seen in the energy level diagrams in Figure 7(c,d), where we compare also the effect of different levels of theory for the embedding MM region. This region includes the atoms of molecules within a pre-defined spherical cutoff radius around the exciplex, as indicated by small spheres in the zoomed-in views in Figure 7. We consider the difference between the case when the MM atoms are reduced to static point charges (GW-BSE/MM s ), and when they are additionally endowed with dipole polarizabilities according to Thole's model [135] (GW-BSE/MM p ). This atomistic polarizable embedding evaluates the self-consistent reaction field to the excitation and feeds that back into the Hamiltonian via an Figure 6. Schematic representation of the state-specific correction due to embdedding effects within the GW-BSE/MM formalism. The GW energy levels are first renormalized by an amount given by the polarization energy for electrons (occupied manifold) and holes (virtual manifold), respectively. Further, the strength of the screened electron-hole interaction, that here we label as W eÀ h , is also reduced due to the additional screening provided by the environment. updated external potential. For a full account, see Refs. [17,131,134]. Figure  7(c,d)shows that the static embedding has a nearly negligible effect on the two Frenkel excitons as well as the CT excitons. This is, however, specific to this particular material combination and morphology, as the C 60 host material has practically no local electric fields in the point-charge representation. Upon the inclusion of polarization in the environment, the energy of the integer CT exciton (CT P ) is significantly lowered. In the CT A configuration, which shows a strong interaction between the two molecules, energetic stabilization is less pronounced but one registers a noticeable increase in the amount of charge transfer upon excitation. The results clearly emphasize the importance of including polarization effects in some form into the embedding method. Using GW-BSE as the QM method is particularly attractive as it can predict in particular the CT excitons without special adjustments, thanks to the balance between direct and exchange interactions in the BSE Hamiltonian.
The bulk morphology with a dilute guest-host mixture of rubrene and C 60 lends itself to a cutoff-based inclusion of the environment due to the lack of local electric fields and the lack of partial structural order. We consider now a related example of another fullerene-based donor-acceptor morphology with interfacial order. An interface between C 60 and DCV5T-Me (3,3), a thiophene derivative with internal acceptor-donor-acceptor architecture [136], is created in MD simulations from crystal structures of the two compounds and then relaxed at room temperature. Figure 8(a), we show for an example DCV5T-Me(3,3)/C 60 exciplex at the interface and the dependence of the Frenkel (FE) and CT exciton energies and the amount of charge transfer on the choice of the cutoff radius in the GW-BSE/MM s setup. Clearly, the CT exciton energies are not even nearly converged for a large cutoff of 40 nm. In the calculation, this behavior can be traced back to an interplay between the partial 2D order of the system with the internal electrostatics of the DCV5T-Me(3,3) molecule, which are dominated by its quadrupole moment. As in the case of the rubrene and fullerene mixture in the previous example, the localized excitations are not significantly affected, due to the much smaller change in molecular dipole moment upon excitation compared to the CT excitons. Although here, unlike in the dilute blend Figure 7, the Frenkel exciton energy exhibits some sensitivity up to a cutoff of 6 nm. When combined with polarization effects in the MM region, the GW-BSE/MM p setup yields an energy level alignment as indicated in Figure  8(b), that is, the CT exciton energies are higher than the Frenkel exciton energy. Such an alignment would not appear very favorable for generation of free charges at such an interface -a process that is known to occur efficiently in photovoltaic devices. This mismatch of the ab initio model with the experimental observation is attributed to the 'unconverged' electrostatic contribution as in Figure 8(a) and highlights emphatically the need for an efficient inclusion of long-range electrostatic interactions in a QM/MM setup for systems with reduced dimensionality and/or partial molecular order. This can be achieved, for instance, with suitably adapted Ewald summation techniques [137].

Biomaterials
Soft biological matter comes in a multitude of forms, and excitonic effects play a variety of roles in them, posing unique challenges for ab initio modeling that stem from the structural complexity of the materials. In contrast to the disordered molecular materials discussed in Section 6, biomaterials often consist of polymeric molecules in liquid solvents environments that can have intricate conformational dynamics. We will consider two examples in the following: the stabilization of a CT exciton in aqueous DNA, and time-resolved emission from polarity sensitive dyes.
In Figure 9(a), we show the properties of charge-transfer excitons between two adenine bases in a double-strand of DNA in an environment of water and solvated ions [109,134], as obtained from GW-BSE with static and polarizable embedding as discussed before. Interest in these systems is related to the experimental observation in an adenine 20-mer [138] of broad low energy absorption in the region of 3:6-4:2teV, below the main adenine absorption peak at [4.8] eV. The low energy absorption is attributed to the formation of CT excitons, which can have significant implications for UV processes causing damage to DNA. A few points are noteworthy about the calculation results: first, upon inclusion of a static environment, the localization of the electron and hole contribution to the charge-transfer exciton switches between the two involved adenine bases. This is accompanied by an energetic stabilization of [0.44] eV. Additional inclusion of polarization not only reduces the excitation energy by a further [0.54] eV and also transforms the partial CT exciton into an integer CT exciton. The strong effect of polarization can be traced back to the aqueous environment, highlighting the nanoscale insights one can obtain from the combination of an accurate ab initio method like GW-BSE with environment models. At the same time, the CT exciton energy is still too high by several [0.1] eV compared to the experimental expectation.
One possible reason for this discrepancy lies in the fact that the polarization included in the GW-BSE/MM p calculation is a purely electronic response, i.e. no structural relaxation of the aqueous environment is included in the model. The process of combined relaxation of the excited state and the solvent environment is exploited in polarity-sensitive dyes, which can be used, e.g. as spectroscopic probes for studies of biological systems. A prototypical dye of this kind is prodan, as seen in Figure 9. In such a molecule, the exciton created upon photo absorption has internal charge-transfer character [139], which manifests itself in a large change in the molecular dipole moment. Depending on their polarity, solvent molecules then react to this enlarged dipole by structural rearrangement as part of the polarization response. The exciton can, in turn, react to the changed external electric field, and the process repeats until an equilibrium is reached. From the point of view of an ab initio model, this situation is challenging because of the need for a technique capable of accurately capturing the details of the internal CT excitation, and the combined electronic and structural relaxation of the solvent-solute system. While the former is a strong point for GW-BSE, as discussed before, the latter is difficult for a lack of analytic force expressions. Motivated by the realization that the main effect of the relaxing CT exciton is its change in dipole moment (or the electrostatic potential felt by solvent molecules), Baral et al. [140] have developed an iterative scheme in which the structural relaxation of a prodan molecule in solution is modeled using classical MD with the dye's partial charges being updated at certain time steps by those extracted from a GW-BSE/MM s calculation. With this simplified scheme, they have been able to simulate time-resolved emission spectra as shown in Figure 9(b) in very good agreement with experiments for a range of solvents of different polarity. Similarly, the Stokes shifts could be predicted at good accuracy as reported in Figure 9(c). Based on these promising findings, they proceeded in investigating the fluorescent emission properties of laurdan [141,142] (a modification or prodan) in lipid bilayer membranes. An example of such a structure with highlighted embedded laurdan molecules is shown in Figure 9(d). It is known that the emission is sensitive to membrane properties like membrane order and packing. In particular, spectra are about 50 nm blue-shifted going from liquid disordered phase to gel phase or cholesterol-rich liquid-ordered phase, linked to the higher exposure to water solvent and the increased polarization response. This behavior could not be seen in the results of the same iterative GW-BSE/MM s plus MD procedure as employed for prodan, as the underlying MD simulations of the input structure before the exciton calculation does not accurately predict the location of the dye inside the membrane. When sampling of the chromophore position is biased to select only configurations in which laurdan is exposed to water, the experimental situation is recovered. These observations emphasize that for an accurate ab initio modeling of excitons in materials, a reliable prediction of structural properties is essential and must not be neglected.

Summary and outlook
In summary, we have discussed in this review the usefulness of the GW and Bethe-Salpeter equation approach for the ab initio modeling of excitons in different material types, from perfect crystals to biomaterials, from hard to soft matter. Examples were chosen to highlight the capabilities of the method as well as its current limitations. Future directions of research and method development should aim at overcoming these limitations, some of which have already been mentioned, including the need for efficient analytic force expressions or the QM/MM embedding including long-range interactions.
In particular concerning the hybrid QM/MM setups, several so far unnamed aspects deserve attention: for instance, the current implementation of the polarized MM relies on the Thole model of distributed atomic polarizabilities. This model is known to have problems for molecules with large anisotropies in their polarizability tensor, e.g. in extended π-systems [143], in which it is difficult to parametrize the atom polarizabilities to reproduce tensor components in-and out-of plane accurately at the same time. Similarly, it is not well-suited to describe polarization induced charge transfer or hydrogen bonding, which might both be important at the boundary between QM and MM regions. The Drude model or chargeequilibration approaches [144,145] could be superior in this regard.
Along similar lines, the use of a polarized MM model in the selfconsistent reaction field approach renders the coupled GW-BSE/MM p Hamiltonian state-dependent. This implies the need to efficient statetracking methods in the iterative solution of the coupled problem, as one has to reliably identify the state of interest in the spectrum of Eq. (27) . This can be achieved with characteristics of the exciton, such as its oscillator strength or the amount of charge transfer, or more detailed information from the overlap of wave functions from different iterations or a suitably defined measure of density matrix differences. All of these procedures have difficulties in tracking states near level crossings, which may happen frequently in dense spectra as the ones shown, e.g. in Figure 7(c,d).