We study the excitonic coupling and homogeneous spectral line width of brick layer J-aggregate films. We begin by analysing the structural information revealed by the two-exciton states probed in two-dimensional spectra. Our first main result is that the relation between the excitonic couplings and the spectral shift in a two-dimensional structure is different (larger shift for the same nearest neighbour coupling) from that in a one-dimensional structure, which leads to an estimation of dipolar coupling in two-dimensional lattices. We next investigate the mechanisms of homogeneous broadening—population relaxation and pure dephasing—and evaluate their relative importance in linear and two-dimensional aggregates. Our second main result is that pure dephasing dominates the line width in two-dimensional systems up to a crossover temperature, which explains the linear temperature dependence of the homogeneous line width. This is directly related to the decreased density of states at the band edge when compared with linear aggregates, thus reducing the contribution of population relaxation to dephasing. Pump-probe experiments are suggested to directly measure the lifetime of the bright state and can therefore support the proposed model.

## I. INTRODUCTION

Organic molecules are promising candidates for the next generation of electronic devices and for solar energy conversion.^{1–6} Among these, assemblies of dye molecules in the form of J-aggregates have attracted attention for their special optical properties.^{7,8} These are understood from delocalization of the exciton formed upon the absorption of light over tens to hundreds of monomers.^{9} Because the dominant resonant transfer interactions between molecules are negative in a J-aggregate, the optically bright state is found at the bottom of the band, leading to a redshift of the absorption peak compared to a single chromophore. Other properties resulting from this exciton delocalization are superradiance and a hidden level structure at the band edge.^{10,11} The exciton delocalization in linear aggregates is reflected in the pump-probe spectrum^{12,13} and in the two-dimensional optical spectrum.^{14}

Most early studies focused on J-aggregates for which the optical properties can be explained with a model of a linear aggregate, which self-assemble in solution and are often studied at low temperature in a glass environment.^{15} Over the past years there has been an intense interest in tubular J-aggregates.^{16–22} It is also possible to manufacture two-dimensional thin film J-aggregates of chromophore molecules, which were found to exhibit a redshift in the absorption.^{23,24} Nonlinear optical experiments produced a two-dimensional spectrum similar to the spectrum of a linear aggregate, consisting of a single pair of positive and negative peaks.^{23}

In order to analyse these findings, a model of a truly two-dimensional aggregate must be used,^{25,26} which goes beyond weakly coupled linear aggregates.^{27} In general, the transfer interactions between molecules depend strongly on their relative orientation. This means that the absorption spectrum is sensitive to the details of the molecular arrangement. Therefore, modeling of the spectrum can help in determining the structure. This is particularly helpful in cases where the structure is not known from other measurements.^{23}

It is clear from experiment, as well as from consideration of the molecular structure, that the excitons in J-aggregates must couple to their environment. This coupling leads to scattering between exciton states and to pure dephasing. The total dephasing process can be studied by measuring the homogeneous line width in experiments such as photon echos, hole burning, or two-dimensional spectroscopy. In particular, the dependence of the homogeneous line width on temperature can be analysed in order to understand the exciton phonon coupling mechanism. Understanding the interaction with the environment is also important to assess how many molecules in the aggregate are entangled upon optical excitation. This entanglement, or delocalization, is the key quantity that causes the interesting properties of these systems. However, interactions with phonons limit the localization size, and determining this quantity as a function of temperature is an important goal. It has been shown that the energy dependent localization size can be extracted from two-dimensional spectra.^{28}

Studies of the homogeneous line width have been performed on J-aggregates for which linear chain models explain the optical properties. Fidder *et al.*^{29} measured the homogeneous width in PIC-Br for temperatures between 1.5 and 190 K. The temperature dependence is clearly nonlinear and was modeled by coupling of the excitons to three harmonic modes with frequencies of 9 cm^{−1}, 305 cm^{−1}, and 973 cm^{−1}. Hirschmann and Friedrich^{30} measured the homogeneous line width in PIC-I for temperatures from 0.35 K to 80 K. The temperature dependence can be fitted by a sum of two exponentials or be explained by a theory that predicts a power law dependence.^{31} The homogeneous line shape is well approximated by a Lorentzian for all temperatures.

In contrast to this work, experiments on thin films using two-dimensional optical spectroscopy have found a *linear* scaling of the homogeneous line width with temperature.^{23} This suggests that a different mechanism is responsible for the line width.

In this work, we use an excitonic model of two-dimensional brick layer J-aggregates and study the homogeneous line width as a function of temperature. Our theoretical model is presented in Section II and the results are presented in Section III. Specifically, we predict the excitonic coupling in molecular aggregates and thus correlate two-dimensional spectra with molecular arrangements. These calculations are the topic of Sections III A–III C. Then we consider population relaxation and the exciton lifetime in Sec. III D and calculate the homogeneous line width as a function of temperature in Sec. III E. Finally, we analyse experimental measurements in Sec. IV and conclude in Section V. In the Appendix, we consider alternative aggregate geometries.

## II. MODEL

The usual Frenkel exciton model of J-aggregates starts from a single bright optical transition on each molecule. The Hamiltonian includes a term which describes local excitation of a molecule with an excitation energy *ϵ _{n}* and a term for the coherent exciton motion from one molecule to the other, and is given in terms of the Pauli creation and annihilation operators

*c*

^{†}and

*c*by

In this Hamiltonian, the sums run over all molecules in the aggregate. If the molecules are far enough apart, the electrostatic interaction between them can be approximated by dipole-dipole coupling, which gives

Here, *μ _{n}* is the transition dipole vector of molecule

*n*, $r\u2192nm=r\u2192n\u2212r\u2192m$ is the relative position vector, $rnm=|r\u2192nm|$ is the distance, and $r\u02c6nm=r\u2192nm/rnm$.

*C*is a constant that scales the magnitude of the coupling and includes possible rescaling effects due to vibrations.

^{32}Here, we will use transition dipole coupling for all pairs of molecules. For nearest neighbours, a better understanding of the coupling can be obtained from quantum chemical calculations.

^{33,34}

Each molecule in the aggregate is influenced by a different local environment. This leads to static disorder in the site energies *ϵ _{n}*, which are different for each aggregate in the ensemble. This, in turn, leads to localization and a distribution of effective sizes of the exciton. Furthermore, the excitons in the aggregate interact with phonons in the surrounding material. Their dynamic effect is usually modeled as a reservoir of harmonic oscillators, which are described by the bath Hamiltonian

*H*. The interaction of these oscillators with the electronic excitations is assumed to be

_{B} where the effective bath coordinate is to be thought of as the sum of couplings to individual bath modes, which can be written as $Xn=\u2212\u2211\alpha gn\alpha x\alpha $. Here, *x*_{α} denote the coordinates of the bath modes, while *g*_{nα} are their coupling constants to the system. The linear dependence on the bath coordinate can be thought of as a lowest order expansion in the coordinate. The properties of the system bath interactions are determined by the correlation functions 〈*X _{n}*(

*t*)

*X*(0)〉. We will make the usual but not completely general assumption that the fluctuations on each site are uncorrelated and that their correlation function is the same on each site, 〈

_{m}*X*(

_{n}*t*)

*X*(0)〉 =

_{m}*δ*(

_{nm}L*t*).

The system Hamiltonian (for each realization of the static disorder) can be diagonalized to give the exciton states *ϕ _{k}*, which we choose to be real, and energies

*E*, such that $HS=\u2211kEkck\u2020ck$. The wave functions relate the exciton basis to the site basis by the equation $ck\u2020=\u2211n\varphi kncn\u2020$. In the same exciton basis, the system bath interaction can be written as $HSB=HSB(0)+HSB\u2032$, with the diagonal fluctuations

_{k}and the off-diagonal fluctuations that couple two different eigenstates

In modified Redfield theory,^{35,36} the diagonal fluctuations are treated exactly, while the off-diagonal fluctuations are included in second order perturbation theory. The zero order Hamiltonian $HS+HSB(0)+HB$ does not couple the exciton states. Therefore, the absorption spectrum for this Hamiltonian is simply the sum of contributions from each exciton state. In this case, where the system Hamiltonian commutes with the system bath interaction, the linear and nonlinear response functions can be calculated analytically with the cumulant expansion. For the linear absorption in the time domain, we find

where *μ _{k}* is the transition dipole from the ground state to exciton state

*k*. The spectrum

*A*

^{(0)}(

*ω*) is given as the Fourier transform of

*A*

^{(0)}(

*t*). The line shape function for each exciton state is given by

*g*(

_{k}*t*) =

*g*(

*t*)/

*N*, where $Nk=1/\u2211n\varphi kn4$ is the inverse participation ratio.

_{k}^{35}The line shape function for a single site is defined as

If the harmonic bath is interpreted in the continuum limit, the correlation function can be expressed in terms of the spectral density *J*(*ω*) (although, in principle, the spectral density can also contain delta functions which describe discrete modes). The quantum correlation function *L*(*t*) is given in terms of the spectral density by the expression

The homogeneous line width, which can be measured in photon echos or two-dimensional spectra can now be explained by two broadening mechanisms. First, there is the pure dephasing contribution contained in Eq. (6). Second, there will be a contribution from $HSB\u2032$. In second order perturbation theory with the Markov and secular approximations, this term will lead to dephasing given as the sum of population relaxation rates. This contribution can be termed dephasing from population relaxation. The perturbative treatment of $HSB\u2032$, which is known from experiment to be weak in certain linear J-aggregates,^{31} leads to scattering between exciton eigenstates. The scattering rate between eigenstate *q* and *k* is given by

where the sum runs over all molecules in the aggregate, *ω _{kq}* =

*E*−

_{k}*E*,

_{q}*J*is the spectral density, and $n(\omega kq)=n\u0304(\omega kq)$ for

*ω*> 0 and $n(\omega kq)=n\u0304(\u2212\omega kq)+1$ for

_{kq}*ω*< 0, with $n\u0304(\omega )=(exp(\omega /kT)\u22121)\u22121$ the Bose-Einstein distribution. The resulting dephasing rate of exciton state

_{kq}*k*is given by $\Gamma k=(1/2)\u2211q\u2260kWqk$. Then, the homogenous absorption line with both pure dephasing and dephasing from population relaxation is given by

Note that we neglect the radiative life time in this expression, which normally gives a negligible contribution to the linewidth. Also, in this work we do not include an ensemble of localization sizes caused by the presence of static disorder. Finally, it should be observed that the homogeneous line width is not equal to the sum of dephasing rates if multiple transitions overlap in the spectrum.

## III. RESULTS

### A. Molecular arrangement and resonant transfer interactions

We assume that the molecules are placed on a brick layer lattice^{25,26} with aspect ratio *A* and slip *s*, see Figure 1. The grid has *N _{x}* molecules in the x-direction and

*N*in the y-direction. These sizes should not be interpreted as the physical size of the aggregate, but as the number of molecules over which an exciton is delocalized. Different geometries can then be obtained by varying the slip, while we assume a constant aspect ratio of

_{y}*A*= 2, which is a reasonable number for molecules typically used to form thin films. Note that with this choice, a slip of 1.0 (half a unit cell) corresponds to a square lattice with dipoles oriented at 45° with respect to the lattice vectors. We note that, for the small aggregates considered in this paper, the choice of boundary conditions is important. We limit our study to the boundary conditions shown in Fig. 1. In the Appendix, we consider the value of

*A*= 3. Other arrangements, for example, herringbone structures with two molecules per unit cell, are outside the scope of this paper. For a constant prefactor

*C*the magnitude of the couplings will change with

*s*. This is shown in Figure 1. We define the nearest neighbour coupling in the x-direction as

*J*. Note that we include all long range couplings in our model as well, and that we do not use periodic boundary conditions.

In Figure 1 we observe that the values of *s* for which an aggregate with negative couplings in both directions is formed are quite limited. In most cases, the coupling in the vertical direction is positive. Because of this combination of negative and positive interactions, the system is not a perfect J-aggregate, in the sense that the bright state is not necessarily at the bottom of the band. We will see that this has observable consequences for the two-dimensional optical spectrum. Negative couplings in both the x and y directions are found around *s* = 1.0, which is the structure close to the one assumed in experimental work on brick layer PTCDA aggregates.^{24}

In Figure 2 we plot the density of states for a linear aggregate and several two-dimensional aggregates. For each geometry, we scale the parameter *J* (or, equivalently *C*), to obtain a shift of approximately 2500 cm^{−1} of the aggregate absorption peak with respect to the absorption peak of the monomer. This choice is made to stay close to the interpretation of experimental results, in which the spectral shift upon aggregation can be measured, but the structure of the aggregate (i.e., a two-dimensional bricklayer lattice or a collection of semi-one-dimensional chains) is not always *a priori* known. In particular, we are interested in comparing with experimental data on BIC aggregates studied in Ref. 23, where this shift was observed to be around 2500 cm^{−1}. In order to obtain this shift, we set *J* = − 1050 cm^{−1} for the linear aggregate, *J* = − 903 cm^{−1} for *s* = 0.5, *J* = − 472 cm^{−1} for *s* = 0.75, and *J* = − 416 cm^{−1} for *s* = 1.0. We note that, because of differences in the electrostatic environment, the shift does not necessarily directly reflect the excitonic coupling.

In the linear case, we observe a strong increase of the density at the band edge, reflecting the $1/E$ scaling in the case of an infinite chain.^{37} The two-dimensional aggregate with *s* = 0.5 mirrors this behaviour, but has additional states at lower energy. For larger values of *s*, however, the density of states decreases with decreasing energy, in line with the expected constant behaviour in the limit of an infinite sheet.^{37} We note that the difference in density of states between a linear and two-dimensional aggregate is quite dramatic. For example, for the energy gap between the two lowest-lying states, we find 9.2 cm^{−1} in a 25 × 25 lattice (*s* = 1.0, *J* = − 416 cm^{−1}), while it is only 0.45 cm^{−1} in a linear chain with 625 molecules (*J* = − 1050 cm^{−1}). (Note that we include long range interactions in the simulations leading to this number.)

### B. Two-dimensional spectra

We now turn our attention to the linear and two-dimensional optical spectra for this model system. For all models the simulation predicts a dominant bright peak in the linear absorption spectra, which, by construction lies about 2500 cm^{−1} below the monomer absorption peak. However, we find that *s* = 0.5 and *s* = 1.0 are clearly distinguishable when the two-dimensional correlation spectrum is considered.

Two-dimensional optical spectroscopy is a third-order nonlinear optical technique which correlates the evolution of the electronic state of the system during two time periods, called *t*_{1} and *t*_{3}.^{38} The signal is plotted as a function of the Fourier transforms of these two time periods, with frequencies labeled *ω*_{1} and *ω*_{3}. The technique can be used to separate homogeneous broadening, which shows up as the anti-diagonal width of peaks in the two-dimensional plot, from inhomogeneous broadening, which contributes to the diagonal width. It therefore provides a tool to measure the homogeneous line width. When one looks at the real value of a two-dimensional spectrum, both negative peaks, colored in blue, and positive peaks, colored in red, are present. Blue peaks arise from interactions where one excitation is created, while red peaks correspond to processes where, during *t*_{3}, coherences between one- and two-quantum states (in which two excitation quanta are present in the system) are present. Positive peaks are blue shifted with respect to negative peaks as a consequence of the Pauli exclusion principle. The vertical distance between positive and negative peaks can be used as a ruler from which the exciton localization size can be determined.

We calculated two-dimensional optical spectra using the sum over states method,^{39} assuming only homogeneous broadening. This simple method will give a good idea of the peak positions and relative intensities, but not of the details of the line shape. Note that in the calculation of the spectra, we have chosen to vary *J* in order to obtain similar peak positions for all values of *s* considered. The reason for this choice is that *J* is not known *a priori* in experiment, but the spectral shift with respect to the monomer can be measured. In our spectra, in the case *s* = 0.5, because of the presence of both positive and negative couplings in the system, the bright state is not at the bottom of the band. Because of the Pauli exclusion principle, two excitons cannot populate the same state.^{13,40} Induced absorption peaks will show up at lower *ω*_{3} than the bleaching and stimulated emission peak. This is most easily understood in the simplified case of a linear aggregate with nearest neighbor interactions only, for which the Hamiltonian can be diagonalized analytically using a Jordan-Wigner transformation.^{41} The two-exciton states are then given as anti-symmetric products of one-exciton states. Because the “first” exciton is not at the bottom of the band, the “second” exciton can go to a lower energy than the first one. This leads to the induced absorption peak. Although this picture of two independent excitons is not strictly valid when long-range interactions are taken into account, the result that an extra induced absorption peak appears at low *ω*_{3} is also found numerically in the full calculation (see Fig. 3).

Note that for a linear aggregate the induced absorption peak is weaker than the bleaching and stimulated emission peak, while the two have approximately equal amplitude for the 2D lattice with *s* = 1.0. Both of these systems have the state with largest oscillator strength at the bottom of the band. The presence of a large coupling in the y-direction makes the 2D lattice different from a product of one-dimensional aggregates, as considered in Arias *et al.*^{23} We will see that the difference is crucial for a correct determination of the couplings in the system. In this case the bright state is at the bottom of the band, and the two-dimensional spectrum shows a dominant pair of a positive and a negative peak (see Fig. 3). Because the occurrence of positive and negative peaks and their relative intensities in this spectrum are close to the observed spectrum for BIC aggregates,^{23} this finding lends support to the model with a slip around *s* = 0.75–1.0 for this system. Finally, we note that in the 2D lattice spectrum, a cross peak due to finite size effects is just visible to the right of the main peak.

### C. Estimating the nearest neighbour couplings

We now consider the difference between a linear (N^{2} × 1) and a brick layer (N × N) aggregate. In both cases, a dominant pair of positive and negative peaks is found in the 2D spectrum.^{23,42}

The relation between the nearest neighbor coupling and the peak shift of the aggregate compared to the monomer, which can be measured by finding the maximum in the linear spectra, is different. This simple method, which relies on the fact that excitonic coupling shifts the peak, is frequently used to determine the coupling from experimentally obtained spectra, even though it neglects the shift in the single molecule transition frequency due to the electrostatic effect of the different environment in both cases. The difference is easily explained from the fact that there are more neighbors and that there is therefore more coupling in the 2D aggregate. Quantitatively, the peak shift in a linear aggregate is 2.4 times the nearest neighbor coupling.^{15} In the 2D system with *s* = 1.0, we find that the shift is several times larger (4.9 *J* for a 6 × 6 bricklayer lattice). Thus, while estimating the couplings from the spectral shift, it is important to take the aggregate geometry into account. Here, we established the rule that can be used to estimate the coupling in a two-dimensional brick layer aggregate from the measured linear absorption spectrum for parameters *s* = 1, *A* = 2. For other parameters a similar rule can be established, which will, in general, be different from the rule derived from calculations on linear aggregates. We note that the magnitude of the couplings and, therefore, spectral shifts strongly depend on the aggregate under study. In particular, in Ref. 24 much smaller couplings were found than in Ref. 23. However, the relation between spectral shift and coupling does not depend on the absolute value of these numbers. We note that the practical rule described here is related to well established results derived from sum rules.^{43}

### D. Population relaxation

We next consider the life time of the bright state, which for both the linear aggregate and the brick layer lattice with *A* = 2 and *s* = 1.0 lies at the bottom of the band. This life time can be measured experimentally using pump-probe spectroscopy. For the well studied system of linear aggregates in a glass, the temperature dependence of the life time of this state is mostly determined by scattering to higher states with the absorption of phonons. The dephasing associated with population relaxation is found to dominate, leading to a *T*^{3.5} dependence of the homogeneous line width.^{31} This result is obtained from perturbation theory, where the scattering rate between two eigenstates is the product of three contributions: the Boltzmann factor of the phonons *n* evaluated at the energy gap between both states, the exciton-phonon spectral density evaluated at the energy gap, which is taken cubic for the glass environment, and the overlap of the wave functions of the two states. We refer to Ref. 31 for further details. In this model, it is clear that the density of states at the band edge is very important for the life time.

As was shown before, this density of states is very different in linear and 2D aggregates. As a result, the established theory for linear aggregates in a glass must be used with caution here. The Boltzmann factor is much larger and more strongly peaked in the linear chain, leading to larger scattering rates, and a smaller life time. Although the spectral density typically increases with energy, this is counteracted by the stronger exponential decay of the Boltzmann factor. We also expect different behavior of the life time as a function of temperature, because the argument used in Heijs *et al.* to arrive at the *T*^{3.5} dependence is valid only for *kT* > *E*_{2} − *E*_{1}. We therefore expect a different behaviour of the temperature dependence of the homogeneous line width, irrespective of the details of the interaction with the phonons.

As the numerical calculations presented in Fig. 4 show, the life time is up to more than an order of magnitude larger in the brick layer system than in the linear aggregate for the same system (delocalization) size. Life times were computed for the same cubic spectral density in all cases and, in contrast to Secs. III A–III C, for a nearest neighbour coupling strength of *J* = − 500 cm^{−1}, which is a typical value for the aggregates considered. All long range couplings were included as well. We note that the population relaxation rate as a function of temperature exhibits a power law dependence for linear aggregates, while a kink is observed for a 2D lattice. The kink can be understood as follows. In the argument leading to a power law dependence of the scattering rates as a function of temperature, the summation over discrete states is replaced by an integration.^{31} This replacement is valid if the temperature is large compared to the energy gap between the relevant exciton states, which are the lowest two states in the band. For the 10 × 10 brick layer lattice used in Fig. 4, the energy gap between these states is found to be 71 cm^{−1}, which corresponds to a temperature of roughly 100 K. Therefore, the power law dependence, which holds for high temperatures compared to the energy gap, breaks down and a kink is observed in the life time. We also note that the temperature dependence of the life time depends very strongly on *A* and *s*. As discussed in the Appendix, there is almost no temperature dependence for certain values, while for other values the life time varies over orders of magnitude.

Note that, to make a direct comparison possible, these calculations assume that the spectral density is the same for linear and two-dimensional aggregates. We will argue later, based on the analysis of the experimentally measured temperature dependence of the homogeneous line width, that there is a linear component in the spectral density for the two-dimensional aggregates. However, the suppression of relaxation rates in the two-dimensional aggregates based on the much smaller Boltzmann factor will occur irrespective of the spectral density.

If, indeed, a linear spectral density is more appropriate for two-dimensional aggregates, we can also obtain population relaxation rates for this case. The result is shown in Figure 5. Parameters were estimated based on the experimental data for a BIC aggregate, see Section IV. For this choice of the spectral density, the pure dephasing can be characterized by a rate if the temperature is low compared to the cut-off frequency of the bath.^{44} The life time is also plotted in the figure. We observe that pure dephasing dominates the homogeneous line width up to a certain temperature, which depends on the system size used in the simulations (i.e., the exciton localization size). For higher temperatures, population relaxation becomes more important. This finding can be used to determine an upper limit for the delocalization size.

### E. Dephasing in two-dimensional aggregates

Comparison with experiment shows that the predominance of dephasing from population relaxation, which is found for linear aggregates, does not explain experiment for the thin film. A perfectly linear relation between the homogeneous line width and the temperature was measured.^{23} Although this finding could be explained from population relaxation with a different *ω* dependence of the spectral density, this is not a plausible explanation of the experimental findings. A sub-Ohmic power^{23} on the order of *ω*^{0.5} in the spectral density would be needed. There is no clear microscopic mechanism that would lead to this behavior. Furthermore, it would be very surprising that the temperature dependence of the homogeneous line width is exactly linear. Therefore, it is more logical to assume that the mechanism leading to homogeneous broadening is different in linear and 2D aggregates.

The different mechanism can be understood by the difference in population relaxation rates found in Sec. III D. In linear aggregates, the density of states at the band edge is large, leading to strong scattering between exciton states and fast relaxation of the population in the bright state. This population relaxation leads to dephasing, which dominates the homogeneous line width. In contrast, in two-dimensional aggregates, the density of states at the band edge is orders of magnitude smaller. Population relaxation, which depends strongly on the difference in energy between eigenstates, is suppressed for the smooth spectral densities that we assume here. The accompanying dephasing is therefore also much weaker than in the case of linear aggregates. Therefore, pure dephasing becomes more important and can dominate the homogeneous line width. This is the main finding of this paper. This model predicts that the life time broadening is only a small part of the homogeneous line width. This prediction could be tested experimentally by measuring the life time of the bright state with pump-probe spectroscopy.

We will now see how the dominance of pure dephasing in two-dimensional aggregates leads to a linear scaling of the homogeneous line width with temperature. We assume that the spectral density is polynomial in frequency up to a cut-off frequency *ω _{C}*,

where *f _{J}*(

*x*) is a cut-off function. Common forms for this function are an exponential, a Lorentzian, or a Heaviside step function. While this spectral density quite generally describes the interaction with the bath, we assume here that intramolecular vibrations do not play a role.

We are now in a position to analyse the pure dephasing term *g*(*t*). In the fast modulation limit, which we assume here because, as we will see, it leads to a linear relation between the homogeneous line width and temperature, one can replace the correlation function by a delta function, *L*(*t*) = Γ^{pure}*δ*(*t*). The pure dephasing contribution to the spectrum is then given by a rate Γ^{pure}, because *g*(*t*) = Γ^{pure}*t*. By comparing with Eq. (8) one finds that the rate is related to the slope of the spectral density at zero frequency,^{45,46} $\Gamma pure=lim\omega \u21920J(\omega )/\beta \omega $. Therefore, there is no pure dephasing rate for a super-Ohmic spectral density (*α* > 1). Note that this statement is also valid outside the fast modulation limit because the linear term in the line shape function is also given by the slope of the spectral density at zero frequency in the more general case. There is still pure dephasing, but *g*(*t*) has no linear term. For an Ohmic spectral density (*α* = 1), which is the most commonly used form, because it corresponds to a linear density of states and a frequency independent exciton phonon coupling, we see that the pure dephasing rate is linearly proportional to temperature.^{47} Thus, for an Ohmic spectral density in the fast modulation limit, we expect a linear scaling of the homogeneous line width with temperature. We note that it would be desirable to measure the time scale of the bath directly to strengthen this argument.

To close this section, we briefly discuss Kubo stochastic line shape theory,^{48} which is valid in the high temperature limit, to estimate the temperature dependence of the pure dephasing in the slow dephasing limit. The correlation function in this case is given by

where the variance of the fluctuations can be expressed in terms of the reorganization energy *λ* and temperature by calculating the correlation function from the Drude-Lorentz spectral density $J(\omega )=2\lambda \omega C\omega /(\omega 2+\omega C2)$. One finds *σ*^{2} = 2*λ*/*β*. The lineshape function for this model is easily calculated to be

One sees that in the fast modulation limit *g*(*t*) = 2*λt*/*βω _{C}*, in agreement with the pure dephasing rate introduced earlier. In the slow modulation limit,

*g*(

*t*) =

*σ*

^{2}

*t*

^{2}/2, where

*σ*is related to the inhomogeneity in the energies of exciton states. The line shape in the frequency domain is a Gaussian with standard deviation proportional to

*σ*. Because the variance

*σ*

^{2}is linearly proportional to temperature, the line width scales as the square root of temperature in this regime. Note that this regime is less relevant for our discussion of the homogeneous line width, which is determined by fast fluctuations, while very slow fluctuation contributes only to the inhomogeneous line width. However, this analysis shows that if the fast modulation limit is not strictly applicable,

^{49}deviations from linear scaling of the line width with temperature are expected.

Before presenting a comparison to experiment, we briefly discuss the role of static disorder. As indicated by Eq. (13), at low temperature, the line-shape and its T-dependence will be dominated by inhomogeneous broadening resulting from static disorder. In fact, a published calculation [see Fig. 7 of Ref. 50] clearly shows the transition from inhomogeneous broadening to homogenous broadening in the disordered chain system. Interestingly, this transition corresponds to optimal diffusion along the chain, suggesting optimization when dynamics and static disorders balance. Another point is that the relative contribution of pure dephasing decreases with the size of localization, as shown in Fig. 5. It is known that the Anderson localization size scales with the disorder and this universal scaling depends critically on the dimensionality. This scaling and its implication on diffusion were reported in a recent study^{27} and will be further explored in two-dimensional J-aggregates.

## IV. EXTRACTING PARAMETERS FROM EXPERIMENT

By comparing the measured 2D spectra in Ref. 23 with our calculated spectra, we observe that a slip around 1.0 is the most plausible structure for these systems. A slip around 0.5 and smaller is ruled out, because an extra induced absorption peak appears below the bleaching and stimulated emission peak, which is not observed in experiment.

As explained in Section III C, the estimate of the resonant transfer interactions in two-dimensional aggregates should be done with a different formula than in linear aggregates. To obtain an estimate of the coupling *J*, one should divide the excitonic peak shift by a factor of 4.9. This factor depends on *s* and *A*. By doing this, we find that the nearest neighbor coupling in the aggregate studied by Arias *et al.* is around −75 meV for BIC (instead of −153 meV based on a linear model) and −104 meV for U3 (instead of −212 meV). In PTCDA aggregates considered by Mueller *et al.*, the shift was found to be 400 cm^{−1}. This would lead to an interaction of −82 cm^{−1}. This value corresponds to −10 meV, and is therefore significantly smaller than for the other two aggregates.

Under the assumption that the homogeneous line width is completely determined by pure dephasing, we can extract system parameters from the experimental data by following the discussion in Sec. III. The slope of the homogeneous line width versus temperature, Γ = *ST*, is given by *S* = 2*λk _{B}*/

*ω*, where

_{C}N_{k}*N*is the localization size of the exciton. Thus, we find that 2

_{k}*λ*/

*ω*= 0.268

_{C}*N*for BIC. To proceed further, it would be desirable to measure

_{k}*ω*independently, for example, from the time dependent Stokes shift.

_{C}From this estimate, by assuming a value for *N _{k}*, we can derive the parameters for the model by Heijs

*et al.*and calculate the population relaxation times as shown in Figure 5. We used the same large value of

*ω*= 10 J, which means that we are consistently in the Markovian regime. We conclude that in order to obtain a linear relationship between homogeneous life time and temperature up to 250 K, the exciton should be localized on a segment smaller than 9 molecules (3 × 3). For larger temperatures, deviations from linear behaviour are expected, as can be seen in the figure. This is consistent with the dynamic localization size estimated from experiments on BIC and U3,

_{C}^{23}and slightly smaller than values found for PTCDA.

^{24}Note that in our modified Redfield approach dynamic localization is not included by construction and that we therefore regard

*N*as a parameter in the calculations.

_{k}## V. CONCLUSION

In conclusion, we have used the standard Frenkel exciton model to study the excitonic properties in brick layer thin film J-aggregates. We have introduced a novel theory to explain the experimentally measured linear temperature dependence of the pure dephasing rate.

We have found that the exciton couplings determined from peak shifts in the linear spectrum depend strongly on the geometry of the system. In particular, the couplings in two-dimensional aggregates are a few times smaller than estimates based on a linear aggregate model.

To explain the linear scaling of the homogeneous line width with temperature, we propose pure dephasing as the main homogeneous broadening mechanism. Because the energy gaps between exciton states at the bottom of the band are much larger in two-dimensional than in linear aggregates, population relaxation is suppressed. This leads to a smaller contribution to the line width from population relaxation in two-dimensional aggregates, as well as to longer exciton life times. From the experimental data, we can extract the product of the reorganization energy and the typical bath reorganization time scale. Pump-probe experiments, which can measure the lifetime of excited states, are suggested to confirm whether dephasing due to population relaxation is indeed less important than pure dephasing. It will also be valuable to determine the time scale of the reorganization of the phonon environment experimentally.

A much simpler model could in principle explain the exact linear temperature dependence. Stochastic line shape theory for a single two-level system coupled to a bath in the fast modulation limit predicts a Lorentzian homogeneous line with a width Δ^{2}*τ*, where Δ is the standard deviation of the fluctuations and *τ* their correlation time. Because the variance of the fluctuations depends linearly on temperature, this would explain the observed temperature dependence irrespective of the form of the spectral density.

A microscopic model that could make the exciton state behave as an effective two-level system is the presence of correlated fluctuations. If the site energy fluctuations are not independent, as in the model of Heijs *et al.*,^{31} but correlated over a distance comparable to the exciton localization size, they will not lead to scattering between eigenstates, and therefore no contribution to the line width from scattering in the exciton band occurs. We believe that the localization size is large enough to make the model of correlated fluctuations not plausible.

Our calculations are limited by the assumption in modified Redfield theory, where off-diagonal fluctuations in the exciton basis are treated perturbatively and therefore dynamic localization is not included. They could be improved by taking both the diagonal and off-diagonal system bath coupling terms into account non-perturbatively. This could be achieved with the cumulant expansion method,^{51,52} stochastic path integrals,^{53} or hierarchy of equation of motion simulations.^{54} These methods can be used to study the effect of dynamic localization,^{55} but are difficult to apply given the large reorganization energy and low temperature compared to the excitonic coupling in the system. They could also be used to calculate the localization size, which is treated as a parameter in our work. We remark that modified Redfield theory becomes more accurate with an increase in the energy gap between the two states involved, which means that the theory is more accurate for systems of modest localization size, as we consider here. For some aggregates, intermolecular vibrations are important,^{20,24} which could be considered in future work. Experimental measurement of the Stokes shift may help to improve our understanding of the role of vibrations. Finally, it would be desirable to get a better estimate of the resonant transfer interactions from quantum chemical simulations and to include radiative decay to the ground state in our calculations. On the experimental side, collecting 2D spectra as a function of the waiting time to elucidate energy transfer mechanisms^{56} would be a possible future extension.

The results presented here contribute to the understanding of the photophysics of two-dimensional J-aggregate thin films.

## Acknowledgments

Work by K.A.N. and A.G.D. was supported as part of the Center for Excitonics, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award No. DE-SC0001088. J.C. was supported by the NSF (Grant No. CHE-1112825). A.G.D. was furthermore supported by a Marie Curie International Incoming Fellowship within the 7th European Community Framework Programme (Grant No. 627864). H.G.D acknowledges financial support by the Joachim-Herz-Stiftung, Hamburg within the PIER Fellowship program.

### APPENDIX: MORE GEOMETRIES

In this Appendix, we consider a wider range of geometries than were presented in the main text, in particular, a brick layer lattice with a larger aspect ratio *A* = 3.

In Figure 6 we plot the resonant transfer interaction for an aspect ratio of 3. Similar as in the case of *A* = 2, we find that there is only a limited range of values of the slip *s* for which interactions are negative in the horizontal as well as the vertical direction. In Figure 7 we plot the density of states for *A* = 3. We observe a clear difference in the density of states at the band edge between linear aggregates and brick layer aggregates. We note that for a slip of *s* = 0.5, the state with maximum oscillator strength is blue shifted with respect to the monomer.

We calculated the life time of the bright state for a 6 × 6 brick layer lattice for different values of *A* and *s*, using the same procedure as outlined in the main text. For *A* = 2, with *s* = 0.25 and *s* = 0.5, the temperature dependence varies dramatically less than for the value of *s* = 1.0 used in the main text (data not shown). In particular, for *s* = 0.25 the lifetime is constant at the very small value of 1.5 ⋅ 10^{−7} ps over the entire range of temperatures considered (20–350 K). For *s* = 0.5, the life time decreases from 8 ⋅ 10^{−4} ps at 20 K to 2 ⋅ 10^{−4} ps at 350 K, still a much smaller variation than was observed in Fig. 4 for *s* = 1.0. We attribute the difference to the presence of states which are lower in energy than the bright state. The system can relax to these states by spontaneous emission of a phonon, which is independent of temperature.

For *A* = 3 and *s* = 0.5, we also find a constant lifetime over the range of temperature from 20 K to 350 K. For *s* = 1.0 and *s* = 1.5, we find a very strong dependence on temperature, as can be seen from the data plotted in Fig. 8. We conclude that measuring the temperature dependence of the life time of the bright state can be used as a tool to probe the molecular arrangement in brick layer aggregates.