Exciton States of Azobenzene Aggregates: A First‐Principles Study

Interaction between azobenzene‐containing molecules in supramolecular structures or self‐assembled monolayers (SAMs) results in the formation of molecular exciton states. These states determine photophysical and photochemical processes in such assemblies. Here, using first‐principles quantum chemical calculations, optical spectra and exciton delocalization of the exciton states in model clusters of azobenzene molecules are studied. Specifically, 1D linear chains and 2D SAM‐like arrangements are considered, and the exciton states are computed by means of time‐dependent long‐range corrected density functional theory (TD‐lc‐DFT) and ab initio configuration interaction singles (CIS), for clusters including up to 18 azobenzene molecules. The nature of the exciton states is analyzed using transition density matrix analysis. In addition, a connection to periodic systems is made applying the Bethe–Salpeter equation (BSE)/Green's function many‐body perturbation theory ( GW$GW$ ) approach to a selected system. It is found that the brightest excitons are dominated by local excitations. The energetic location of charge transfer states in the electronic spectra of aggregates depends to a large extent on a given method and distance between nearest neighbors. Furthermore, it is analized how an excitonic delocalization pattern changes with varying molecular orientation in the unit cell of SAMs.


Introduction
Azobenzene is a molecular photoswitch of trans ↔ cis type. Shifting between the isomers with light enables various transformations and control over numerous properties of systems ranging from macromolecules [1][2][3] through interfaces [4,5] to solids. [6][7][8] In (soft) condensed matter systems, such as azobenzene-containing self-assembled monolayers (SAMs), [9,10] micelles, [11][12][13] liquid . The rectangular unit cell (black outline) with lattice parameters a = 3.5 Å and b = 6.0 Å contains one molecule (marked by an arrow). Gray outline shows a smaller model, dodecamer, which was also studied. d) Orientation of an azobenzene unit with respect to the surface (the latter shown as gray bar) as determined in ref. [9]. Red arrow is the normal to the surface, green arrow is the normal to the molecular plane (the plane of the phenyl ring, marked by light green color), and blue arrow is the molecular axis defined to be parallel to the C−H(4') bond of the upper phenyl ring. The orientation is described by the following angles: = 76 • , = 30 • (the tilt angle), = 61 • (the twist angle). Note that at = 0 • the green arrow lies in the plane spanned by the red and blue arrows. e) SAM model constructed from the molecules oriented as shown in (d). [9] The rectangular unit cell (black outline) with lattice parameters a = 7.80 Å and b = 6.05 Å contains two molecules (marked by arrows). The second molecule is placed in the middle of the unit cell. f) Same as (e) but with the second molecule of the unit cell rotated by 60 • about the surface normal and shifted to avoid clashes. The , , and angles are the same for the both molecules of the unit cell. g) Octamer model (selected from the SAM model (e)) used as a starting point for geometry optimization (of two molecules in the middle marked by arrows). h) Optimized octamer structure (only two middle molecules were optimized).
Furthermore, several works have addressed exciton formation in covalently linked multiazobenzenes. [38][39][40][41][42][43][44][45][46] Here, we expand the previous studies by treating large 1D and 2D clusters of azobenzenes (up to 18 molecules) by means of first-principles quantum chemical methods such as timedependent long-range corrected density functional theory (TD-lc-DFT) and ab initio configuration interaction singles (CIS). Nowadays, thanks to the progress in computational resources, these methods can be applied to large molecules/clusters of molecules, and, importantly, they do not suffer from the so-called charge transfer (CT) problem of TD-DFT. [47,48] Furthermore, we analyze the nature of the exciton states using the TDM analysis, which allows us to judge on exciton delocalization pattern and importance of charge transfer contributions for the studied azobenzene aggregates. We also perform BSE/GW calculations for a selected chain model to reveal the effect of infinite, periodic arrangement on azobenzene optical absorption.

Models
The models considered in the present work are shown in Figure 1. The geometry of the azobenzene monomer (Figure 1a) was optimized in the electronic ground state using the B3LYP functional [31,32] and the def2-TZVP basis set. [49] Using this geometry, we constructed linear chains of face-to-face stacked azobenzenes (H-type chains, Figure 1b). The number of molecules in the chain, N, was varied from 1 to 10, and the distance between neighboring molecules, d, from 3.0 to 5.5 Å. Then we considered 2D, SAM-like arrangements. First, a model geometry of upright standing molecules (perpendicular to an imaginary surface) was studied (the top view is shown in Figure 1c). In this model, the rectangular unit cell with lattice parameters a = 3.5 Å and b = 6.0 Å contains one molecule. Octadecamer (6 molecules in a-direction × 3 molecules in b-direction) and dodecamer (4 molecules in a-direction × 3 molecules in b-direction) clusters of this type were studied. The dodecamer model was used in our earlier work on nonadiabatic dynamics in SAMs. [21,22] Further, in earlier experimental work, the molecular orientation of the azobenzene units within SAMs has been probed experimentally by means of near-edge X-ray absorption fine structure (NEXAFS) spectroscopy, and the corresponding tilt and twist angles describing molecular orientation with respect to a surface of a substrate have been reported. [9,50] Specifically, in ref. [9] the following angles defining the molecular orientation were obtained: = 76 • , = 30 • (the tilt angle), = 61 • (the twist angle); see Figure 1d for the definition of the angles. It was also reported that there are two molecules per unit cell in the investigated azobenzene-containing SAMs. [9,50] Using lattice parameters of ref. [9], a = 7.80 Å and b = 6.05 Å, we first constructed the SAM model shown in Figure 1e. In this model, the second molecule of the rectangular unit cell is obtained by translation of the first molecule to the middle of the unit cell, that is, by vector 0.5a + 0.5b. The a-direction, in turn, is chosen such that the projection of the normal to the molecular plane (green arrow in Figure 1d) onto the surface (ab-plane) builds angle of 30 • with the a-direction. Overall there are 18 molecules in this model and we will refer to it as 18a. A similar model but containing five molecules has been studied computationally by Utecht et al. [29] We note that the mutual orientation of the molecules within the unit cell is not directly extractable from the experiment. [9] A possible herringbone arrangement of molecules in the SAMs has been proposed based on analysis of scanning tunneling microscope (STM) and atomic force microscope (AFM) images. [9,50] Therefore, we also considered a herringbone structure. As described in ref. [9], the two angles and (see Figure 1d) fix the orientation of the molecule with respect to the surface but the second molecule of the unit cell can still be rotated about the surface normal while preserving , , and angles. (In passing we note that the three angles are related as cos = sin × cos .) The second molecule was obtained by first rotating the first molecule by 60 • about the surface normal, and then translating it by vector 0.735a + 0.735b. This vector was chosen such to avoid clashes with neighboring molecules. The herringbone model constructed this way is shown in Figure 1f. It consists of 18 molecules in total and we will refer to it as 18b. Finally, we performed a partial optimization starting from an octamer structure shown in Figure 1g (which was cut from model 18a) and optimizing two molecules in the middle (the "second" molecules of the unit cells) while keeping the six perimeter molecules (the "first" molecules of the unit cells) frozen. The optimization was performed at the B3LYP+D3(BJ) [51] /def2-SV(P) [49] level, that is, with a dispersion corrected DFT method. The optimized structure is shown in Figure 1h. It reveals a herringbone arrangement, but the optimized molecules change the orientation to the surface. Namely, ≈ 71 • , ≈ 19 • , ≈ 16 • for the lower optimized molecule, and ≈ 73 • , ≈ 18 • , ≈ 14 • for the upper optimized molecule. The geometries of the studied models are provided in the Supporting Information (in the geometries.xyz file).
Finally, we also considered a periodic chain model with d = 3.5 Å (like in Figure 1b, but an infinitely long chain) with a vacuum layer of 8.9 Å in a -and 15.7 Å in b-direction (see Figure S1, Supporting Information for details). For comparison we also calculated absorption of the "periodic" monomer, setting the intermolecular separation to 8 Å, similarly to ref. [34].
We note that all models studied herein do not include a surface, in contrast to SAMs studied experimentally. Often, to weaken the effect of a metallic surface, the azobenzene chromophores are attached to the surface using rather long molecular linkers, for example, alkanethiols. [52] On the other hand, there are examples of a more direct linkage to, for example, a silicon surface. [5,53] The surface may affect molecular states and lead to further spectral changes in addition to the excitonic effects. [20,33] In the present work, we address the effects arising solely from the intermolecular interactions of azobenzene chromophores.

Methods
The excitation energies and oscillator strengths for the models described above were computed using linear response timedependent (TD) methods (Hartree-Fock, DFT, and semiempirical) and the configuration interaction singles (CIS) scheme. [54,55] We used two global hybrid functionals, B3LYP [31,32] (20% of exact exchange) and BHandHLYP (50% of exact exchange; as im-plemented in the Gaussian package, for "half-and-half" theory see ref. [56]), and five range-separated hybrid functionals, CAM-B3LYP, [57] M11, [58] B97, [59] B97X, [59] and B97X-D. [60] We also used TD-HF, ab initio CIS, as well as TD and CIS methods in combination with AM1 [61] and PM7 [62] semiempirical methods. In addition, the long-range corrected density functional tight binding method (TD-lc-DFTB) [63,64] was tested. For calculations with non-semiempirical methods the def2-SV(P) basis set was used. [49] Absorption spectra were calculated using excitation energies and oscillator strengths as Here, I is intensity, E is excitation energy, E i and f i are the calculated excitation energy and oscillator strength, respectively, for the S 0 → S i transition, and is a broadening parameter ( = 0.05 eV in this work). Exciton (de)localization and charge transfer contributions were analyzed using the "fraction of transition density matrix" (FTDM) matrix, obtained for a given transition as [23] Here, P [AO] is the TDM in atomic orbital (AO) basis (computed with Multiwfn 3.8 [65] ) and S is the AO overlap matrix. Diagonal elements F XX quantify contributions of local excitations and off-diagonal elements F XY , Y ≠ X charge-transfer excitations (X and Y denote monomers of an aggregate). In what follows, we will express the FTDM elements in %, that is as F XY × 100%. The overall CT contribution for a given transition was quantified as the sum of all off-diagonal FTDM elements, ∑ X≠Y F XY .
All aforementioned calculations besides the TD-lc-DFTB ones were done using Gaussian 16. [66] The TD-lc-DFTB calculations were done using DFTBaby (release January 11, 2021). [64,67] We also performed cluster (molecular) and periodic BSE@G 0 W 0 @DFT calculations. The cluster BSE@G 0 W 0 @DFT calculations were done with the FHI-aims package. [68,69] First, a standard self-consistent DFT calculation with the PBE functional [70] was performed with a numeric atom-centered tier 2 basis to give the Kohn-Sham eigenvalues and the corresponding orbitals. A perturbative one shot G 0 W 0 calculation [36] was then applied on top of the PBE reference, to get the self energy and the corresponding quasiparticle energies. A two pole fitting for the analytical continuation [71] was applied with a 200 grid point Gauss-Legendre grid for the numerical integration of the integral over the imaginary frequency axis. The Tamm-Dancoff approximation (TDA) [72] was applied to decrease computational cost.
The optical properties in the periodic chain and periodic monomer case were calculated with the linearized augmented plane waves plus local orbitals (LAPW+LO) based code exciting. [73,74] First, azobenezene units were put in a periodic box of a dimension a = 20 Å, b = 20 Å, and c = 3.5 Å or c = 8 Å for the periodic chain and monomer case, respectively Figure 2. Spectra of the monomer. 50 excited states were requested in Gaussian calculations, and 30 excited states in the DFTBaby calculation. Vertical gray lines show the experimental, gas-phase peak position. [75] ( Figure S1, Supporting Information). The DFT ground state calculations were performed with the PBE functional by sampling the Brillouin zone with a homogeneous 1 × 1 × 4 Monkhorst-Pack k-point grid for the chain and a Γ point calculation for the monomer case. Muffin tin orbitals with radii of 0.8 bohr for H, 1.2 bohr for C, and 1.2 bohr for N have been used with |G + k| max of 3.5 and |G| of 14 including 100 empty states for the ground state properties. Quasiparticle (QP) energies were calculated in the G 0 W 0 approximation including 100 empty states in the confined space of ten occupied and ten unoccupied states.
On top of the QP energies BSE calculations with 200 empty states for the screened Coulomb potential have been performed with a slightly shifted k-grid. The BSE Hamiltonian included ten occupied and ten unoccupied bands. The eigenvalue problem was simplified through the TDA approximation.
An optical absorption stick spectrum is given by the imaginary part of the macroscopic dielectric tensor M : [34] Here, |t i | 2 are oscillator strengths of transitions with excitation energies E i , and Ω is the unit cell volume. The exciting code performs a Lorentzian broadening with a value of 0.2 eV in order to mimic the finite excitation lifetimes. Band contributions to the excitons were calculated from BSE eigenvectors as described in ref. [34] (see Equation (5) there).

Monomer
As a first step of our study we calculated absorption spectra of trans azobenzene using TD-DFT, TD-HF, and CIS with the def2-SV(P) basis set. A moderate (double-zeta) basis set was chosen to enable direct comparison with calculations for large aggregates (see below), for which better basis sets (e.g., triple-zeta) are impractical. In addition, we tested several semiempirical methods. The calculated absorption spectra of the monomer are shown in Figure 2. The lowest energy peak (for each method) corresponds to the S 0 → S 2 , * transition. The S 0 → S 1 , n * transition is optically forbidden. The calculated * absorption peak position ranges from 2.97 eV (TD-PM7) to 4.79 eV (CIS). The experimental, gas-phase value is about 4.12 eV. [75] The ab initio CIS yields the most blue-shifted band (4.79 eV), followed by TD-HF (4.53). Notably, the difference between TD-HF and CIS is quite large, 0.26 eV, which may point to a ground-state RHF→UHF (triplet) instability. [76] Indeed, the stability analysis confirmed the presence of the triplet instability on the HF/def2-SV(P) level, for the used monomer geometry (optimized with B3LYP/def2-TZVP). Further to the red follow the long-range corrected hybrids B97 (4.44 eV), M11 (4.39 eV), and B97X (4.35 eV), see spectra in shades of blue in Figure 2. Next to them (to the red) go the global hybrid BHandHLYP, with 50% of exact exchange (4.21 eV), long-range corrected functional B97X-D (4.19 eV), and range-separated functional CAM-B3LYP (4.17 eV), see spectra in shades of green. These three methods yield the * peak positions closest to the experimental value (4.12 eV). In this respect, we should note that we compare here only vertical excitation energies, neglecting vibrational fine structure. Further to the red follows the semiempirical TD-lc-DFTB method (3.90 eV) and global hybrid B3LYP, with 20% of exact exchange (3.85 eV). Then follow semiempirical methods AM1/CIS (3.50 eV), TD-AM1 (3.44 eV), PM7/CIS (3.12 eV), and the most red-shifted TD-PM7 (2.97 eV). We also note that the oscillator strength of the S 0 → S 2 transition is much smaller for the used semiempirical methods (f < 0.5) than for TD-DFT, TD-HF, and ab initio CIS (f > 0.75).
The dependence of spectra on basis set and geometry as well as TDM analysis for the S 1 and S 2 states are given in ref. [23]. The excited states of azobenzene have also been investigated with higher level wave function based methods (such as CC2, ADC(2), and CASPT2), [77][78][79] but these methods are too costly for large aggregates considered in what follows; therefore, we do not discuss them here.

H-Type Linear Chains
We first consider predictions of a simple exciton model based on PDA of the intermolecular coupling. For a linear chain model of N equidistantly distributed identical molecules, assuming periodic boundary conditions, the exciton states are |k⟩ = 1 √ N ∑ N n=1 e i 2 k N n |n⟩ with k = 0, 1, 2, ⋯, N − 1 and |n⟩ being the states with excitation on molecule n and all other molecules staying in their ground states. [29,80] The transition dipole moments from the aggregate ground state to the exciton states are ⃗ k=0 = √ N ⃗ mon and ⃗ k≠0 = 0, where ⃗ mon is the transition dipole moment of the monomer. [19,29] Further, for the H-type arrangement (e.g., cofacially stacked molecules with parallel transition dipole moments), the |k = 0⟩ state is the upper state of the exciton band [29,80] and, thus, the upper state is the only bright state in this model. The transition to the |k = 0⟩ state is expected to be blue-shifted with respect to the monomeric transition if the exciton coupling shift surpasses the van der Waals shift. [20] Now, without periodic boundary conditions, the exciton states can be calculated by diagonalizing the Frenkel-exciton Hamilto-nianĤ = ∑ n E mon |n⟩⟨n| + ∑ m ∑ n ≠ m V mn |m⟩⟨n| (see, e.g., refs. [81,82]). E mon is the monomer excitation energy and V mn are the exciton couplings. Further, in PDM, for an H-aggregate, the exciton coupling between molecules m and n is given by V mn = 2 mon R 3 mn (R mn is the distance between point transition dipoles of molecules m and n, which build an angle of 90 • with R mn ). [17,18] For the S 0 → S 2 transition of azobenzene, 2 mon = 8.13 (in atomic units) at the TD-B97X-D/def2-SV(P) level. In Figure 3 the spectra of the linear decamer (with d = 5.5 Å) calculated using PDM and TD-DFT (with the B97X-D functional) are shown. First of all, we find, that not only the upper state is bright, but also four other states have nonzero transition dipole moments, in contrast to the model with periodic boundary conditions. The bright states alternate with dark states. It is interesting to note that this alternation has been observed for other systems (brick-wall-type PTCDA aggregates) as well. [83] Second, the exciton band width (or exciton splitting, which we define as the difference of the energies of the highest-and the lowest-energy exciton states) is much larger for PDM than for TD-DFT, 780 meV versus 307 meV, respectively. This situation is also realized for shorter chains, see Table 1, where we list PDM, B97X-D, and CIS exciton band widths.
Here and in what follows we select B97X-D and CIS methods since: i) they do not suffer from spurious, low energy CT problem; ii) B97X-D performs well for the monomer (see Figure 2); and iii) ab initio CIS is an affordable wave function based method for excited states of large aggregates. Further, we note that it was also reported for other systems (acenes) that exciton couplings calculated with PDM are larger than those calculated using the TD-DFT supramolecular approach. [84] We note that for azobenzene the molecular length, calculated as the para-H-para-H distance, is about 11 Å, that is, two times larger than the considered intermolecular separation d = 5.5 Å. Therefore, it is expected that PDM may fail short in this case. In general, the total electronic coupling is composed of the Coulomb coupling (also referred to as the long-range coupling) and the shortrange coupling, which includes exchange and intermolecular orbital overlap effects. [85] The coupling calculated using PDM is an approximation for the Coulomb coupling. It has been reported that for short intermolecular distances the PDM Coulomb coupling is (much) larger than the "exact" Coulomb coupling (which, in turn, may be calculated from transition densities [86] or specially fitted transition charges [87] ). [88] This is one source of the PDM error in the calculation of excitation energies. The other source is the short-range coupling which is completely neglected in PDM. Overall, the total coupling may be smaller than the PDM Coulomb coupling. [89,90] Further, we analyzed characteristics of the brightest exciton state for linear chains of varying length and different intermolecular separations. The state label, excitation energy, and oscillator strength of the brightest state calculated with TD-B97X-D/def2-SV(P) and CIS/def2-SV(P) are shown in Figure 4. For relatively large separation distances (4-5.5 Å), the brightest state is excited state 2N, as expected in the molecular exciton model introduced above. [91] For smaller distances, deviations from this expectation are observed. In particular, the brightest state carries a label larger than 2N (for longer chains) starting at 3.7 Å for TD-B97X-D calculations and 3.5 Å for CIS calculations. The reason for this deviation is that: i) exciton states arising from the monomer states higher than S 2 may appear below the brightest exciton (which, in turn, arises from monomeric S 2 ); and ii) charge-resonance states [28] (or states with admixture of charge transfer contributions) stemming from S 2 (or S 1 ) may also have lower energies than the brightest exciton. In turn, (i) is caused by large dispersive red shift and, moreover, the corresponding states may exhibit some appreciable CT character. [23] As a simple example, the brighest exciton of the trimer 3.5 Å is S 8 at the TD-B97X-D level (instead of S 6 expected from Kasha's exciton model). Inspection of NTOs of the S 6 and S 7 states of the trimer reveals that these states correspond to the S 3∕4 states of the monomer (see the NTOs in Figures S2 and S3, Supporting Information). We also note that mixed transitions, involving n → * and → * excitations, may also take place (see the NTOs of the S 0 → S 10 transition in Figure S3, Supporting Information).
Further, in the middle row of Figure 4 we see that the excitation energy of the brightest transition monotonically rises with the number of molecules, with an exception of d = 3 Å. The excitation energy rises more strongly for shorter nearest neighbor distances d, as expected from PDM. However, quantitatively the first-principles excitation energies are smaller than those predicted with PDM (compare TD-DFT/CIS and PDM curves for d = 5.5 Å in Figure 4, middle row). For d = 3 Å the excitation energy curves go down below the shown energy window since there is no bright state within the first calculated 50 states, starting at N = 7 for TD-B97X-D and at N = 9 for CIS, as can be seen in the lower row of Figure 4 where the oscillator strength of the brightest transition is shown. There we also see that for fixed N the oscillator strength is larger for larger d. It is opposite to expectation from the simple exciton model. In the latter, the oscillator strength is given (in atomic units) by f = 2 3 E| ∑ n c n | 2 2 mon , where E is the excitation energy of an electronic transition for the aggregate, c n are components of the eigenvector ofĤ, and mon is the transition dipole moment of the monomer. Since c n does not depend on separation, [92] the only difference in the oscillator strength comes from E, which, in turn, is larger for smaller d. We also note that PDM oscillator strengths are larger than those calculated by first-principles methods (Figure 4, lower row). There are two sources of the discrepancy: i) the aforementioned error in the excitation energy E contributes to the error in the oscillator strength; and ii) the transition dipole moment, calculated from the monomer transition dipole and eigenvectors of the PDM exciton-Hamiltonian, is also not exact. For example, for the linear decamer with d = 5.5 Å, for the brightest transition, we obtain the following transition dipole strengths ( 2 , in atomic units): 73.5 using PDM; 55.4 at the TD-B97X-D level.
Further, for the decamer d = 3.5 Å, we calculated absorption spectra using the methods applied to the monomer. The spectra are shown in Figure 5. The position of the * band follows the same trend (with respect to the used methods) as for the monomer. We note that the B3LYP curve (the black line y ≈ 0 in the upper panel) does not show the intense * peak, since the requested 50 lowest excited states turn out to be not enough to reach it (the TD-B3LYP 50th excited state is located at 3.13 eV, which is below the S 0 → S 2 transition of the monomer at 3.85 eV, compare with Figure 2).
The brightest state labels and energies calculated with different methods as well as the blue shifts with respect to the monomeric S 0 → S 2 transition are listed in Table 2. The semiempirical methods and TD-HF predict the brightest transition to be the S 0 → S 20 transition, agreeing with the simple exciton model. Ab initio CIS and TD-DFT, on contrary, predict the brightest state label to be larger than 20, thus placing some additional states (which are not expected in the simple exciton model) below it. In particular, CIS yields S 0 → S 25 as the most intense transition, whereas the TD-DFT brightest state labels vary from 28 to 44, depending on a functional. In particular, for the used TD-DFT methods, the more blue-shifted is the * absorption band the smaller is the brightest state label. For larger intermolecular distances, however, TD-DFT labels become smaller, for example, 20 for B97X-D, for the decamer with d = 5.5 Å. Further, we analyzed the charge transfer contribution to the calculated excited states for TD-B97X-D and CIS spectra, for decamers d = 3.5 Å and d = 5.5 Å. The sum of off-diagonal (CT) elements of FTDM matrices for the lowest 50 excited states of these decamers is shown in Figure 6. The brightest transitions are marked with stars there. It is seen that for d = 5.5 Å, the CT contribution is virtually zero for transitions below the brightest one. In fact, it is zero for all 50 calculated transitions at the CIS level, whereas for B97X-D the first 30 transitions are charge transfer free and the next 20 transitions are entirely of the CT type. For d = 3.5 Å, on contrary, the observed picture is more complex. First ten transitions (the n * transitions) show (very) minor CT, which is more pronounced for B97X-D than for CIS. For the higher energy transitions the CT contribution becomes larger, for example, 37% ( B97X-D) or 21% (CIS) for the S 0 → S 11 transition. And some of the transitions with even higher excitation energies are the CT-dominated transitions. In this respect, we note that the low-lying CT states are present in some organic solids, for example, in the pentacene crystall. [93] For the brightest transitions the CT contribution is small (≈8% for B97X-D and ≈2% for CIS).
The FTDM matrices for the brightest transitions of the decamers with increasing d are shown in Figure 7. First, the matrices are virtually diagonal, with some small nearest neighbor Adv. Theory Simul. 2023, 6, 2200907 Figure 5. Spectra of the linear decamer with the intermolecular separation of 3.5 Å. 50 excited states were requested in Gaussian calculations and 30 excited states in the DFTBaby calculation. Vertical gray lines show the experimental, gas-phase monomer peak position. [75] Table 2. Brightest state labels (Label dec ) for the linear decamer with 3.5 Å nearest-neighbor distance, and spectral shifts (E dec − E mon ) with respect to the monomer. Excitation energies of the brightest transitions of the decamer (E dec ) and the monomer (E mon ) are also provided. For the monomer, the brightest transition is S 0 → S 2 .

Method
Label CT contributions at smaller d. Next, we see that while the exciton is delocalized over the decamer, the contributions of local excitations of central monomers are larger than those of the terminal molecules-this is the effect of the finite size of the chain. Furthermore, we see that the diagonal values change with d (the heat map normalization in Figure 7 is done with respect to the maximum element among all shown matrices). We note that it is not expected from the simple PDM exciton theory, in which the wave function coefficients c n (see above) do not change with varying separation of chromophores. The diagonal values F XX as a function of monomer X are plotted in Figure 8, for TD-B97X-D and CIS calculations. The F XX curves become more concave  with increasing d. This effect is more pronounced at the TD-B97X-D level than at the CIS level. The overall CT contributions (shown with dotted lines in Figure 8) are higher for TD-B97X-D. Namely, at the TD-B97X-D level they are 8%, 4%, 1%, 0.3% and 0.05% for d = 3.5, 3.7, 4.0, 4.5, and 5.5 Å, respectively. At the CIS level, the corresponding total CT contributions are 1.6%, 1.1%, 0.6%, 0.2% and 0.02%.
Having examined the TD-DFT/CIS results for the chains consisting of up to ten molecules, we turn now to the BSE@G 0 W 0 @PBE calculations. We were particularly interested in the absorption spectrum of an infinite (periodic) chain. But at first we performed BSE@G 0 W 0 @PBE calculations for the monomer. Here, two approaches were used-i) BSE@G 0 W 0 @PBE calculations with numeric atom-centered orbitals as implemented in FHI-aims; and ii) BSE@G 0 W 0 @PBE calculations with linearized augmented plane waves as implemented in exciting. The monomer spectra calculated using these two approaches are shown in orange in Figure 9. The lowest energy absorption maxima are at 3.43 eV (FHI-aims) and 3.50 eV (exciting), so there is a remarkable agreement between the two approaches. We also note that these values are slightly smaller than the excitation energy calculated at the BSE@evGW@PBE level (evGW is an iterative GW method, in which the quasiparticle gap is converged within 0.1 eV), which is ≈3.7 eV. [37] We also recall that the experimental excitation energy of the gas-phase monomer is about 4.12 eV, as mentioned above. [75] The higher energy parts of the monomer spectra, however, show a larger discrepancy (compare orange curves in Figure 9), which we attribute to differences between the used approaches (molecular vs periodic). We next calculated the spectra of the dimer and the trimer, both with d = 3.5 Å, using FHI-aims (Figure 9, upper panel). The monomer-to-dimer blue shift is ≈0.21 eV and the monomer-totrimer blue shift is ≈0.29 eV. These shifts are larger in comparison to those calculated with TD-DFT (≈0.13 eV and ≈0.20 eV, respectively, at the TD-B97X-D level) and with CIS (≈0.16 eV and ≈0.23 eV, respectively). exciting calculations for the monomer and H-type chain 3.5 Å. Vertical gray line shows the experimental, gas-phase monomer peak position. [75] Furthermore, we computed the absorption spectrum of the infinite linear chain with d = 3.5 Å, using exciting (see Figure 9, lower panel). Here, we observe a very pronounced monomer-tochain blue shift, of ≈1.6 eV. This blue shift is much larger than the TD-DFT/CIS monomer-to-decamer blue shifts of about 0.3-0.4 eV (see Table 2). The strong chain peak at ≈5.1 eV is only ≈1.7 times intenser than the monomer peak at ≈3.5 eV. We note, however, that the spectrum is normalized to the unit cell in this case. We also note weak absorption (in the BSE@G 0 W 0 @PBE spectrum) at energies near and below the monomer peak, which, in principle, can also be used to excite the chain.
A deeper analysis shows that the weak transition at ≈3.7 eV in the periodic chain arises mainly from the VBM−1 to CBM excitations at Γ and K. (Here, VBM stands for the valence-band max-imum and CBM for the conduction-band minimum.) The same states contribute to the monomer exciton at 3.5 eV. For the intense exciton transition at ≈5.1 eV (for the chain), apart from the dominant VBM−1 to CBM contribution, there is also a sizable contribution from VBM−1 to CBM+1 at K. The contributions of the bands to these excitons are shown in Figure 10.

2D Arrangements
We consider first the dodecamer 2D model with arrangement of molecules shown in Figure 1c. The spectra of the dodecamer calculated with various methods are shown in Figure 11. Again, the order of the methods regarding the position of the most intense band is the same as for the monomer (compare with Figure 2), and this order agrees with that for the linear decamer (see Figure 5). The B3LYP "spectrum" again shows no bands (the black line y ≈ 0 in the upper panel) as in the case of the decamer, since 80 excited states requested in the present case are not enough to reveal the intense dodecamer * transition with B3LYP (the 80th state is located at 3.16 eV).
Comparing to the linear decamer, the monomer-to-aggregate blue shift (of the * transition) is slightly larger, ≈0.4 eV for the 2D dodecamer versus ≈0.3 eV for the 1D decamer at the TD-DFT levels used here (compare Tables 3 and 2). The dodecamer brightest state labels are higher than those for the decamer, for example, 51 versus 43 at the TD-B97X-D level. TD-HF and the used semiempirical methods (with exception of TD-AM1) again agree with the simple exciton model, predicting the brightest state to be 24 (see Table 3). TD-AM1 predicts it to be 36, but we should note that the TD-AM1 density of states is very large in the region of excitation energies close to the brightest state-TD-AM1 state 24 is at ≈3.66 eV and state 48 is at ≈3.70 eV.
Furthermore, we computed the spectra of a larger model, octadecamer shown in Figure 1c. Additional blue shift and absorbance enhancement of the most intense peak are observed (Figure 12). Specifically, in comparison to the dodecamer the * band is blue-shifted by ≈0.06 (0.07) eV and the maximum absorbance is enhanced in ≈1. 3 (1.3) times at the TD-B97X-D (CIS) level. We also note some weak absorbance near and below   the monomer S 0 → S 2 transition. The brightest state labels are 97 and 48 for TD-B97X-D and CIS calculations, respectively. The FTDM matrices for the transitions to these states are shown in Figure 13, left column. The dominant contributions are on the diagonal. For B97X-D, there are also some nonzero CT contribu- tions between -stacked nearest neighbors (see molecule numbering in Figure 13, right). The spatial distribution of local excitations (F XX elements) is shown in Figure 13, middle column. It is seen that the brightest exciton exhibits the largest contribution  from the central dimer, similarly to the case of the linear decamer (compare to Figure 7). Next, we calculated the absorption spectra of the SAM models shown in Figure 1e,f, which are constructed using the experimental information on lattice parameters and orientation of azobenzenes with respect to a surface, as reported in ref. [9]. The spectra calculated with TD-B97X-D and CIS are shown in Figure 14.
The most intense band is more blue-shifted for the 18a model ( Figure 1e) than for the 18b model (Figure 1f). Specifically, the 18b-to-18a shift is ≈0.10 eV for TD-B97X-D and ≈0.05 eV for CIS. Another remarkable difference is a more pronounced low-energy absorbance for the 18b model. In this respect we note that an appearance of the allowed, red-shifted transition is expected already for the unit cell dimer, since the two molecules of the unit cell are not parallel to each other. [18] Indeed, we find for the isolated dimer, corresponding to the unit cell of 18b, the S 0 → S 3 transition at 3.89 eV with oscillator strength f = 0.20, in addition to the S 0 → S 4 transition at 4.17 eV with f = 1.05 (at the TD-B97X-D level). For the isolated dimer from 18a, on the contrary, the S 0 → S 3 transition at 4.08 eV is dark (f = 0.00), and only the S 0 → S 4 transition at 4.27 eV is bright (f = 1.49), see Figure S4, Supporting Information. Interestingly, for the dimer from model 18b, the S 3 and S 4 states are somewhat localized on either monomer ( Figure S5, Supporting Information), in contrast to the simple exciton model, in which each monomer contributes 50% (c n = 1∕ √ 2) to the overall excitation. This demonstrates the importance of specific molecular orientation at short distances. We also note that the localization of excitons on a single monomer was reported for edge-to-face benzene dimers. [94] The brightest transitions for 18a are S 0 → S 44 (with B97X-D) and S 0 → S 36 (with CIS). For 18b, it is S 0 → S 36 with both methods. The FTDM matrices for the brightest transitions are shown in Figures 15 and 16 for 18a and 18b, respectively. For 18a (Figure 15), TD-B97X-D and CIS matrices look similar. The exciton is delocalized over the entire aggregate with the largest contribution from the two central molecules located on the spatial antidiagonal (Figure 15, middle column). There are also minor CT contributions between molecules in close contact, for example, 2 ↔ 3, 2 ↔ 7, etc. We also note that the TD-B97X-D FTDM matrix possesses some slight asymmetry, for example, F 9,9 = 6.4% ≠ 6.3% = F 10,10 . We relate it to the fact that S 44 is degenerate with S 43 at the TD-B97X-D level (both states are 4.3864 eV above the ground state; the corresponding oscillator strengths are 0.001 for S 43 and 3.903 for S 44 ). Previously we encountered such a situation for an azobenzene dimer. [23]   In the case of 18b, the difference between TD-B97X-D and CIS FTDM matrices is large (Figure 16). The CIS brightest exciton is more delocalized than the TD-B97X-D exciton. The largest diagonal element of the CIS FTDM matrix is F 10,10 = 8.4%, whereas for TD-B97X-D it is F 11,11 = 21.4%. We note, however, that the stick spectra obtained with the two methods are qualitatively different. While CIS predicts one strong transition (S 0 → S 36 at 5.03 eV with f = 8.11), TD-B97X-D predicts a series of closely lying transitions with smaller oscillator strengths (the most intese among them are S 0 → S 36 at 4.29 eV with f = 2.00 and S 0 → S 38 at 4.30 eV with f = 1.49), see Figure S6, Supporting Information.
Finally, we performed a partial geometry optimization for an octamer model selected from 18a (see Figure 1g,h), and calculated the spectra for the initial and optimized structures (Figure 17). The intense * band for the optimized structure is blue-shifted by ≈0.06 eV in comparison to the band for the initial structure. In both cases, the brightest transition is S 0 → S 16 .
As described in Section 2.1, the optimized structure adopts a herringbone arrangement (see Figure 1h). The corresponding FTDM matrices are shown in Figure 18. While both excitons are delocalized over the octamers, the delocalization pattern is different, with dominant contributions from molecules 1 and 8 in the case of the inital structure, and from molecules 3 and 7 in the case of the optimized structure. It is also seen that the CT elements are reduced for the optimized structure.

Conclusions
In this work, we studied exciton states of azobenzene aggregates by means of first-principles calculations. In particular, we applied time-dependent long-range corrected density functional theory (TD-lc-DFT) and ab initio configuration interaction singles (CIS) to large clusters including up to 18 azobenzene molecules. The character of the exciton states was analyzed computing "fraction of transition density matrix" (FTDM) matrices allowing one to quantify exciton delocalization pattern and charge transfer contributions to electronic transitions. The brightest excitons of the studied aggregates are dominated by local excitations of monomers composing the aggregates. Moreover, for the investigated model clusters, these excitons are delocalized over the en- Figure 18. FTDM matrices (left column) and spatial distribution of their diagonal elements (middle column) for the brightest exciton of the initial (upper row) and optimized (lower row) 2D octamer models. Numbering of molecules is shown on the right. Matrix elements larger than 0.1% are printed. Calculations are performed at the TD-B97X-D/def2-SV(P) level. tire cluster (with exception of the TD-B97X-D calculation for 18b). The charge transfer (CT) contributions to the brightest transitions are small, decrease with increasing nearest neighbor distance, and are larger for TD-lc-DFT (specifically, B97X-D) than for CIS. The CT states appear earlier in TD-lc-DFT spectra than in the CIS spectra. This leads to large differences in the brightest state labels at short intermolecular distances. Comparison of TD-lc-DFT/CIS results with point-dipole model (PDM) calculations shows that PDM overestimates excitation energies and oscillator strengths in the present case. Moreover, we found that treating a periodic chain at the BSE/GW level of theory leads to a larger blue shift in comparison to the calculations on cluster models. Finally, it should be mentioned that the (de)localization of the exciton states is affected by conformational disorder, as was demonstrated for an azobenzene tetramer in our recent work. [95] Ultimately, the photoinduced dynamics may lead to further changes in (de)localization and result in exciton transfer between monomers, as was reported very recently by Sangiogo Gil et al., based on surface hopping calculations employing an exciton model. [96]

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.