We consider theoretically near-field absorption spectra of molecular aggregates stemming from a scattering scanning near-field optical microscopy type setup. Our focus is on the dependence on the direction and polarization of the incoming electromagnetic radiation, which induces a Hertz dipole with a specific orientation at the tip-apex. Within a simple description, which is based on the eigenstates of the aggregate, absorption spectra are calculated for the near field created by this dipole. We find that the spatial patterns of the spectra have a strong dependence on the orientation of this tip-dipole, which can be understood by considering three basic functions that only depend on the arrangement of the aggregate and the molecule tip distance, but not on the orientation of the tip-dipole. This allows direct access to spatial dependence of the aggregate eigenstates. For the important cases of one- and two-dimensional systems with parallel molecules, we discuss these spectra in detail. The simple numerically efficient approach is validated by a more detailed description where the incoming radiation and the interaction between the tip and molecules are explicitly taken into account.

## I. INTRODUCTION

In recent years, there has been a growing interest in self-assembled molecular aggregates on dielectric surfaces.^{1–6} The surface degrees of freedom barely couple to molecular excitations, and fluorescence is only weakly quenched.^{1,7,8} Strong interactions between the transition dipoles of the molecules lead to delocalized excitonic eigenstates where an electronic excitation is *coherently* shared by many molecules.^{9–11} These eigenstates depend sensitively on the arrangement of the molecules. There have been many different arrangements observed.^{12–15} Common arrangements are one- and two-dimensional arrays where all molecules are parallel to each other.^{11,12,16,17}

Much effort has been devoted to obtain information about these delocalized eigenstates because they are crucial to understand the optical and transport properties of the molecular aggregate. It has been shown theoretically that by using optical fields that are *spatially inhomogeneous* on the scale of a few nanometers, one can not only access essentially all eigenstates, but also reconstruct the corresponding wavefunctions by moving the near field spatially and recording several spectra with the near field centered at different positions.^{18,19} Similar ideas to circumvent far-field selection rules and obtain information on the eigenstates have been discussed for various systems.^{20–23} There are many possibilities to realize such near fields experimentally and to detect the optical response.^{24–34}

In our previous investigations on extracting the aggregate wavefunctions from near-field spectra,^{18,19} a situation has been considered where the electromagnetic field that acts on the molecular aggregate stems from a source that is described by a point dipole. Such light sources can be approximately realized using metallic tips and considering the response caused by incoming electromagnetic far-field radiation. In Refs. 18 and 19, the field distribution was modeled as the field from a point dipole and the direction of the dipole was chosen perpendicular to the surface.

In the present work, we investigate the dependence of the spatially resolved spectra on the direction of the dipole that creates the near field. We see that for different orientations, the resulting spectra can be qualitatively quite different, which can be traced back to the underlying collective wavefunctions. The basic setup that we consider is sketched in Fig. 1. A metal tip (which we model by a small spherical particle) is located at some distance (∼1 nm) above the molecules that are arranged on a plane. In an experimental setup, near-field radiation could, for example, be created by using a localized surface plasmon^{28} or by a scanning near-field optical microscopy (SNOM) setup^{24,35} where a tip is irradiated with an external electromagnetic far-field (as indicated in Fig. 1). In this type of setup, also the molecules are irradiated by the external field.^{24,35} As we show below, a simple model where only the near field from a dipole induced in the tip is considered gives in many cases good results. We model this near field as stemming from a dipole whose direction $d\u20d7$ depends on the direction $k\u20d7inc$ and polarization $E\u20d7inc$ of the incoming external radiation. We see that the resulting spectra can be decomposed into three basic contributions that are weighted by simple functions that depend on the orientation of $d\u20d7$. These three contributions contain mutually complementary information about the underlying wavefunction with information that goes beyond that of our previous studies where only a single function was accessible and used for the wavefunction reconstruction.^{19} As an instructive example, we focus on the experimentally relevant situation where the molecules are all parallel to each other.

The paper is organized as follows: In Sec. II, we provide details on our model of the aggregate and the tip. We introduce a simple quantum mechanical model and define the absorption strength for the eigenstates of the aggregate for the near field of a tip. Here, we also show that the spectra can be understood from three basic functions that are independent of the orientation of the excitation-dipole. In Sec. III, we discuss as examples a one-dimensional regular chain and a two-dimensional array. More cases are presented in the supplementary material. In Sec. IV, we demonstrate that our simple model agrees well with a more complete description for the parameter regime considered. We conclude in Sec. V by summarizing our results and giving an outlook of the usefulness of the present study in the context of wavefunction reconstruction.

## II. SETUP AND THEORETICAL MODELLING

To show the relevant features most clearly, we use first a simple model where the molecules are treated as electronic two-level systems that experience an electromagnetic field that comes from a Hertz dipole. Hereby, we do not treat the incoming radiation explicitly. We describe this model in detail in the present section. In Sec. IV, we show that the spectra calculated via this simple approach are very similar to spectra that are obtained taking the incoming radiation, which shines on both the tip and molecules, into account explicitly.

### A. The molecular aggregate

For modeling of the aggregate, we use a widely adopted description.^{18,36} For each of the *N* individual molecules comprising the aggregate, we take two electronic states into account: the ground state $gn$ and the first excited state $en$, where the index *n* labels the monomers. The transition dipole between these two states is denoted by $\mu \u20d7n$. Initially, the aggregate is in the global ground state $gagg=g1\u2026gN$. For linear absorption, we are interested in states with one excitation. Using basis states $m=em\u220fn\u2260mNgn$, the excited state Hamiltonian for the system is written as

Here, *ɛ*_{m} is the transition energy of molecule *m* and *V*_{mn} is the transition dipole–dipole interaction, which we take as

with $R\u20d7mn$ being the distance vector from molecule *m* to molecule *n* and $Rmn=|R\u20d7mn|$. This form of interaction is sufficient for the purpose of the present work. The general features that we observe will not change when using an interaction that takes the transition charge density of the molecules into account in more detail or when considering modifications caused by the substrate.

From the Schrödinger equation $Hex\varphi (\u2113)=E\u2113\varphi (\u2113)$, one obtains *N* eigenenergies *E*_{ℓ} with corresponding eigenstates (which are labeled by *ℓ*),

In general, the coefficients $cm(\u2113)$ are complex numbers, but in the present case of a real Hamiltonian, they can be chosen to be real.

### B. Form of the near field and absorption

In our theoretical modeling, we take the essential features of the polarization of the incoming light into account. To this end, we assume that the created near field can be described by a field distribution that corresponds to the radiation stemming from a dipole. We characterize the dipole by its strength and orientation. We also ignore for the moment back-coupling of the molecular transition dipole to the tip. In a previous study,^{18} we have found that this should be fine for the typical strength of transition dipoles when the apex diameter is smaller than 5 nm and the tip is in the order of 1 nm away from the molecules. In Sec. IV, we confirm that this is true also for the present setup. We model the electromagnetic field stemming from the tip as a Hertzian dipole located at position $R\u20d7dip$ a few nanometers (molecule-tip distance + radius of the tip) above the aggregate. The resulting field has a strong spatial variation on the scale of a few nanometers in the plane of the aggregate. The resulting electromagnetic radiation of frequency *ω* with an electric field component that at time *t* and position $r\u20d7$ is given by

With the argument $R\u20d7dip$, we indicate the explicit dependence of the electric field on the position of the Hertz dipole, to which we also refer to as the “position of the tip.” In the near-field zone, the field of the Hertzian dipole with dipole moment $d\u20d7$ located at $R\u20d7dip$ can be written as

We have here not explicitly included the factor 1/4*πϵ*_{0}, which contains the dielectric constant *ϵ*_{0}. Note that we can write $E\u20d7(r\u20d7;R\u20d7dip)=E\u20d7(r\u20d7\u2212R\u20d7dip)$. The direction $d\u20d7$ of the tip-dipole can be affected by the polarization and propagation direction of the incoming light beam. It is convenient to use spherical coordinates and write [see also Fig. 2(b)]

with

The absorption strength from the ground state to the eigenstate $\varphi (\u2113)$ can be written as^{18}

This expression is valid for fields that vary strongly over the extent of the aggregate but only very weakly over the extent of a single molecule (see the supplementary material of Ref. 18).

From Eq. (10), the two-fold importance of the molecular transition dipoles is visible: First, they determine the interaction with the electromagnetic radiation (the strength is given by the projection of the electric field on the transition dipole); second, the mutual interaction is caused by the transition dipoles (this interaction has a strong dependence on the mutual orientations of the transition dipoles and on the orientations with respect to the distance vector between the molecules).

### C. Separating angular and spatial dependence

with functions $Gj(\u2113)$ that are independent of the angles *θ*_{dip} and *ϕ*_{dip} and contain the molecular orientation and the wavefunction coefficients,

with $R\u20d7dipm=R\u20d7m\u2212R\u20d7dip$, and $(a\u20d7)j$ denotes the *j*th component of the vector $a\u20d7$.

Since the functions Φ_{j} are linearly independent, the prefactors $Gj(\u2113)$ are the basic functions from which the spectra are composed. We will show this exemplarily in Sec. III. Sometimes, we drop the eigenstate-index *ℓ* for brevity.

## III. EXAMPLES

### A. Arrangement and overview

In this work, we consider molecular arrangements where the transition dipole moments of all molecules are parallel and have equal amplitudes (i.e., $\mu \u20d7n=\mu \mu \u0302$, with $\mu \u0302$ being the direction of the transition dipoles). Such arrangements have been found in several experiments.^{11,12,16,17} For parallel transition dipoles, we write [see also Figs. 2(a) and 2(b)]

The interaction simplifies to $Vmn=\mu 2Rmn31\u22123\u2061cos\alpha mn$ with *α*_{mn} being the angle between the molecular dipole $\mu \u20d7$ and the connecting vector $R\u20d7mn$ between molecules *m* and *n*. We use the following notation for the components of the position vectors of the monomers: $R\u20d7m=(xm,ym,zm)$. Similarly, the components of the tip-dipole are written as $R\u20d7dip=(xdip,ydip,zdip)$. In the following, we focus exemplarily on a one-dimensional chain of regularly spaced parallel molecules, as sketched in Fig. 2(a). Experimentally, such one-dimensional chains have been realized, e.g., in Refs. 12 and 16. For definiteness, we place *N* molecules on the *x*-axis of the coordinate system such that they are symmetrically arranged around zero: $R\u20d7m=(\Delta (m\u2212N/2),0,0)$. Here, Δ is the center-to-center distance between two neighboring molecules. For this regular one-dimensional chain, the eigenfunctions $\varphi (\u2113)$ and eigenenergies *E*^{(ℓ)} are known analytically (see Appendix B). While the eigenenergies depend on the angle *ϕ*_{Mon} of the molecules with respect to the *x*-axis, the eigenstates do not depend on it, i.e., in this case, the wavefunction coefficients from which the spectra are derived are independent of the orientation of the molecules. However, the orientation enters directly the spectra in a pronounced way as can be seen in Eq. (12). For each tip position, the absorption spectra consist of several peaks located at the eigenenergies *E*_{ℓ}. When the individual absorption peaks can be well resolved (as has been assumed in the formulas given above), one can then plot the spectra for each eigenenergy (i.e., for each *ℓ*) as a function of the position of the tip. In the following, we thus present the spatially resolved spectra for individual eigenstates $\varphi (\u2113)$, which are obtained by numerical diagonalization of the Hamiltonian (1). Here, we keep the distance from the tip to the surface constant and vary the tip position in a plane parallel to the surface. In the following, we will investigate how these position dependent spectra depend on the direction of $d\u20d7$. For all examples presented in this work, we use *z*_{dip} = 3Δ. Within the simple approach of Sec. II, one can conveniently introduce dimensionless units. We choose the distance Δ between neighboring molecules as the unit of distance and the magnitude of the molecular transition dipole *μ* as the unit for dipoles. Then, for example, energy is expressed in units of *μ*^{2}/Δ^{3} and the *G*-functions are in units of *μ*/Δ^{3}.

In Fig. 2(c), we show the dependence of the absorption strength [Eq. (10)] on the tip position and tip-dipole direction for the eigenstate that has no nodes for a linear chain with five molecules and *ϕ*_{Mon} = 45°. Each panel corresponds to a different orientation of the dipole $d\u20d7$. The respective angles *θ*_{dip} and *ϕ*_{dip} are indicated in Fig. 2 [all angles are defined in Fig. 2(b)]. Because we are interested primarily in the spatial variation of the spectra, we have normalized all spectra to have the same maximal value. In the upper row, *θ*_{dip} = 0, which corresponds to $d\u20d7$ perpendicular to the surface. Thus, for all angles *ϕ*_{dip}, we have the same spectrum.

For the other eigenstates, similar patterns as those in Fig. 2 can be observed (see Sec. 2.2 of the supplementary material). The main difference is the appearance of additional nodal lines because of the increasing number of sign changes of the eigenstate coefficients. The spectra for the state considered in Fig. 2 are actually very similar to the ones of a single molecule; the main difference is that the features are more elongated along the *x*-axis compared to the single molecule case (see Sec. 1 of the supplementary material, where a single molecule is discussed).

We note that for many cases the maximal intensity of the spectra is not located on the *x*-axis where the molecules are located; thus, measurements at finite *y*_{dip} values will provide relevant information.

In the supplementary material (Sec. 3), we also show plots similar to Fig. 2 for the case of a two-dimensional arrangement of parallel molecules. We observe a similar pattern as for the one-dimensional chain.

### B. Separating angular and spatial dependence

Since the absorption spectrum for an arbitrary orientation of $d\u20d7$ can be constructed from the three *G*-functions, it is instructive to take a closer look at these functions. Let us first make a few general remarks. The function $Gz(\u2113)$ is the only function appearing for *θ*_{dip} = 0, i.e., for an excitation-dipole perpendicular to the *x*–*y*-plane. This corresponds to the situation that we considered in Refs. 18 and 19. The function $Gx(\u2113)$ represents the cases *θ*_{dip} = 90° and *ϕ*_{dip} = 0, i.e., for an excitation-dipole parallel to the *x*-axis. Finally, $Gy(\u2113)$ is the only function appearing for *θ*_{dip} = 90° and *ϕ*_{dip} = 90°, i.e., for an excitation-dipole parallel to the *y*-axis.

For the case discussed in Fig. 2, the functions *G*_{x}, *G*_{y}, and *G*_{z} are shown in panel (d). In contrast to the spectra shown in panel (c), which are positive because of the modulus squared in Eq. (10), the three basic functions still have positive and negative parts. After taking the modulus square, according to Eq. (11), the corresponding spectra of panel (c) are recovered. In the supplementary material (Sec. 2.2), the three *G*-functions of the other eigenstates are also shown.

For a linear chain with parallel molecules, the wavefunction coefficients do not depend on the orientation of the transition dipoles *ϕ*_{Mon}, which only scales the eigenenergies. However, the *G*-functions and the position dependent absorption spectra have still the explicit dependence from the term $(\mu \u20d7)i$ entering in the definition Eq. (12). We have for this arrangement $(\mu \u20d7)x=cos\varphi Mon$, $(\mu \u20d7)y=sin\varphi Mon$, and $(\mu \u20d7)z=0$. Although this linear chain with parallel transition dipoles is a special case, it is nevertheless instructive to look at the dependence of *G*_{x}, *G*_{y}, and *G*_{z} on the angle *ϕ*_{Mon}. In Fig. 3, this dependence is shown for the eigenstate with the smallest number of sign changes of the coefficients (which we label by *ℓ* = 1); for the considered case of a linear chain, this state is located either at the highest or lowest energy and has the same sign for all the coefficients $cn\u2113$.

For all orientations, the pattern is point symmetric (or anti-symmetric) with respect to the center of the chain. Upon increasing *ϕ*_{Mon}, the pattern starts to rotate and to deform slightly. The pattern of *G*_{z} is rotated by the angle *ϕ*_{Mon}, and the pattern of *G*_{x} and *G*_{y} is rotated by *ϕ*_{Mon}/2. For the other eigenstates, we observe a similar behavior (see Sec. 2.3 of the supplementary material). The main difference is that the spectra of the other eigenstates possess more nodal lines because of the sign changes of the eigenstate coefficients.

The general features observed for the one-dimensional chain can also be observed for a two-dimensional lattice. As an example, we show the case of a 5 × 5 aggregate in Fig. 4. Here, the molecules are placed on positions (*n*_{x}Δ, *n*_{y}Δ, 0), where Δ is the distance between neighboring molecules and *n*_{x} and *n*_{y} take the values −2, −1, 0, 1, 2. The molecular transition dipoles are oriented by 45° with respect to the *x*-axis. This arrangement is also indicated by the bars. The results for a few eigenstates are shown in different rows of Fig. 4, with the corresponding eigenstate coefficients present in the right panel of each row (the squares correspond to the respective monomers of the lattice). In the upper row, we show the eigenstate that is closest to one where all coefficients have the same sign. The three *G*-functions then look indeed very similar to the ones of the chain (for *ϕ*_{Mon} = 45°). For eigenstates that have more pronounced sign changes of the eigenstate coefficients (shown in the other rows), the pattern of the *G*-functions becomes more complicated; in particular, they now have many nodal lines. Nevertheless, one still sees similarities to the first row. For example, the pattern of *G*_{x} can be obtained from that of *G*_{y} by a rotation of 90° and then mirroring along the *x*-axis.

## IV. VALIDATION OF RESULTS BY TAKING THE INCOMING FIELD EXPLICITLY INTO ACCOUNT

In this section, we now take the incoming electromagnetic field explicitly into account. To do so, we apply a self-consistent field approach, similar to what has been used in Ref. 18. To this end, we model the tip as a polarizable sphere. Both the tip and molecules are described by point dipoles and characterized by their respective polarizability tensors. In this respect, we treat the molecules and tip on the same footing. We calculate the extinction coefficient by^{37}

Here, the index *n* = 0 refers to the tip-dipole and the indices *n* = 1, …, *N* refer to the molecules. The polarizations $P\u20d7n$ are obtained from solving the coupled set of equations,

with the transfer tensor

As before, $R\u20d7mn$ is the separation vector between *m* and *n*; with $I\u2194$, we denote the identity tensor. The use of Eq. (16) for the interaction between the particles is equivalent to using the point dipole–dipole interaction for *V*_{mn} as given in Eq. (2). The polarizability tensor for the m-th molecule is taken as $\alpha \u2194m(\omega )=\mu \u20d7m\u2297\mu \u20d7m\omega \u2212\omega m+i\gamma m,$ where decoherence is taken into account via the parameter *γ*_{m}. As before, $\mu \u20d7m$ and *ω*_{m} = *ϵ*_{m}/*ℏ* are the transition dipole and the transition frequency of molecule *m*. To establish a connection to the approach used above, the molecular polarizability has to be chosen as an infinitely sharp resonance, i.e., *γ*_{m} → 0. For more extended discussions, see, for example, Refs. 38 and 39. For the nanotip-dipole, we use the polarizability of a sphere with radius *a*_{r}, which is^{40},$\alpha \u21940(\omega )=\alpha \u2194tip(\omega )=ar3\u03f5(\omega )\u2212\u03f5env\u03f5(\omega )+2\u03f5envI\u2194,$ where *ϵ*_{env} is the dielectric constant of the surrounding medium (we use *ϵ*_{env} = 1). The complex dielectric function, *ϵ*(*ω*), is evaluated according to a generalized Drude model with finite-size effect correction for particles with small radius (below ∼5 nm), $\u03f5(\omega )=\u03f5b+\omega p2\omega (\omega \u2212i\gamma p)\u2212\omega p2\omega (\omega \u2212i\gamma p\u2212ivF/ar).$ Here, *ϵ*_{b} is constant, and *ω*_{p}, *γ*_{p}, and *v*_{F} are the plasma frequency, Ohmic damping constant, and Fermi velocity in the bulk material of the tip, respectively. For the parameters appearing in this formula, we use the same values as in Ref. 18, which are *ϵ*_{b} = 9.0, *ω*_{p} = 7.26 · 10^{4} cm^{−1}, *γ*_{p} = 400.0 cm^{−1}, and *v*_{F} = 1.39 · 10^{8} cm/s. Different from the treatment of Sec. III B, now there are no simple scaling relations and we present results for specific choices of parameters; in particular, we fix the molecular transition frequency, the magnitude of the transition dipole, and the distance between molecules.

To connect to the results of Secs. II and III, we apply the following procedure. We scan along both the frequency and the position. From these scans, we can determine the positions and width of the molecular resonances. For the parameters used in the present article, where the sphere has a small diameter (∼5 nm) and is rather far away from the molecules (∼1 nm), these resonances are located very closely to the transition energies obtained from diagonalizing Eq. (1) and are quite narrow (a few wavenumbers). In Fig. 5(a), we show such scans for various positions of the tip. Once clearly sees well separated resonances with intensities that depend strongly on the position of the tip. The vertical dashed gray lines indicate the eigenenergies obtained from the diagonalization of the Hamiltonian (1). We then evaluate, for these eigenfrequencies, the extinction coefficient as a function of the *x* and *y* position of the tip. We average for each position over a frequency range of 6 cm^{−1}. Since, in Eq. (10), the absorption strength is calculated with a field coming from the tip alone, in the next step, we aim at removing the direct contribution of the incoming light. We do this by subtracting from *C*_{ex} that is obtained with a tip (∼1 nm) above the surface a reference extinction $Cexref$ where the tip is far away from the surface (in our calculation ∼10 nm). For the two energetically lowest resonances, such plots are shown in the left column of Fig. 5(b). As comparison in the right column, we show the position dependent absorption for the corresponding eigenstate calculated using Eq. (10) where we associate the direction of the Hertz dipole $d\u20d7$ with the direction of the polarization of the incoming electromagnetic field $E\u20d7inc$. As one can see, there is nearly perfect agreement. We have found similar agreement also for the other angle combinations and for the other eigenstates. This demonstrates that the simple formalism of Sec. II is indeed meaningful despite its many approximations. However, one has to keep in mind that the quality of the results of Sec. II starts to break down for stronger coupling between the tip and molecules, e.g., for smaller tip-molecule separation or larger tips.

## V. CONCLUSIONS

We have investigated the dependence of near-field spectra of molecular aggregates on the direction of an “excitation-dipole.” A particular focus was spatially resolved spectra for the fixed frequency of the electromagnetic radiation. These frequencies were chosen such that they correspond to transitions between the ground state and single exciton-states of the molecular aggregate. We found that there is a strong dependence of these spectra on the direction of the “excitation-dipole.” For general arrangements of the molecules that form the aggregate, we have seen that the spectra can be understood from three (in general complex) functions that only depend on the position, but not on the direction of the “excitation-dipole.” The complete spectrum is then obtained from a sum of these functions where each function is weighted by a factor that depends on the orientation of the “excitation-dipole.” Exemplarily, we have then focused on one- and two-dimensional aggregates where the molecules are all parallel to each other and are placed on a square lattice. Here, one finds very distinct spatial spectra for different angles of the molecular orientation with respect to the axis of the lattice. It would be interesting to investigate in the future how these spectra look for more complicate arrangements, such as “herringbone” structures.^{41} We have seen that from the near-field spectra one can on the one hand access states that are dark in far-field spectroscopy; on the other hand, one can obtain information on the position and orientation of the molecules inside the aggregate.

Furthermore, in Ref. 19, we have shown that for the special case of $d\u20d7$ along the *z*-axis (*θ*_{dip} = 0), one can obtain the wavefunctions (i.e., the coefficients $cn(\u2113)$) from the position dependent near-field spectra, for the case of one-dimensional and two-dimensional aggregates, even when there is considerable noise on top of the spectra. Note that the spectra in this case are given by $|Gz(\u2113)|2$. As we have seen above, for different orientation of the tip-dipole, the spectra look quite different, which is caused by the different weighting of the coefficients $cn(\u2113)$ in Eq. (12). Thus, each spectrum contains additional information that can be used to improve the reconstruction. In Appendix A, we discuss an example of a minimal set of dipole orientations $d\u20d7$ that allows one to obtain at each tip position $R\u20d7dip$ the values $Gx(\u2113)(R\u20d7dip)$, $Gy(\u2113)(R\u20d7dip)$, and $Gz(\u2113)(R\u20d7dip)$ up to a global sign. Thus, we expect that the reconstruction of the wavefunction works even better than for the case *θ*_{dip} = 0 alone, as considered in Ref. 19. Therefore, the results presented in the present work will be useful for the understanding of the collective eigenstates of molecular arrangements. We want to point out that in several cases there is a nodal line directly on the molecular chain. For example, for *G*_{z}, *ϕ*_{Mon} = 90°. Thus, even in the case of a one-dimensional arrangement, it is useful to scan along a two-dimensional plane with the tip. Finally, let us comment on approximations made in the present work: First, we did not explicitly consider the tip but simply assumed that the tip creates a dipolar field that acts on the molecules. The main reason for this was that it allows a clearer discussion of the relation between wavefunction coefficients and spectra. We have verified our results by using local field theory as in Refs. 18 and 19, where we model the tip as a single isotropic dipole that is irradiated by the external field. Clearly, when one wants to treat a specific experimental situation, one would have to model the tip in more detail. Then, also one should take into account that the incoming light not only hits the tip but also the molecular aggregate and the substrate. We expect that in this case the spectra still show the basic features discussed in the present work. To extract information about the molecular system, one could, for example, use experimentally modulation techniques.^{24,35} From the theoretical side, one could use machine learning approaches, such as the one used in Ref. 19, to extract the wavefunction directly from these spectra. Another aspect that we did not touch in the present work is the coupling to molecular vibrations and the environment. We treated the molecules simply as a two-level system, where the two states belong to the ground and first excited electronic molecular states. Typically, however, in both electronic states should be described by Born–Oppenheimer surfaces that host vibrational states. We neglected these vibrations because of two reasons: First, in the experimental work on molecular aggregates on a dielectric surface,^{10,11} it was found that excited vibrations are energetically well separated from the dominant 0–0 transition, and the dipole–dipole coupling is so small that to a good approximation, one can model the system as dipole–dipole interacting two-level systems. Second, before tackling the much more complicated case of exciton vibrational coupling, one first needs a good understanding of the purely electronic case, which is also numerically fast to calculate. Similar to internal vibrations, there will also be coupling to environmental degrees of freedom. At low temperature, the resulting effect on the molecular spectra is very small as has been found experimentally.^{10,11} However, at higher temperature, pronounced broadening of the spectra^{10} and decoherence^{11} has been observed. In such situations, it will also no longer be obvious how to assign wavefunctions to the spatially dependent spectra. However, we expect that for a broad range of parameters, the basic pattern obtained from our simple model will be visible and will help us to understand the more complicated observations.

## SUPPLEMENTARY MATERIAL

See the supplementary material for the other eigenstates and additional arrangements.

## ACKNOWLEDGMENTS

A.E. acknowledges support from the DFG via a Heisenberg fellowship (Grant No. EI 872/5-1). We thank M. T. Eiles for many helpful comments.

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: EXTRACTION OF THE *G*-FUNCTIONS

Because of the modulus appearing in Eq. (11), the *G*-functions cannot be simply extracted from an absorption spectrum. Writing Eq. (11) as $A(R\u20d7dip;\theta dip,\varphi dip)=d2\u2211ij(Gi(\u2113))*Gj(\u2113)\Phi i\Phi j$, one sees that from suitable choices of the angles *θ*_{dip} and *ϕ*_{dip}, it is possible to find the products $(Gi(\u2113))*Gj(\u2113)$. For example, *θ*_{dip} = 0 provides $|Gz|=A(0,\varphi dip)/d$. Similarly, *θ*_{dip} = *π*/2 and *ϕ*_{dip} = 0 provide directly $|Gx|=A(\pi /2,0)/d\theta dip=\pi /2$, and *ϕ*_{dip} = *π*/2 provides directly $|Gy|=A(\pi /2,\pi /2)/d$. Thus, at each tip position $R\u20d7dip$, the magnitude of *G*_{x,y,z} is determined, but not the sign.^{42} The relative signs can be obtained from $GxGz+|Gz|2=2A(\pi /4,0)/d2$ and $GyGz+|Gz|2=2A(\pi /4,\pi /2)/d2$.

While it is possible to obtain the magnitude and the relative signs of the *G*-functions at each point $R\u20d7dip$, it is not directly possible to obtain the relative signs at different tip positions. One has to keep in mind that the sign of the *G*-functions depends on the chosen overall sign of the wavefunction. If this is fixed, then all signs of the *G*-functions are also fixed.

### APPENDIX B: EIGENFUNCTION OF THE CHAIN

In this appendix, we discuss some properties of the 1D chain, which are helpful to better understand the numerical observations provided above. We first note that for the considered chains, the dipole–dipole interaction simplifies to

To obtain simple analytical results, we restrict to only nearest neighbors’ interaction that is $V0=\mu 2\Delta 3(1\u22123\u2061cos2\varphi Mon)$.

The eigenenergies and eigenstates are given by

and

From Eq. (B3), one sees that the eigenstates are independent of the value of the interaction and in particular independent from the orientation of the molecules. Therefore, it is convenient to label the eigenstates according to the number of sign changes of the coefficients (B3), i.e., for *ℓ* = 1, there is no sign change, and for general *ℓ*, there are *ℓ* − 1 sign changes. After we have thus fixed the labeling of the wavefunctions, we see that the eigenenergies are for positive values of *V*_{0} decreasing with increasing *ℓ*, and for negative values of *V*_{0}, they are increasing with increasing values of *ℓ*. The sign of *V*_{0} depends on the angle *ϕ*_{Mon}. *V*_{0} is positive for *ϕ*_{Mon} ≳ 54° and negative for *ϕ*_{Mon} ≲ 54°.

## REFERENCES

Note that in our case the functions *G* are real valued because the wavefunction coefficients can be chosen to be real valued.