In this paper, we present dyadic adaptive HOPS (DadHOPS), a new method for calculating linear absorption spectra for large molecular aggregates. This method combines the adaptive HOPS (adHOPS) framework, which uses locality to improve computational scaling, with the dyadic HOPS method previously developed to calculate linear and nonlinear spectroscopic signals. To construct a local representation of dyadic HOPS, we introduce an initial state decomposition that reconstructs the linear absorption spectra from a sum over locally excited initial conditions. We demonstrate the sum over initial conditions can be efficiently Monte Carlo sampled and that the corresponding calculations achieve size-invariant [i.e., $O(1)$] scaling for sufficiently large aggregates while trivially incorporating static disorder in the Hamiltonian. We present calculations on the photosystem I core complex to explore the behavior of the initial state decomposition in complex molecular aggregates as well as proof-of-concept DadHOPS calculations on an artificial molecular aggregate inspired by perylene bis-imide to demonstrate the size-invariance of the method.

## I. INTRODUCTION

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.^{6} 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,^{38} 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^{50} and nonlinear^{51} 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.^{5} 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],^{52} 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.

## II. THEORETICAL BACKGROUND

### A. Hamiltonian

*N*chromophores with an open quantum system Hamiltonian given by

*n*⟩) with vertical excitation energy (

*E*

_{n}), and an electronic coupling between pigments

*V*

_{nm}. The excited state of each pigment is linearly coupled,

*J*

_{n}(

*ω*) =

*∑*

_{q}|

*κ*

_{nq}|

^{2}

*δ*(

*ω*−

*ω*

_{nq}) and the inverse temperature $\beta =1kBT$. We decompose the bath correlation function into a sum of exponentials indexed by

*j*

_{n},

### B. Light–matter interaction

**) and the transition dipole operator of the individual pigments (**

*ɛ*

*μ*_{n}). The action of the collective dipole moment operator on the ground state of the aggregate is to excite the superposition state,

### C. Hierarchy of pure states (HOPS)

**z*** is a stochastic trajectory with components associated with individual thermal environments

*z*

_{n,t}defined by $E[zn,t]=0$, $E[zn,tzn,s]=0$, and $E[zn,tzn,s*]=\alpha n(t\u2212s)$. The time-evolution of the system wave functions can be calculated using the nonlinear HOPS equation

*γ*

_{n}), and

*k*

_{max}using the triangular truncation condition, where we only include auxiliary wave functions satisfying the condition ${\psi (k\u20d7)\u2208A:\u2211iki\u2264kmax}$, though other static filtering approaches have been explored previously.

^{53}

### D. Dyadic HOPS

*et al.*

^{54}and the subsequent work of Chen

*et al.*,

^{50}the dipole autocorrelation function that defines the linear absorption spectrum can be written as

*v*

_{η}(

*t*;

**z***)⟩ is the initial state

*η*is the coefficient associated with each of the four pure state initial conditions (

*η*∈ {±1, ±

*i*}). This correlation function reduces (accounting for cancellation between

*η*terms) to

## III. DYADIC ADAPTIVE HOPS

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.^{52} 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\u2295St)$ 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 $\u2202\u0303\Phi \u0303$ is less than a user-specified threshold *δ* (i.e., $||\u2202\Phi \u2212\u2202\u0303\Phi \u0303||<\delta $). 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 $\delta 2=\delta A2+\delta 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\u2295St)$ 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.

### A. Sparsity in electronic coupling

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*^{3} 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.

### B. Spectral density

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.

### C. Normalization

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 $\Delta A/S(t)=\delta A/S\u22c5||\psi (0\u20d7)(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.

## IV. MONTE CARLO SAMPLING LOCAL TRAJECTORIES

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.

*N*

_{pig}molecules with parallel dipole moments. The electronic coupling is assumed to be nearest neighbor (

*V*= −100 cm

^{−1}) and we describe the bath correlation function as

*λ*= 35 cm

^{−1}, Γ = 50 cm

^{−1},

*ν*= 10 cm

^{−1}, Γ

_{±}= Γ ±

*iν*, and

*T*= 295 K.

### A. Initial state decomposition

**d**) given by

*n*), i.e. $\sigma n=(|n\u232a\u2329g|+|g\u232a\u2329n|)$ and $An=\mu n\u22c5\epsilon $, Eq. (25) reduces to

#### 1. Example

*N*

_{pig}= 4) homogeneous chain (

*E*

_{n}= 0). In this case, the total correlation function can be decomposed into a sum of four single-site initial conditions [

*C*

_{n}(

*t*)]

*μ*

_{0}=

*μ*_{n}·

**. The symmetry of the Hamiltonian gives rise to only two unique initial conditions**

*ɛ*### B. Monte Carlo sampling dyadic HOPS

**z***) for each initial condition is independent of the sum over the

*N*

_{D}initial conditions (

**d**). If we assume unbiased sampling with equally sized clusters, then the number of noise trajectories per initial condition is

*N*

_{d,ens}=

*N*

_{ens}/

*N*

_{D}and the correlation function can be calculated as

*N*

_{ens}pairs (

**d**,

**z***) to calculate the full correlation function.

#### 1. Example

*n*,

**z***) space. We find that the statistical convergence of the total correlation function is equivalent to that of the edge (green circles) and inner (blue squares) spectral components calculated independently using Eq. (27).

One advantage of Monte Carlo sampling is the trivial incorporation of an additional sampling over disorder in the system Hamiltonian $(H\u0302s)$. 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^{−1}.

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\u22c5$ 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.

### C. Monte Carlo sampling dyadic adaptive HOPS

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^{4} trajectories [Fig. 4(a), gray horizontal line]. In Fig. 4(b), the edge spectral contribution calculated with converged DadHOPS parameters (thin line, *δ*_{A} = 10^{−3}) 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^{6} elements. Figure 4(c) plots the total absorption spectrum calculated with DadHOPS (thick line, *δ*_{A} = 10^{−3}, *δ*_{S} = 0) using the full initial condition $(|\psi 0\u232a=112\u2211n=112|n\u232a)$. 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^{−3}]. 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.

## V. APPLICATION

### A. Photosystem I (PSI)

^{59}The complicated spatial arrangement of chlorophyll makes this an ideal test system for exploring the behavior of the initial state decomposition with respect to different definitions of clusters. Here, we use a previously reported Hamiltonian

^{60}where chlorophyll excitation energies (

*E*

_{n}) and excitonic couplings (

*V*

_{n,m}) were calculated using time-dependent density functional theory (TDDFT) evaluated by Gaussian16

^{61}with the CAM-B3LYP density functional

^{62}and the 6-31G* basis set;

^{63}this approach is known to provide quality descriptions of both vibrational and vibronic processes apparent from the high-resolution absorption and emission spectra of chlorophyll-type molecules.

^{64}The Hamiltonian is provided in the supplementary material. The system–bath coupling is described by a Drude–Lorentz spectral density given by

*α*

_{n}(

*t*) is continuous at time

*t*= 0. To agree with previous HEOM calculations,

^{58}we introduce a Gaussian distributed static disorder on the chlorophyll site energies with a standard deviation (SD) of 150 cm

^{−1}, we approximate the true chlorophyll Qy transition dipole moment [Fig. 5(b)] by the vector defined by the position of the NA and NC nitrogen atoms (IUPAC notation for chlorophyll

*a*),

^{65}and we set the spectral density parameters to be

*λ*= 35 cm

^{−1}, Γ = 50 cm

^{−1},

*T*= 300 K, and Γ

_{mark}= 500 cm

^{−1}.

**),**

*ɛ***to the $Ad$ defined in Eq. (22) to indicate the polarization dependence. Figure 5(c) shows that, as expected, the HOPS calculations (green line) reproduce previous HEOM calculations from Ref. 58 (gray line).**

*ɛ*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.

### B. Perylene bis-imide (PBI)

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^{−1}) and 28 intramolecular vibrational modes with appreciable Huang–Rhys factors.^{70} 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),^{70} multilayer MCTDH (ML-MCTDH),^{71} and time-dependent density matrix renormalization group theory (TD-DMRG).^{72} 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}

Mode . | g (cm^{−2})
. | γ (cm^{−1})
. |
---|---|---|

1 | 1.2 × 10^{5} | 50 + 170i |

2 | 1.6 × 10^{6} | 100 + 1550i |

Mode . | g (cm^{−2})
. | γ (cm^{−1})
. |
---|---|---|

1 | 1.2 × 10^{5} | 50 + 170i |

2 | 1.6 × 10^{6} | 100 + 1550i |

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^{−1} to account for the group of strongly coupled vibrations ranging from 1370 to 1630 cm^{−1} previously reported from TDDFT calculations.^{70} 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.

^{−1}(blue line) compared with the SD = 300 cm

^{−1}(green line) we use for the aggregate spectrum. We can further improve the agreement between the simulated and experimental monomer spectra by introducing additional moderate frequency modes (data not shown), but this additional complexity was not required for the reproduction of the experimental aggregate spectrum. We smooth the aggregate spectra reported here by weighting the time-domain dipole autocorrelation function with a cosine appodization window,

^{75}

^{,}

*t*

_{max}). The use of this window function suppresses noise in the calculated spectra that arises from the combination of zero-padding

^{75}and the incomplete cancellation of the correlation function at long times due to finite sampling of the trajectory ensembles.

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.^{76} 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.

## VI. CONCLUSIONS

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.^{77}

## SUPPLEMENTARY MATERIAL

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

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

The data that support the findings of this study are openly available in Zenodo.^{78}

### APPENDIX A: DYADIC HOPS WITH GENERALIZED INITIAL STATE DECOMPOSITION

*A*

_{a}) defined to ensure that

*v*

_{η,a}⟩ is time-evolved according to the nonlinear HOPS equation [Eq. (11)]. As written, the time-evolution of the pure state |

*v*

_{η,a}⟩ is performed using a system Hamiltonian $(H\u0302S)$ that contains both the ground and excited electronic states. However, since the system Hamiltonian does not couple the electronic ground and excited states, the time-evolution of the two components of |

*v*

_{η,a}⟩ can be decoupled, giving

*ψ*

_{a}(

*t*;

**z***)⟩ is propagated using the nonlinear HOPS equation in which $H\u0302S$ includes only the first excitation manifold and the expectation value of the system–bath coupling operator in Eq. (13) is redefined as

*η*, the cancellation in the numerator due to summation over

*η*∈ {±1, ±

*i*} gives

#### 1. Cluster decomposition

*ψ*

_{d}(

*t*;

**z***)⟩ is the state $|\psi d\u232a=\sigma \u0302d|g\u232a$ time-evolved according to the dyadic HOPS equation.

#### 2. Pigment decomposition

*n*(

*t*;

**z***)⟩ is the single-pigment excitation state |

*n*⟩ time-evolved according to the dyadic HOPS equation.

### APPENDIX B: DEFINING STRONGLY INTERACTING CLUSTERS

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.

**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.**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*).**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*).**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.

## REFERENCES

*Photosynthetic Excitons*

*Charge and Energy Transfer Dynamics in Molecular Systems: A Theoretical Introduction*

*Principles of Nonlinear Optical Spectroscopy*

*Exciton Dynamics in Molecular Crystals and Aggregates*

Note 1, the most recent version of MesoHOPS can be found here: https://github.com/MesoscienceLab/mesohops.