Simulating optical linear absorption for mesoscale molecular aggregates: An adaptive hierarchy of pure states approach

Molecular aggregates, where the constituents are bound together by non-covalent interactions, are ubiquitous, appearing in contexts ranging from self-assembly in solution to photosynthetic pigment–protein complexes.^1–5 The long-range Coulomb interaction between individual components of an aggregate can lead to excited electronic states where the excitation is coherently delocalized over many molecules. These delocalized electronic states, which depend on the molecular arrangement and coupling to nuclear degrees of freedom, influence the remarkable optical and excitation transfer properties that make molecular aggregates appealing for technological applications, including solar energy harvesting, sensors, and organic light emitting diodes (OLEDs).

Optical spectroscopy provides essential insight into excited-state processes of molecular aggregates. The interpretation of spectroscopic observables in the presence of both structural heterogeneity and broad homogeneous line shapes typically requires detailed numerical simulations. However, established quantum dynamics methods can struggle with the combination of strong electron-vibrational coupling and spatially extended structures characteristic of molecular aggregates.⁶ As a result, the development of efficient algorithms for simulating time-resolved optical spectroscopy measurements of extended molecular aggregates remains a persistent challenge.

From a molecular perspective, simulating optical spectroscopy requires solving the time-evolution of nuclear wave packets on electronic potential energy surfaces. A variety of methods have been developed to solve this time-evolution within a formally exact framework, such as multiconfiguration time-dependent Hartree (MCTDH),^7–9 multilayer MCTDH (ML-MCTDH),^10–13 multiconfiguration Ehrenfest,^14–16 and ab initio multiple spawning.^17,18 However, the simultaneous description of electronic, intramolecular vibrational, and environmental degrees of freedom on equal footing introduces intractable computational scaling for large molecular aggregates: a problem colloquially known as the curse of dimensionality.

Open quantum system approaches provide a powerful coarse-graining where the relevant electronic degrees of freedom are time-evolved in a reduced density matrix while coupled to an effective thermal environment. Most commonly, an open quantum system Hamiltonian incorporates a linear coupling to an infinite collection of harmonic oscillators parameterized to mimic the influence of both the intramolecular vibrations and the surrounding (condensed phase) environment within a linear response approximation. In this framework, several formally exact methods are available for modeling molecular excitons, including hierarchical equations of motion (HEOM),^19,20 quasi-adiabatic path integrals,^21,22 time-dependent density matrix renormalization group theory (TD-DMRG),^23,24 time-evolving density operator with orthogonal polynomials,^25,26 and the multi-D1 Davydov ansatz.^27,28 While these methods are capable of scaling to much larger aggregates than are usually achievable with wave function propagation techniques, they still suffer from effectively exponential scaling with the number of electronic states (i.e., molecules) included in the calculations. Recent developments in reduced scaling techniques—including modular path integrals^29,30 and tensor contraction^31–35—suggest new paths forward may continue to extend the size of molecular aggregates that are tractable with these methods.

Recently, there has been a growing interest in solving open quantum systems using stochastic wave functions. In these approaches, the reduced density matrix is unraveled into an ensemble of wave functions that can be time-evolved independently. A variety of formally exact stochastic Schrödinger equation (SSE) methods^36,37 have been developed, including the hierarchy of forward–backward SSEs,³⁸ the stochastic Liouville–von Neumann (SLN) equation,^39,40 quantum jumps,^37,41 and non-Markovian quantum state diffusion (NMQSD).^42,43 In addition, there are a wide variety of stochastic equations that provide an approximate solution to the excited-state dynamics, including the well-known Haken–Strobl model.^44–46 All of these methods share the advantage of reduced memory scaling arising from propagating wave functions in place of a density matrix but introduce the need for many individual realizations to produce an ensemble of wave functions.

The hierarchy of pure states (HOPS) is a numerically convenient formulation of the NMQSD equations^47–49 that has been recently extended into the dyadic HOPS equation to simulate both linear⁵⁰ and nonlinear⁵¹ optical spectroscopy. Dyadic HOPS propagates the bra- and ket-side of the reduced density matrix separately according to the HOPS equation and provides a clear connection to the nonlinear response function formalism.⁵ Dyadic HOPS is limited to small molecular aggregates by the same poor computational scaling that plagues HOPS and other formally exact methods. The recent development of adaptive HOPS (adHOPS), which leverages the locality of physical wave functions to achieve computational costs that stop increasing with system size for sufficiently large aggregates [i.e., $O(1)$ scaling],⁵² raises the possibility of dramatically expanding the reach of dyadic HOPS within an adaptive framework. In this paper, we demonstrate dyadic adaptive HOPS (DadHOPS) calculations of linear absorption for large molecular aggregates. We present an initial state decomposition for the dipole autocorrelation function (which characterizes linear absorption) that provides a local construction of the spectroscopic observable. By combining DadHOPS with an initial state decomposition, we are able to perform size-invariant scaling simulations of linear absorption spectra for mesoscale molecular aggregates and easily incorporate static disorder in the system Hamiltonian.

This paper is organized as follows: In Sec. II, we provide a theoretical background for the readers, beginning with discussion of the Hamiltonian for an open quantum system and followed by a brief introduction for the HOPS and dyadic HOPS methods. In Sec. III, we describe an algorithm to connect adaptivity with dyadic HOPS method and provide conditions for relative error bounds that are necessary for DadHOPS. In Sec. IV, we apply the initial state decomposition method to decompose the dipole autocorrelation function into a sum over local correlation functions. We discuss an effective method for Monte Carlo sampling over stochastic noise trajectories and initial states with both dyadic HOPS and DadHOPS. Using 4-site and 12-site chains as examples, we demonstrate the applicability of the DadHOPS method. In Sec. V, we apply the initial state decomposition to photosystem I and use DadHOPS to simulate the absorption spectrum of a realistic model for artificial molecular aggregates inspired by Perylene bis-imide (PBI). Finally, we conclude in Sec. VI with a summary and a brief outlook. In  Appendix A, we provide a detailed derivation of the dipole autocorrelation function decomposition into a sum of local correlation functions based on a generalized initial state decomposition.

Here, we present an algorithm for constructing an adaptive basis during the propagation of a dyadic HOPS trajectory to substantially reduce computational requirements when simulating large systems. The algorithm presented here builds on previous work implementing adaptive HOPS for a density matrix calculation.⁵² Because HOPS trajectories are localized by the presence of their thermal environments (i.e., bath Hamiltonians), the time-evolution of their dynamics can be captured by a local basis that moves with the trajectory. To ensure the adaptive calculations retain the formally exact structure of HOPS, the local basis must guarantee a controllable error bound on the calculation.

Previously, Ref. 52 has demonstrated that it is possible to construct an adaptive HOPS algorithm where computational expense does not scale with system size for sufficiently large aggregates [i.e., $O(1)$ scaling]. Figure 1 presents a sketch of this algorithm, where the local basis is constructed as a direct sum $(At⊕St)$ of the set of auxiliary wave functions (⁠$At$⁠, called the “auxiliary basis”) and the set of molecular states (⁠$St$⁠, called the “state basis”). At each time point, the algorithm builds a new auxiliary and state basis that ensures the difference between the full derivative of all wave functions (represented by ∂Φ) and the derivative constructed in the reduced basis $∂̃Φ̃$ is less than a user-specified threshold δ (i.e., $||∂Φ−∂̃Φ̃||<δ$⁠). For convenience, we split the user-specified derivative error bound (δ) into two components, viz., the auxiliary derivative error bound (δ_A) and the state derivative error bound (δ_S), such that $δ2=δA2+δS2$⁠. Splitting the error bound allows us to independently control the precision of the calculation when constructing the auxiliary and state bases.

For adaptive HOPS, both memory requirements and computational expense scale with size of the current basis $(At⊕St)$ rather than the full basis. The decreased memory requirements arise directly from the structure of the basis itself. Computational expense, however, also depends on the algorithm used to construct the adaptive basis. Figure S3 from Ref. 52 shows that the computational expense of adHOPS as a function of size of a linear aggregate does not increase between calculations containing 16 and 1000 molecules. In other words, the adaptive HOPS calculations show size-invariant [i.e., $O(1)$] scaling of computational expense for sufficiently large aggregates.

All calculations reported here are run with MesoHOPS v1.2.1.^55,56 Below, we provide a few caveats and commentaries on how the current calculations compare to the previous description of the adaptive algorithm in Ref. 52.

First, to achieve size-invariant scaling, the adaptive HOPS algorithm requires that the system Hamiltonian be sparse. Physically, the electronic coupling between neutral molecules decays as 1/R³ where R is the distance between molecules, leading to a natural length scale beyond which coupling can be neglected. All adaptive HOPS calculations presented here are performed with system Hamiltonians that have nearest neighbor coupling and hence a natural sparsity.

Second, the original derivation of adaptive HOPS assumed the vibrational environment described by the spectral density contained only over-damped oscillators (e.g., Drude–Lorentz) for which, by construction, the amplitude of any associated auxiliary wave function saturates near unity. The same guarantee does not hold for under-damped oscillators (such as we use in Sec. V B), where the amplitude of the associated auxiliaries can become very large. As a result, in the calculations presented below, when we include under-damped vibrations, the adaptive HOPS scheme provides an inefficient truncation for the associated auxiliary wave functions because their large amplitude makes them disproportionately unlikely to be truncated. The result is a formally exact calculation that is presumably more expensive than an optimized adaptive scheme would require. Future work will focus on developing alternative adaptive approaches optimized for a more general spectral density.

Finally, the original derivation of the adaptive algorithm used a normalized nonlinear HOPS equation, but the dyadic HOPS construction uses a non-normalized (though still nonlinear) HOPS. As a result, dyadic adaptive HOPS (DadHOPS) requires an extension of our previous adaptive algorithm for an equation of motion that does not guarantee a normalized physical wave function. We do this by using the same adaptive algorithm while introducing new error bounds that are scaled relative to the magnitude of the physical wave function $ΔA/S(t)=δA/S⋅||ψ(0⃗)(t;z*)||$⁠. The result is that we account for the relative importance of different basis elements compared to the magnitude of the physical wave function. For example, as the magnitude of the physical wave function decays in time, the corresponding error bound becomes more stringent to ensure the relative accuracy is maintained.

Here, we present a method for combining the DadHOPS algorithm with a local representation of the dipole autocorrelation function to enable efficient simulation of absorption spectra for mesoscale molecular aggregates. For the DadHOPS algorithm to be efficient, the spatial extent of delocalization must be substantially smaller than the full size of the molecular aggregate. Simultaneously, when implemented directly, the dyadic HOPS expression for the dipole autocorrelation function has an initial condition that can be arbitrarily delocalized. In the following, we first demonstrate that the dipole autocorrelation function can be decomposed into contributions arising from a set of local initial conditions. Second, we demonstrate that the total correlation function can be reproduced by a simultaneous Monte Carlo sampling of noise trajectories $(zt*)$ and local initial conditions. Finally, we combine the initial state decomposition with the DadHOPS algorithm to demonstrate a local construction of a linear absorption spectrum that can achieve size-invariant [i.e., $O(1)$] scaling.

One advantage of Monte Carlo sampling is the trivial incorporation of an additional sampling over disorder in the system Hamiltonian $(Ĥs)$⁠. Figure 3(a) shows the statistical convergence of the total spectrum calculated with single-site initial conditions in the presence of Gaussian distributed disorder on site energies (black open triangles) is almost unchanged compared to the homogeneous case (black filled triangles). For the disordered calculations, site energies form a Gaussian distribution with mean value of zero and a standard deviation (SD) of 100 cm⁻¹.

Finally, we find that the statistical error of the Monte Carlo sample for a given number of trajectories is inversely proportional to the size of the clusters (N_c = N_pig/N_D) used for the initial state decomposition. In Fig. 3(b), we compare the normalized statistical error (⁠$Nc⋅$ Error) as a function of the number of independent trajectories (N_traj) when using different initial conditions: single site (triangles, N_c = 1), pairs (circles, N_c = 2), and all four sites (squares, N_c = 4). The equivalence of these three normalized results shows that increasing the size of the clusters decreases the number of trajectories required to reach a given statistical error. Thus, when using an equation of motion where the delocalization extent of the initial state does not change the computational expense, the initial state decomposition provides no advantage.

The local correlation functions arising from the initial state decomposition can be efficiently simulated using DadHOPS. Figure 4(a) shows how, for a 4-site linear chain, the error of the edge component decreases with the derivative error bound δ_A. Here, we consider DadHOPS calculations converged when the adaptive error is lower than the statistical error for the associated ensemble of 10⁴ trajectories [Fig. 4(a), gray horizontal line]. In Fig. 4(b), the edge spectral contribution calculated with converged DadHOPS parameters (thin line, δ_A = 10⁻³) reproduces the corresponding dyadic HOPS calculation (thick transparent line).

While dyadic HOPS is limited to relatively small system sizes, combining the local construction of the dipole correlation functions arising from the initial state decomposition with DadHOPS allows us to calculate linear absorption for much larger aggregates. As a first demonstration of scaling, we consider a 12-site linear chain where the basis for the full HOPS calculation has more than 10⁶ elements. Figure 4(c) plots the total absorption spectrum calculated with DadHOPS (thick line, δ_A = 10⁻³, δ_S = 0) using the full initial condition $(|ψ0〉=112∑n=112|n〉)$⁠. We can reproduce this spectrum by Monte Carlo sampling over single-site initial conditions and incorporating adaptivity in both the auxiliary and state basis [Fig. 4(c), thin line, δ_A = δ_S = 10⁻³]. Combining adaptivity in both the auxiliary and state basis with localized initial states can help achieve size-invariant scaling for sufficiently large aggregates. Figure 4(d) plots the number of state basis elements (bottom) and the number of auxiliary wave functions (top) required for both the full dyadic HOPS (thick gray lines) and the DadHOPS algorithm (thin green lines) as a function of the number of sites in the linear chain.

Here, we have demonstrated that the advantage of localized wave functions introduced by the initial state decomposition is realized by combining it with the DadHOPS equation of motion.

To what extent does the specific choice of clusters matter for a given cluster size? In the case of a linear chain, the choice of clusters to be sets of adjacent pigments may appear obvious, but there is no equivalent choice for the heterogeneously coupled pigments in PSI. Figure 5(d) compares the statistical convergence of the linear absorption spectrum for two different definitions of clusters: The first randomly assigns four pigments to each cluster (black solid circles), while the second constructs clusters of four strongly interacting pigments (green open circles) using an algorithm described in  Appendix B. We find that there is equivalent convergence either randomly assigning pigments to clusters or using clusters defined by strong coupling. Combined with the results presented in Fig. 3(b) showing square-root-scaling of error with cluster size, our calculations suggest that using a cluster initial condition only acts to increase the effective number of combined (n, z*) pairs that are being sampled.

Perylene bis-imide (PBI) is a family of pigments that can form both J- and H-type aggregates in solution.^66–69 J-aggregates formed by linear aggregates of PBI-1 have been the object of particular study and show a strong vibronic progression in their linear absorption spectra. Previous time-dependent density functional theory (TDDFT) calculations of PBI-1 have characterized electronic coupling between adjacent monomers (V ≈ −500 cm⁻¹) and 28 intramolecular vibrational modes with appreciable Huang–Rhys factors.⁷⁰ This vibrational structure has served as a starting point for different calculations of linear absorption spectra for PBI-1 aggregates using techniques such as multiconfiguration time-dependent Hartree (MCTDH),⁷⁰ multilayer MCTDH (ML-MCTDH),⁷¹ and time-dependent density matrix renormalization group theory (TD-DMRG).⁷² Calculations of the exciton dynamics in aggregates ranging in size from 2 to 25 monomer units have also been reported with varying degrees of complexity in their description of the molecular vibrations.^69,73,74

We have constructed a minimal bath correlation function composed of two modes to describe the vibrational environment of PBI. We include a vibration with a central frequency near 1500 cm⁻¹ to account for the group of strongly coupled vibrations ranging from 1370 to 1630 cm⁻¹ previously reported from TDDFT calculations.⁷⁰ We also include a low-frequency vibration to account for the broad dissipative environment. The final parameters (Table I) were selected to agree with the available experimental data.

The advantage of the adaptive dyadic HOPS formalism is the ability to calculate even very large molecular aggregates efficiently. The size of a molecular aggregate is often important for capturing the influence of exciton delocalization on both the 0-0 peak position and the relative magnitude of vibronic side-bands.⁷⁶  Figure 6(c) compares the monomer (black line), dimer (green line), trimer (blue line), and heptamer (cyan line) line shapes with that calculated for a 1000-site linear chain (gray line). The magnitude of the vibronic side-band decreases rapidly from monomer to trimer and then more slowly with increasing chain length. A similar trend is seen for the red-shift of the 0-0 peak. Figure 6(d) quantifies this effect by plotting the central position of the 0-0 peak (ΔE, top) and the ratio of the intensity of the 0–1 peak to the intensity of the 0-0 peak (I_0,1/I_0,0, bottom) as a function of the number of pigments in the aggregate. In both cases, the gray lines represent the asymptotic limit of a 1000-site linear chain.

Finally, let us consider the relationship between aggregate size, locality, and the onset of size-invariant scaling in the DadHOPS equation of motion. When the adaptive basis size stops increasing with the length of the aggregate, the extent of exciton delocalization is sufficiently small that most trajectories will never sample the edge. This implies that, for aggregates sufficiently large that the basis is size-invariant, the spectroscopic signatures should also be size-invariant. Figure 7 shows that by N_pig = 10, the average size of the adaptive site basis begins to plateau, before reaching a size-invariant value of 13 for chains of 100 PBI molecules. This is consistent with the convergence of the spectral features close to the asymptotic values at N_pig = 7 [Fig. 6(d)] and the expectation that complete size-invariance will occur when only a relatively small set of trajectories sample pigments on the edge of the chain. We suggest, then, that DadHOPS should be thought of as introducing a dynamic separation of length scales where the basis size required to capture the spectral features is solved on the fly rather than being asserted a priori. As a result, DadHOPS is advantageous for simulating realistic molecular aggregates, particularly in the presence of structural disorder, where separations of length scales can be obscure.

Here, we have presented a new algorithm, combining the dyadic adaptive HOPS (DadHOPS) equations with an initial state decomposition that is capable of calculating optical linear absorption spectra for mesoscale molecular aggregates. Our approach introduces a dynamic separation of length scales by adaptively constructing a reduced basis to describe the time-dependence of local contributions to the dipole autocorrelation function. The adaptive basis construction is capable of achieving size-invariant scaling [i.e., $O(1)$] with the number of molecules for sufficiently large aggregates. We have applied the initial state decomposition to the 96 chlorophyll photosystem I core complex and found that the specific choice of pigment clusters does not affect the statistical convergence of the calculations. We simulated a 1000-molecule J-aggregate of perylene bis-imide (PBI) with DadHOPS and characterized how the absorption spectra changes with aggregate size. Because locality cannot reduce the computational expense of describing individual pigments, DadHOPS remains limited to calculations where each thermal reservoir has a small number of exponential modes. In the future, it may be possible to describe a more complex thermal environment by incorporating techniques such as tensor contraction.⁷⁷ 

The Hamiltonian for the monomer of the photosystem I complex is available in the supplementary material.

The authors thank Jacob K. Lynd for editing and reviewing the manuscript and Harsh Dayal for the help with visual design. A.E. acknowledges support from the DFG via a Heisenberg fellowship (Grant No. EI 872/10–1). T.G., E.J.T., and D.I.G.B.R. acknowledge support from the Robert A. Welch Foundation (Grant No. N-2026-20200401) and start-up funds from Southern Methodist University. D.I.G.B.R. acknowledges the support received through US National Science Foundation CAREER Award under Grant No. CHE-2145358.

The authors have no conflicts to disclose.

Tarun Gera: Formal analysis (equal); Investigation (lead); Validation (equal); Visualization (lead); Writing – original draft (equal); Writing – review & editing (equal). Lipeng Chen: Conceptualization (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Alexander Eisfeld: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Jeffrey R. Reimers: Investigation (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Elliot J. Taffet: Investigation (equal); Methodology (equal); Writing – original draft (supporting); Writing – review & editing (supporting). Doran I. G. B. Raccah: Conceptualization (lead); Funding acquisition (equal); Investigation (equal); Methodology (equal); Software (lead); Supervision (lead); Validation (equal); Writing – original draft (lead); Writing – review & editing (equal).

The data that support the findings of this study are openly available in Zenodo.⁷⁸ 

To explore the influence of cluster definitions, we constructed clusters between four strongly coupled pigments in photosystem I (PSI). We use an iterative, three-step algorithm designed to ensure that strongly interacting pigments are preferentially included inside the same cluster.

1. Step 1: Locate the largest coupling element (V_i,j) among the pigments not yet assigned to a cluster. Pigments (i, j) form the nucleus of a new cluster.

2. Step 2: Locate the largest coupling element to either pigment in the new cluster to another pigment not yet assigned to a cluster (V_m,j or V_i,m). Add this pigment to the new cluster: (i, j, m).

3. Step 3: Locate the largest coupling element involving a pigment in the new cluster with another pigment (n) not yet assigned to a cluster. Add this pigment to the new cluster: (i, j, m, n).

4. End Condition: If any pigment remains unassigned, return to step 1 and define a new cluster.

This algorithm does not guarantee that all strong couplings are contained within a cluster, but it does ensure that clusters nucleate around strong coupling elements.