Collective Strong Coupling Modifies Aggregation and Solvation

Intermolecular (Coulombic) interactions are pivotal for aggregation, solvation, and crystallization. We demonstrate that the collective strong coupling of several molecules to a single optical mode results in notable changes in the molecular excitations around a single perturbed molecule, thus representing an impurity in an otherwise ordered system. A competition between short-range coulombic and long-range photonic correlations inverts the local transition density in a polaritonic state, suggesting notable changes in the polarizability of the solvation shell. Our results provide an alternative perspective on recent work in polaritonic chemistry and pave the way for the rigorous treatment of cooperative effects in aggregation, solvation, and crystallization.

The coherent interaction of molecules with confined optical modes leads to hybrid light-matter states called polaritons [1][2][3].In the strong coupling regime, usually achieved by coupling several molecules to the optical cavity, experiments show significant modifications of chemical properties [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19].Due to their delocalized nature, polaritons in quantum optics are studied from a collective perspective where the molecules are modeled as few-level systems that indirectly interact solely through the photon field [20][21][22][23][24][25][26][27].Yet, chemistry is governed by local interactions, and the molecular complexity requires refined chemical methods to analyze processes like reactions.The behavior of molecules is hence susceptible to their immediate surroundings.For example, their spectral absorption and emission can be altered by their interplay with the solvent (solvatochromic) or by close proximity to identical molecules (concentration-dependent). Experiments under strong coupling often use organic molecules that are sensitive to their surroundings, show intense excitations, or even form aggregates [6-9, 11, 28-41].The chemical environment can then play an active role and exert significant influence, as emphasized by several vibrational strong coupling (VSC) experiments showing modification of assembly [42][43][44] and reactivity [14,45,46].For a clear chemical understanding of polaritonic processes, we must then investigate the coexisting roles of the molecule, the solvent (chemical environment), and the cavity (optical environment).To this end, ab initio quantum electrodynamics (QED) merges the knowledge of quantum optics and quantum chemistry to explore single-molecule effects retaining the chemical complexity [47][48][49][50][51][52][53][54][55][56][57][58][59].Even if collective states can show larger contributions from selected molecules [54,60], local changes tend to reduce due to collective delocalization and an increasing number of quasidark states [61].While recent work indicates a larger relevance of dynamic electronic polarization [59,62,63], it remains puzzling how the polaritons' collective nature enters chemical processes.The role of changes in the chemical environment (solvents, aggregates) under strong coupling has been widely disregarded up to this point.
In this letter, we add another facet to the understanding of polaritonic chemistry by reintroducing the chemical environment.Using QED coupled cluster (QED-CC) [55], which is at present the most accurate approach for medium-sized molecular polaritons, we study the yet unexplored competition between photon and intermolecular interactions.We extract the single-molecule response in collective ensembles and investigate the microscopic changes of aggregates or solvation systems in polaritonic chemistry.We demonstrate local response modifications unambiguously distinguishable from collective delocalization and arising from the interplay between Coulomb and transverse fields.Moreover, our results show a slower decrease than conventional local effects in polaritonics, i.e., not with the total number of coupled molecules but with the number of affected solvation shells.As a result, the immediate surroundings of the impurity, representing a solute, nucleation, or reaction center, undergo notable changes.Structural changes in the chemical environment can then influence the dynamics of the impurity.In the following, solvation, aggregation, and nucleation are considered synonymously due to their conceptual similarity.
Modelling Solvation and Collective Coupling -Molecular properties are significantly modified by chemical environments such as polymer matrices, solvents, or aggregates.The interaction of a solute with its surrounding molecules leads to changes in the ground and excited electronic densities, which affect energy levels, molecular geometries, and even chemical reactivity.In addition, the solvent also directly contributes to properties such as spectroscopic signals.To study the impact of polaritons on the complex chemical response, we must simultaneously include the photon coupling and the Coulomb interactions in the solute-solvent system.Therefore, we focus on a system with non-negligible intermolecular interactions as a prototype for the Coulomb forces in supramolec- ular structures.Our model system is a chain of N hydrogen molecules that form an H-aggregate illustrated in Figure 1.An impurity (P), representing a solute or reactive molecule, is introduced by stretching the bond of the central hydrogen, and we investigate different intermolecular distances to alter the Coulomb forces.Figure 1 shows the absorption spectrum of the undressed (out-ofcavity) chain for different N and 5 Å intermolecular separation.The global response of the system is governed by the transition density ϱ j0 (x, y, z), which is for instance connected to the optical properties via the transition dipole moment µ j0 = R 3 rϱ j0 (x, y, z)d 3 r.To characterize the local behavior of each dimer, we integrate ϱ j0 (x, y, z) perpendicularly to the excitations transition dipole moments (that is, in the Y and Z directions) around P, its nearestneighbor A1 (first solvation shell), and its next-nearest neighbor A2.The integrated local transition densities I dy ϱ j0 (x, y, z) (I = P, A1, A2) for the excited states E1 and E2 are depicted in Figure 1.The excitation E1 is mainly localized on P and A1 with anti-aligned local transition dipoles, while E2 is delocalized and shows an all-parallel transition moment alignment.In Figure 2, we show the polaritonic absorption spectrum with the photon resonant to the (isolated) unperturbed H 2 excitation (with light-matter coupling strength λ = 1/(ε 0 V ) = 0.005 a.u.).Three polaritonic branches, the lower (LP), middle (MP), and upper polaritons (UP) emerge.The electron and electron-photon correlation is accurately described by the ab initio QED-CC ansatz wave function, thus capturing any modification of inter-(dispersion) and intra-molecular forces [46,55,56,64,65].Nevertheless, due to the low coupling strength, the molecular ground state is here essentially unaffected by the embedding in the optical environment, and polaritonic effects ef-fectively arise from collective coupling.The dimers are indirectly coupled through the cavity, establishing an interplay between longitudinal (Coulomb) and transverse (photon) fields.The simplicity and engineerability of our model allow us to separate intermolecular and cavity-induced effects, providing a satisfactory proof-of-concept model for investigating microscopic molecular responses in complex chemical and photonic environments.
Modification of Solvation -In Figure 2, we plot the LP local transition densities for P, A1, and A2.As N increases, the A1 local transition density changes sign, i.e., the transition dipole of A1 aligns with the other molecules.This effect requires that both the (undressed) excitations E1 and E2 contribute to the LP state, and hence occurs around the avoided crossing domain illustrated in Figure 2.While the alignment clearly originates from collective effects, its impact is evidently localized in the first solvation shell and requires an explicit treatment of the solute-solvent interface.The local change arises from competition between the short-range intermolecular forces, which pattern E1 and E2, and the collective interaction to the optical mode, which tends to align the molecular transition moments (see also Sec.III of the SI).Indeed, for larger intermolecular separations, the LP local moments are aligned (see Sec. V of the SI), confirming that the sign change stems from the interplay between Coulomb and photon fields.Which effect prevails depends on the setup, specifically on the excitation energies (solute-solvent system), the chain length N (collective coupling), and the light-matter coupling strength λ = 1/(ε 0 V ) (optical device features), and the intermolecular forces (solute-solvent system and intermolecular separation).
The number of solvent molecules heavily outweighs the concentration of solutes.It is therefore instructive to ver- ify how our observation scales in the thermodynamic limit N e , V → ∞, N e /V = const., i.e., with fixed N e /V but increasing particle number N e and quantization volume V (controlling the fundamental coupling strength).To this end, we resort to a simplified Tavis-Cummings (TC)-Kasha [66][67][68][69] model (parameterized according to the electronic CC calculations), with a nearest-neighbor transition dipole-dipole coupling between the molecular excitations (see Sec. III-C and V-A of the SI).The total number of molecules N e = N rep × N is given by the number of molecules N per aggregate (i.e., the same number used in Figure 1 and 2), and how many times the system appears in the cavity N rep .In Figure 3(a), we report the eigenvector coefficients for a single P and A1 to quantify how the individual molecules contribute to the collective excitations of the TC-Kasha model.The ratio N e /V is fixed, i.e., the Rabi-splitting is kept constant while increasing N .The simulations for N rep = 1 and 4 show analogous results and thus prove that the described behavior is qualitatively resistant to the thermodynamic limit.Still, the single-molecule contribution (i.e., the TC-Kasha eigenvector coefficients of a single dimer) decreases by approximately 1/ N rep .Similar conclusions are drawn from Figure 3(b-c), where for each panel λ and N rep are fixed while increasing N .This allows us to investigate the effect of changing the number of coupled molecules and thus the collective effects.Increasing the number N − 1 of coupled solvent molecules at fixed N rep and λ shortens the critical intermolecular distance at which the A1 coefficient shows a change.The panels 3(b1) and (b2), computed at different couplings, show similar behavior and order of magnitude, but the N for which the transition moment of A1 changes sign is shifted to longer chains.We also confirm this trend by simulating two H 2 chains in our QED-CC calculations (see Sec. V-A of the SI).Comparing Figure 3(b1) and (c) shows the effect of increasing N rep /N while rescaling λ by 1/ N rep , which thus retains the overall collective coupling.The panels show the same qualitative behavior, with coefficients rescaled by 1/ N rep , as in 3(a).Therefore, Figure 3 demonstrates that the magnitude of the change in the environment response for fixed N rep /V decreases with the number of activated species P, which is significantly smaller than the number of coupled molecules N e , and P does not need to be strongly coupled to the device.In other words, if the solvent exists in vast excess as in experiments showing modifications of crystallization [42][43][44] and ionic conductivity [70], each solvation shell will experience a considerably larger effect from the cavity than the bulk of the solvent molecules.Moreover, keeping the solute-solvent ratio fixed (increasing N rep ) shortens the critical length N as seen from Figure 3(b) and (c).Therefore, the intermolecular distance and chain length N for which the transition moment of A1 changes sign depends solely on the solvent strong coupling defined by the ratio N rep × (N − 1)/V .
So far, we considered all the molecules at the same distance.However, analogous results are found for a small H 2 cluster, with the other hydrogens at larger distances.This is illustrated in Figure 4 for (H 2 ) 3 with intermolecular distance 3.7 Å, and slightly increased λ = 0.01 a.u. to compensate for the larger blueshift of E1.The remaining N -3 hydrogens are placed at 20 Å distance.The results show the same behavior as in Figure 2 and emphasize the locality induced by the short-range Coulomb interactions, clearly highlighting the interplay between the chemical surroundings and the collective optical dressing.
The TC-Kasha model successfully reproduces the modifications in the dynamic response of the first solvation shell.Since the model is parametrized with bare (electronic) transitions, it stresses the importance of collective behavior, which alters the LP microscopic response near P once a critical cooperative coupling strength is reached.This is consistent with experimental setups where the same cavity (i.e., λ fixed) is employed at different molecular concentrations (here, the solvent molecules N − 1).c) Same as (b) but for Nrep = 4. Compared to (b1), the overall collective coupling is maintained, but due to the larger number of solute molecules, the coefficients are rescaled by 1/ Nrep.The impurity concentration P then defines the decay of the wavefunction coefficients.Notice that P does not need to be in strong coupling.See Sec.V-A of the SI for more results.
Nevertheless, the TC-Kasha model has apparent limitations in that it will perform poorly in modelling the electronic dynamic of more complex systems and it cannot account for non-linear effects included in the QED-CC calculations, i.e., a self-consistent and more realistic treatment is likely to give similar or more favorable scaling.It should also be noted that the strength and structure of solvation vary widely with the type of intermolecular interaction, such that the onset of this effect can change dramatically with the specific choice of solvent and solute.Quantum chemistry then provides a flexible tool-set to explore such solvent effects beyond the heavily simplified dipole-dipole picture [71][72][73][74][75][76].
Disorder and Resonance dependence -While aggregation, solvation, and crystallization imply a somewhat structured environment, disorder (in the form of orientation or inhomogeneous broadening) is inevitable in chemistry under standard ambient conditions.Inhomogeneous broadening and the molecular orientation relative to the cavity polarization, for all solvent molecules besides A1 (long-range disorder), is included in our TC-Kasha model by sampling the excitation energies from a Gaussian distribution and assigning random H 2 orientations (see SI Sec.V-A for details).As demonstrated in Figure 5 and Sec.V-A in the SI, the long-range disorder merely increases the required number of molecules N , but the overall solute-solvent effect remains unaffected.At the same time, local disorder severely impacts the intermolecular interactions since the Coulomb forces are strongly directional, i.e., other forms of aggregation will impact the solute differently.In the SI, we report the local transition densities for a J-aggregate (headtail) (H 2 ) N configuration and show a sign flip for the local transition dipole of A1 in the MP instead of the LP.The modification of the solvent-solute response is conceptually identical and merely moves to a different polaritonic state.This stems from the different transition moment patterns in the undressed excitations of the J-aggregate.Hence, we can expect that a realistic solvation system will exhibit the discussed effect partially in the LP and the MP.
Lastly, the solvent-solute polarization features a clear resonant behavior.When the cavity is detuned from the E2 excitation (entirely off-resonant or tuned to E1) the local changes are either negligible or indistinguishable from collective delocalization effects.If the field is tuned to E2 (instead of the excitation of an isolated unperturbed H 2 ), the results are in line with the foregoing discussion, although require a slightly higher N to achieve the sign flip in A1.These results are consistent with the position of the avoided crossing between the LP and MP, which depends on the photon frequency and determines the mixing of E1 and E2 (and thus the local impact of the strong coupling regime).
Conclusion -Using state-of-the-art ab initio QED-CC and quantum-optics models, we illustrate how the intermolecular forces interplay with collective coupling in cooperative systems and induce considerable changes in the dynamic polarizability of the first solvation shell.Such changes arise from a competition between the collective interaction with the optical mode and the local Coulomb interactions, which favors aligned (J-aggregate) or antialigned (H-aggregate) patterns in the transition moments.Which effect prevails depends on the local chemical environment around the solute and the collective coupling strength.Cooperative coupling to a solvent thus induces strong modifications in the dynamic polarization of the solute's first solvation shell.A change in the dynamic polarization of the first solvation shell is then likely to affect the solute dynamics, possibly leading to changes in solvent rearrangement, chemical reactivity, nucleation, aggregation, and ionic conductivity.Such a mechanism could be experimentally investigated using ultrafast UV electronic spectroscopy [40,41,77,78].Moreover, our results suggest that these effects feature a clear resonance dependence and will be relevant even in the thermodynamic limit and for disordered systems.Hence, our work provides an alternative perspective on cooperative strong coupling and the underlying mechanism that modifies chemistry via changes in the interaction within solvents and aggregates.Several VSC experiments have shown that the resonant coupling can induce modifications in solvent effects [14,45,46] and assembly [42][43][44], a feature proposed to emerge from collective coupling and structural reorganization rather than significant individual molecular changes [79].Longitudinal interactions are ubiquitous and fundamental in understanding chemical effects, and the concepts developed in this letter can help and inspire further discussions on both electronic and VSC.While we fixed the nuclear positions in this work, motivated by a separation of time scales, it is apparent that vibrational strong coupling will require the inclusion of nuclear motion.We can therefore expect that changes in the first solvation shell might result in the reorganization of the solvation structure.These results also encourage the development and refinement of multiscale approaches [54,80] to extend our study to experimentally investigated systems.
The model system we focus on in this letter is a chain of N parallel equidistant molecules.The dimers are arranged either in an H-aggregate parallel-bond or a Jaggregate head-tail configuration as depicted in Figure 1, with constant intermolecular separation d.The central dimer, labeled P, is slightly stretched with a bond length of 0.78 Å, while the other chain dimers have a bond length of 0.76 Å.The chain lies in the XY plane, with the hydrogens aligned on X.The first and second adjacent dimers of P are labeled A1 and A2.The system is modeled at CCSD and QED-CCSD-2 level with the aug-cc-pVDZ basis set (see II for a brief description of the methods), and the light-matter coupling constant for the QED calculations is set to 0.005 a.u. with polarization along X.The stretched dimer has an (isolated and undressed) excitation energy of 12.337 eV, slightly lower than the unperturbed excitation of 12.498 eV.The transition dipole moment is along the bonds (X-axis).For each excitation, the transition densities ϱ j0 (x, y, z) are computed and saved in a cube file with a grid of 0.1 Å spacing.The local transition densities along X ρ I j0 (x) (direction of the transition moment) are computed by numerically integrating the cube on Z and Y around the hydrogen I = P, A1, A2 (within 2.0 Å from each atom) Figure 1.Configurations of the hydrogen aggregates.* Electronic address: henrik.koch@ntnu.no† Electronic address: christian.schaefer.physics@gmail.com All the calculations are computed on a local branch of the e T program [1].The Hamiltonian in the dipole approximation and second quantization is where b † (b) is the creation (annihilation) operator for the photon mode of frequency ω, p and q are indices for the one-electron orbitals, E pq and e pqrs are the standard one and two electron singlet excitation operators [2], λ is the light-matter coupling vector along the photon polarization, h pq and g pqrs are the one and two electron integrals, and d is the dipole operator.

II. QED-CC
The QED-CC ansatz is where T is the cluster operator and |HF, 0⟩ is the direct product of the QED-HF Slater determinant and the electromagnetic vacuum [3].The electronic cluster T e is an electronic excitation operator [2] where |µ⟩ is an excited electronic determinant.The cluster operator T p includes pure photonic excitations arXiv:2312.08814v1 [quant-ph] 14 Dec 2023 where we considered a single cavity mode.Notice that the bosonic space must be truncated up to a finite number n of photon occupations.Finally, the operator T int contains excitations of both matter and photons where for instance Notice that QED-CC is defined for the Hamiltonian transformed with the QED-HF coherent-state transformation U QED-HF [3] where ⟨d⟩ QED-HF is the mean QED-HF dipole moment.The QED-HF orbitals are obtained by minimizing the energy of the transformed Hamiltonian in Equation 11, with respect to the orbital coefficients.This leads to a modified Fock matrix such that Brillouin's theorem F ia = 0 determines the coefficients where F e pq is the standard HF matrix [2], and a, i refer to virtual and occupied orbitals.The QED-CC state is then obtained by projection of Equation 11 onto the space spanned by the electron-photon excitations |µ, n⟩ [3] The QED-CC dual state is defined as where tµn are the Lagrangian multipliers [2,3] obtained by the equation QED-CC is based on a hierarchy of approximations defined by a truncation of the cluster operators, such that the equations are projected only onto a subspace of the full Fock space [2,3].In this letter, we make use of the QED-CCSD-2 method, in which the electronic cluster includes only single and double excitations and the photonic operator includes single and double photonic creation operators The interaction cluster T int includes all the products between the operators in T p and T e .The excited states are defined in the equation of motion (EOM) framework via a CI-like parametrization, which includes the same states as in the cluster operator The left and right excited state coefficients are therefore obtained by diagonalization of the non-hermitian CCtransformed Hamiltonian [2,3] where From the left ⟨L k | and right |R k ⟩ eigenvectors of Equation 22 we define the transition and state CC densities from which we can compute excited and transition properties such as the oscillator strengths.Notice that since the CC theory is not hermitian, the left and right transition densities are not adjoints.

III. SIMPLE MODELS A. Jaynes-Cummings model
The coupling of an optical mode with a single excitation is often described by the Jaynes-Cummings Hamiltonian H JC , where the molecule is described as a two-level system coupled to the photon by an electric dipole interaction in the rotating wave approximation [4] (28) In Equation 28, b † (b) is the creation (annihilation) operator for the photons, |g⟩ and |e⟩ are the ground and excited molecular states, with energy difference ω mol (the ground state energy is arbitrarily set to zero), and g0 is the coupling strength where λ is the fundamental light-matter coupling strength directed along the polarization of the field and d ge is the transition dipole moment.For the Hamiltonian in Equation 28, because of the rotating wave approximation, the excitation number operator η where |n⟩ is the electromagnetic n-photon number state, commutes with H JC and therefore the Hamiltonian is block diagonal in the subspaces {|e⟩ |n − 1⟩ , |g⟩ |n⟩}.We then focus on the "vacuum Rabi splitting," with the subspace containing at most one photon excitation (excitation number η = 1).Then, in the resonant case ω mol = ω ph , the eigenstates of H JC are the cavity polaritons arising from the symmetric and antisymmetric combinations of the photon and the molecular excitation where |1⟩ and |0⟩ are respectively the one-photon and vacuum electromagnetic states.The upper and lower polaritons are separated by the vacuum Rabi splitting Ω while the ground state is the unchanged molecular ground state in the electromagnetic vacuum |g⟩ |0⟩.In the general subspace of excitation number η = η, the Rabi splitting scales as √ η.

B. Tavis-Cummings model
In the case of N identical molecules coupled to the same photonic mode, the TC or Dicke Hamiltonian reads [5,6] where the molecules are assumed to have the same excitation energy ω mol , the same coupling strength g0 , and the index k runs over the different replicas.We now assume the system to be resonant ω ph = ω mol = ω, and the Hamiltonian in the first excited manifold is The system is more easily diagonalized by using the collective state basis for matter The states |D j ⟩ |0⟩ are eigenvectors of H T C with eigenvalues ω and are referred to as dark states since they have no contributions from the excited photon field.The upper and lower polaritons |±⟩ are and have contributions from both the photon field and the N molecules, and are separated by a Rabi splitting that is, the Rabi splitting scales as the squared root of the number of coupled replicas.Notice that the sign of the transition dipole moment is arbitrary.If we flip the sign for one of the replicas, which we can assume without loss of generality to be the first one, the Hamiltonian reads The matrices in Equation 34 and Equation 40 are similarly equivalent, connected by the transformation Therefore, the new polaritonic states have the coefficients of |e 1 , g 2 , g 3 , . ..⟩ with a flipped sign compared to the other representation.This ensures that all of the transition dipoles constructively sum and not cancel with each other.Therefore, in the polaritons, all of the transition moments must be aligned to maximize the electric-dipole interaction with the cavity field.On the other hand, in the dark states, the phases of the molecular excitations are such that the transition moment sum vanishes.

C. Disordered TC
The molecular disorder affects both the molecular energies (inhomogeneous broadening) and the coupling strengths.Following Houdré et al. [7], the disordered TC Hamiltonian is and the polaritonic energies E are defined by the equation Each molecule is then endowed with a different excitation energy (inhomogeneous broadening) and coupling strength (orientational and positional disorder).If we disregard the inhomogeneous broadening of the excitations and set ω k = ω ph = ω, the orientational and positional disorder only affects the coupling strengths gk via the scalar product λ • d ge and the system have again N − 1 degenerate states of energy ω and two polaritonic excitations with a Rabi splitting Therefore, the Rabi splitting is determined by the quadratic average of the different single-molecule couplings, and the contribution of each excited molecule to the polaritons is weighted by its coupling [7].Notice that from Equation 44 we see that if all the molecules have the same module coupling |g k | ≡ g each contribution is weighted by gk /g = sgn(g k ), that is, it is weighted by the sign of the coupling strength so that all of the contributions from the transition moments add up.The Rabi splitting is smaller than in the absence of disorder, but the overall polariton effects are the same.
If the inhomogeneous broadening of the molecular excitation is introduced such that we can resort to a perturbative solution for the polaritonic energies as discussed in Ref. [7] ω The remaining N − 1 states are such that Therefore, for the photon coupling to an inhomogeneously broadened set of oscillators, the TC predicts that the Rabi splitting is unchanged and, to first-order, the eigenstates are the same as in Equation 44 [7].

D. TC analysis of the hydrogen chain
For the (H 2 ) N system under study, the TC Hamiltonian reads where the primed quantities refer to the stretched dimer and the photon is resonant with the unperturbed H 2 molecules.The eigenstates of H ′ T C are obtained by diagonalizing the matrix In Figure 2, we plot the coefficient of the photon and H 2 excited states that contribute to the LP and MP from equation Equation 50for 4 ≤ N ≤ 11.Notice that different scales are applied to each panel.In the LP, the perturbed and unperturbed H 2 excitations are in phase (same sign of the CI coefficients), and are in antiphase with the photon contribution.The coefficient of the impurity decreases with N and the photon character becomes more and more relevant, while the unperturbed dimer contribution reaches a maximum around N = 9 and then decreases again.Notice that, since in the TC the intermolecular forces are completely disregarded, each unperturbed H 2 will give the same contribution to the polaritonic states.The differences observed in the dimers A1 and A2 shown in the letter are due to the intermolecular forces.In the MP, the perturbed dimers are in antiphase with the photon and the unperturbed H 2 molecules.The polariton becomes more localized on the impurity while the photon and the impurity contribution approach zero.
IV. HYDROGEN AGGREGATES ELECTRONIC STRUCTURE PROPERTIES (NO QED)

A. H-aggregates
In Figure 3, we show the electronic (undressed) absorption spectra of the (H 2 ) N ≥4 H-aggregate (see Figure 1) for different intermolecular separations.As the system enlarges, the upper excitation E2 is blueshifted due to the intermolecular forces; therefore, the shift is less significant for large intermolecular distances.Notice that E2 gains oscillator strength with N , suggesting that all the unperturbed dimers contribute to this excitation.The excitation E1 is almost independent of the chain length, and it redshifts and loses intensity for larger longitudinal couplings (shorter intermolecular separations).Therefore, this suggests that E1 is mainly associated with the perturbed dimers P, and the nearest dimers produce an effective screening of its transition dipole moment.The system also possesses N −2 quasi-dark states, which gain some oscillator strength due to the intermolecular forces.The contribution of each dimer to E1 and E2 can be extracted from the transition densities integrated around each H 2 as explained in section I and section II.
In Figure 4 and Figure 5, we show the left and right local transition densities of E1 for P, A1, and A2 for 5 Å and 6 Å intermolecular separation, respectively.As for the absorption spectra of Figure 3, the local densities for E1 are independent of N and are mainly localized on P. Due to the intermolecular forces, an excitation of P induces a response on its neighbor dimers, with a less important effect for more distant dimers.As the intermolecular distances become larger, the excitation becomes more localized on P and less on the other dimers, as expected when decreasing the intermolecular couplings.Adjacent molecular transition dipoles are anti-aligned, partially screening the excitation, resulting in smaller oscillator strength than the isolated P. Notice that for A2, we only plot the local densities for N ≥ 6 since for N = 4, 5, its "first solvation shell" is incomplete.Therefore, its properties have significant variations from N = 4 to 6 due to local intermolecular forces.
In Figure 6 and Figure 7, we show the left and right local transition densities of E2 for P, A1, and A2 for 5 Å and 6 Å intermolecular separation, respectively.In this case, the local transition dipoles are all aligned, and all the dimers in the chain give a similar contribution.Therefore, the properties of E2 strongly depend on the length of the whole chain and the intermolecular forces acting on each dimer.With larger intermolecular distances, the contribution of P to E2 decreases, while the other dimers maintain a similar local transition density.1).The excitation is almost independent of N and mainly localized on P, A1, and A2 with decreasing contributions.Compared to Figure 4, the intermolecular forces are less relevant, and therefore, the excitation is more localized on P and less on A1 and A2.Still, the overall anti-alignment pattern of the molecular transition dipoles is maintained.Figure 6.Local left and right transition densities for E2 for an intermolecular separation of 5 Å in an H-aggregate configuration (parallel-bond arrangement, see Figure 1).The excitation blueshifts and gains oscillator strength with the chain length N and is delocalized over the whole aggregate with a larger contribution from the unperturbed molecules.The local transition dipole moments are all aligned.Figure 7. Local left and right transition densities for E2 for an intermolecular separation of 6 Å in an H-aggregate configuration (parallel-bond arrangement, see Figure 1).The excitation blueshifts and gains oscillator strength with the chain length N .It is delocalized over the whole aggregate, with a more significant contribution from the unperturbed molecules.The local transition dipole moments are all aligned, and compared to Figure 6, there is a smaller contribution from P and a similar contribution from the rest of the chain.For larger distances, that is, smaller intermolecular forces, the blueshift is less pronounced.

B. J-aggregates
Since the Coulomb forces are strongly directional, modifying the local arrangement severely affects the system's properties.In Figure 8, we show the absorption spectra for a head-tail configuration (J-aggregate, see Figure 1) for different intermolecular distances.Compared to the H-aggregate of Figure 3, the effects of the intermolecular forces are more pronounced, as seen from the more complicated structure of the spectra for 5 Å and 6 Å.For d = 5 Å and 6 Å, the N -2 quasi-dark state gained relevant oscillator strength, and the appreciable changes of E1 from N = 4 to N = 6 suggest that the dimer A2 and even father away give a significant contribution to E1.For larger intermolecular separations, E1 is almost independent of N as in Figure 3.The excitations E1 and E2 redshift for increasing chain length and decreasing intermolecular separation, while E2 blueshifts for the H-aggregate.Unlike the H-aggregate case, E1 gains oscillator strength from embedding in the H 2 chain, while E2 loses intensity with increasing intermolecular forces.Since the excitation energy E1 shows little dependence on the chain length, we deduce that it is still mainly localized on P, A1, and A2.E2 becomes more intense with N , suggesting that all dimers contribute to this excitation.The differences between the H-and J-aggregates result from the different effects of the local Coulomb (longitudinal) forces on the dimers.In Figure 9 and Figure 10, we report the left and right local transition densities of E1 for P, A1, and A2 for 6 Å and 8 Å intermolecular separation, respectively.Contrary to the H-aggregate arrangement, in this configuration, the local transition moments are all aligned, increasing the excitation's intensity.The densities are almost independent of the chain length N , and the contribution of the neighbor dimers A1 and A2 increases with decreasing intermolecular separation.
In Figure 11 and Figure 12, we show the left and right local transition densities of E2 for P, A1, and A2 for 6 Å and 8 Å intermolecular separation, respectively.In this configuration, the transition dipole moment of P is anti-aligned with the rest of the chain.As for the Haggregates, E2 is delocalized over the whole chain with a more relevant contribution from the unperturbed dimers, and the intermolecular forces tend to increase the contribution of P.However, P's local transition dipole moment is anti-aligned with the rest of the hydrogens.Therefore, the alignment pattern of P and A1 in E1 and E2 is swapped compared to the H-aggregate arrangement.This reflects the severe consequences that intermolecular forces can have on the local properties of the system.9. Local left and right transition densities for E1 for an intermolecular separation of 6 Å in a J-aggregate configuration (head-tail arrangement, see Figure 1).The excitation is independent of the chain length N and mainly localized on P, A1, and A2 with decreasing contributions.The local transition dipole moments are aligned, increasing the oscillator strength compared to the isolated P.
Figure 10.Local left and right transition densities for E1 for an intermolecular separation of 8 Å in a J-aggregate configuration (head-tail arrangement, see Figure 1).The excitation is almost independent of N and mainly localized on P, A1, and A2 with decreasing contributions.Compared to Figure 9, the intermolecular forces are less relevant, and therefore, the excitation is more localized on P and less on A1 and A2, but the overall alignment pattern of the molecular transition dipoles is maintained.Figure 11.Local left and right transition densities for E2 for an intermolecular separation of 6 Å in a J-aggregate configuration (head-tail arrangement, see Figure 1).The excitation delocalized over the whole chain with a larger contribution from the unperturbed dimers.The local transition dipole moments are anti-aligned, contrary to E2 for the H-aggregate configuration.
Figure 12.Local left and right transition densities for E2 for an intermolecular separation of 8 Å in a J-aggregate configuration (head-tail arrangement, see Figure 1).The excitation delocalized over the whole chain with a larger contribution from the unperturbed dimers.The local transition dipole moments are anti-aligned, contrary to E2 for the H-aggregate configuration, and compared to Figure 11 there is a smaller contribution from P and a similar contribution from the rest of the chain.

C. Kasha-Frenkel model for H-and J-aggregates
The spectral properties of the H-and J-aggregates can be qualitatively understood using the simple Kasha-Fernkel model [8][9][10][11][12].The interplay between the excitation of different molecules is modeled by a dipole-dipole interaction V AB between their transition dipoles where ϵ is the medium dielectric constant, A and B label the two interacting molecules with isolated transition dipole moments µ A and µ A (usually associated with the S 0 → S 1 excitation), and R AB is the vector, of length R AB , connecting their center of mass.The simple Frenkel exciton Hamiltonian, which describes the first-excitation manifold, is where the states |n⟩ correspond to the n-th molecule excited while the others are in the ground state, and E S0→S1 is the excitation energy of the molecular transition in the supramolecular aggregate, which includes an energy shift from the gas-phase excitation due to the molecular environment [12].Diagonalizing the Hamiltonian H ex for a fixed number N of molecules (or imposing periodic boundary conditions), we can reproduce the blueshift (redshift) of the excitation energy of the Haggregates (J-aggregates).The difference in the energy shift is given by the sign of the intermolecular coupling V AB in the different configurations, which is negative for J-aggregates and positive for H-aggregates.Such a simple model qualitatively accounts for the delocalized nature of the excitation and for the photoluminescence properties of these aggregates [10][11][12], and can be easily generalized including a different molecule (impurity) in the excitation manifold.
While this simple model is straightforward and qualitatively correct, several relevant chemical aspects are disregarded.First, we approximate the excitation manifold, including only the first excited states, while more excited states should be included to obtain a more realistic description.In addition, an a priori knowledge of the "environment shifted" excitation energy E S0→S1 might not be known and could require an ab initio computation, for instance, via polarizable continuum models (PCM) [13][14][15].The effect of molecular vibrations is also fundamental for an appropriate description of the photoluminescence properties of aggregates [12].Moreover, the point-dipole coupling is a crude approximation and significant improvement is achieved by using extended charge distributions [16] or approaches based on quantum chemistry [17][18][19][20][21] which can capture the full com-plexity of the short-range interaction.Indeed, in Figure 13, we show the changes in the ground state density in the XY plane for an (H 2 ) 2 dimer, compared to the isolated molecules.As in the previous calculations, the two dimers have a bond length of 0.76 Å and 0.78 Å and form an H-aggregate.While at large distances, the intermolecular interaction is described by London forces and therefore decays fast with the distance R −6 , at smaller intermolecular distances, the intermolecular coupling has a complex form and depends on the wavefunction overlap.Similar modifications also occur in the excited state densities.In Figure 14, we show the difference between the first excited density of the H 2 H' 2 dimer (mainly associated with the excitation of the perturbed dimer) and the sum of the isolated ground state density of the unperturbed H 2 and the excited state density of the impurity.As can be seen, the differences are more relevant than the ground state changes, and due to intermolecular forces, the excitation of one dimer polarized the other H 2 .Analogous results are seen in Figure 15, where we plot the second excited state density of H 2 H ′ 2 dimer, mainly associated with the excitation of the unperturbed dimer, and subtract the excited density of the unperturbed dimer and the ground state density of the impurity.Such density modifications include rescaling due to collective delocalization and polarization due to intermolecular forces.While for the simple H 2 system the intermolecular forces are small, for different molecules with permanent dipole moments or carrying specific interaction sites such as hydrogen bonds or π − π interactions, the density modifications can be more relevant and even influence chemical processes.These considerations apply to any supramolecular structure, including solute-solvent systems in which the solvation can significantly alter the molecular properties via hydrogen bonds, dipole-dipole, van der Waals, or other specific interactions, in addition to the effects of a different environment dielectric constant ϵ.Indeed, it is well known that solvation can modify chemical reactivity, and a suitable description of both the reactive molecule and the molecular surroundings is necessary to obtain qualitative and quantitative correct predictions.Therefore, a detailed expression of the interactions with the chemical environment is generally necessary in experimental setups such as solid matrices or solutions.
Figure 13.Variations in the ground-state densities of two parallel H2 dimers (with bond length 0.76 Å and 0.78 Å placed on the XY plane, the perturbed dimer is placed at y = 0) with respect to the isolated dimers for different intermolecular distances.While the long-range interaction via the London dispersion force decays fast, there is a significant charge redistribution at small intermolecular distances.At small distances, the modifications are due to both Coulomb coupling and wavefunction overlaps, and therefore, their computation requires a proper quantum description.Similar modifications occur for the excited state densities.Density modifications also occur for larger distances for systems carrying permanent dipole moments.Moreover, specific interactions such as hydrogen bonds can significantly alter molecular densities.In Figure 16, we show the polaritonic (field-dressed) absorption spectra of the (H 2 ) N ≥4 H-aggregate (see Figure 1) for different intermolecular separations.The cavity is tuned to the undressed excitation of the (isolated) unperturbed dimers, and from the interaction between the photon field and the H-aggregate electronic excitations E1 and E2, three polaritonic branches, the lower, middle, and upper polariton emerge.As the intermolecular forces become more relevant (shorter distances), the MP gains oscillator strength while the LP loses it.While the upper and lower polaritonic branches gain intensity with the chain length N , the MP shows a more involved trend.In Figure 17, we show the LP and MP photon characters computed from the transition densities for an intermolecular separation of 5 Å.The photon character of the MP steadily decreases with the chain length N , while the opposite happens for the LP.The trends are consistent with the TC analysis of Figure 2, and analogous trends are found for other intermolecular distances.Therefore, the matter excitations E1 and E2 gain weight in the MP properties as N increases.In Figure 18, we report the energy of the LP and MP as a function of the number of dimers N in the chain for 5 Å intermolecular distance.The avoided crossing between the lower and upper polaritonic branches suggests relevant mixing between the E1 and E2 excitations.Analogous avoided crossings are obtained in the TC model and for different H 2 separations, but they occur at different chain lengths since the intermolecular forces modify the undressed excitations.We can then analyze the local effects of such mixing by extracting the integrated local transition densities, as explained in section I and section II.
In Figure 19 and Figure 20, we plot the LP left and right integrated local transition densities for the Haggregate configuration with an intermolecular distance of 5 Å and 6 Å, respectively.The upper panels of Figure 19 are the same results as 2 in the main text.We see that the responses of A1 and A2 are different because of the intermolecular forces.In the undressed excitation E1, the local transition dipole moments are antialigned, while for E2, they are aligned.In the polaritonic branches, the coupling to the photon mixes E1 and E2 and tends to align the local transition moments as discussed in section III for the TC model.Therefore, a competition between longitudinal and transverse fields results in a change in the alignment pattern of the transition dipole moments in Figure 19, for 5 Å intermolecular distance.For smaller intermolecular forces, at 6 Å distance, the photon effects prevail, and the transition moments are always aligned in the LP.This confirms that the pattern change results from the competition between intermolecular forces and collective strong coupling.Nevertheless, for 6 Å separation, we see significant changes in the response of A1.
In Figure 21 and Figure 22, we plot the MP left and right integrated local transition densities for the H-aggregate configuration with an intermolecular distance of 5 Å and 6 Å, respectively.Compared to the LP branch, the perturbed dimer P transition moment is anti-aligned with the rest of the chain for both distances.
Notice that the MP is mainly localized on P, and the relation between the transition moments of A1 and P resembles E1.Nevertheless, the dimer A2 is now aligned to A1, contrary to E1.This is likely to be due to the mixing of E1 and E2 since, in E2, all of the dimers have a similar contribution in an aligned pattern.In contrast, the contribution of A2 in E1 is tiny and, therefore, easily overcome by E2 and the transverse field effects.Compared to Figure 21, for 6 Å intermolecular separation, the MP is more delocalized with similar contributions from P and A1, while for 5 Å, it is more localized on P. Notice also that for 6 Å, the contribution of P increases with N , while it decreases for A1 and A2 ( Figure 22).Instead, for 5 Å intermolecular separation ( Figure 21), the changes with N are less evident, and the contribution of P decreases with N , behaving similarly to A1 and A2.
The cavity tuning is a fundamental experimental knob, and empirical observations show that detuning the photon disrupts the polaritonic modifications.Therefore, we investigate the effect of the photon frequency on our sys-tem.Notice that the definition of the resonance condition is nontrivial since the excitation E2 blueshifts with increasing N .We then consider four cavity frequencies for 5 Å intermolecular separation: i) tuned to the (undressed) excitation of the isolated perturbed dimer P, ii) tuned to the excitation E1 (basically independent of N ), iii) tuned to E2 for the aggregate (H 2 ) 7 , and iv) completely detuned at higher energies.Resonance with E1 or P provides a model for the photon tuned with the impurity excitation.For the resonance with E2, we expect effects similar to the above results.Finally, when the photon is completely detuned, we expect no modifications compared to the undressed case.In Figure 23, we show the LP local left transition densities of P, A1, and A2 (Haggregate, intermolecular separation of 5 Å).When the photon is tuned to the impurity (either resonant with the undressed excitation of the isolated P molecule or resonant to the excitation E1, first two columns), we observe small local changes that are more pronounced for P. The transition moments are aligned as in E2, and the changes in A1 and A2 are modest.When the cavity is tuned to E2 of (H 2 ) 7 (third column), we see significant changes in the local transition moment of A1, which resemble the results for the photon tuned to the undressed excitation of E1.Finally, when the photon is completely detuned to the excitation (rightmost panels of Figure 23), there is no local change in the transition densities, and the LP shows the same pattern as E1.Since the cavity has a much larger energy than the matter excitation, the hybridization is almost negligible, and the states are effectively the same as the undressed excitations.However, the collective strong coupling for such N is not enough to flip the sign of the A1 transition densities, and longer chains would be necessary.We then recover the resonance condition as a fundamental knob to observe local modifications of molecular properties under cooperative strong coupling.In our model, the cavity frequency determines the efficiency in the hybridization of E1, E2, and the cavity photon, therefore determining the changes in the alignment pattern of the transition dipoles.In Figure 24, we show the LP and MP energies as a function of N for the cavity resonant to A1 (left), E1 (center), and E2 in (H 2 ) 7 (right).In all three cases, the energies display an avoided crossing between the polaritonic branches.Still, when the photon is tuned to the undressed excitation of the isolated A1 (or to E2, first and last columns), the avoided crossing is more pronounced compared to the tuning to E1 (central column).This suggests a more efficient mixing of the undressed excitation and the one-photon state, resulting in the observed local modifications.On the other hand, when the cavity is tuned to E1, the avoided crossing is less pronounced and already overcome for small chain lengths, suggesting a less efficient mixing in the polaritonic branches.This is consistent with the changes in the patterns of the local transition moments.
Finally, we investigate the spectral effect of decreasing the coupling strength while increasing the number of dimers to keep the ratio N/V el fixed, where V el is the  1).As we increase the number of molecules coupled to the photon field, the local transition dipole moment of A1 changes sign.This change in the alignment patterns of the transition moments is a consequence of the competition between the longitudinal intermolecular forces and the collective coupling to the optical mode.This points out the local role of the concentration in polaritonic chemistry and hints at a modification of the solvation response under cooperative strong coupling.
quantization volume (connected to the light-matter couplings strength λ = 1 ε0V ).In Figure 25, we show the absorption spectra (LP and MP branches) for the 5 Å intermolecular separation H-aggregate.The left panel is the same as Figure 16, while the right panel is for a system with twice as much dimers and coupling strength rescaled by 1/ √ 2 to maintain the N/V el ratio.The LP and MP excitation energies are the same in the two realizations, which is consistent with the TC predictions in the thermodynamic limit.At the same time, the peaks are more intense due to the larger number of oscillators contributing to the excitations.1).Compared to Figure 19, the local transition dipole moments are already aligned even for short chains, so there is no sign-flip in the transition moment of A1.Nevertheless, we see that as N increases, the response of A1 becomes more and more relevant and more similar to A2.This behavior is still a consequence of the intermolecular forces between the dimers, but the collective photon effects here prevail since, at larger distances, the longitudinal coupling is smaller.Figure 21.Local left and right transition densities for the middle polariton for an intermolecular separation of 5 Å in an H-aggregate configuration (parallel-bond arrangement, see Figure 1).The P and A1 local transition dipole moments are antialigned as in E1, but, on the opposite, A2 is aligned to A1.The alignment A1-A2 is due to photon field effects from mixing with E2 and the favored alignment of the dimers due to the photon field.Since the contribution of A2 to E1 is small, it is easily overcome by the transverse field effects.The MP is mainly localized on P but also shows significant contributions from A1 and A2. Figure 22.Local left and right transition densities for the middle polariton for an intermolecular separation of 6 Å in an H-aggregate configuration (parallel-bond arrangement, see Figure 1).The P and A1 local transition dipole moments are antialigned as in E1, but, on the opposite, A2 is aligned to A1.The alignment A1-A2 is due to photon field effects from mixing with E2 and the favored alignment of the dimers due to the photon field.Since the contribution of A2 to E1 is small, it is easily overcome by the transverse field effects.The pattern is the same as Figure 21, but the MP is more delocalized over different dimers, with more similar contributions from A1 and P. When the photon is resonant with the impurity (either E1 or the undressed P excitation, first two columns), there are some modifications in the P transition density and smaller changes in A1 and A2.The excitations resemble E2, and no change in the alignment pattern of the transition moments is predicted.When the excitation is resonant with E2 (third column), we observe a trend similar to the cavity tuned to the undressed excitation of A1.Still, the collective strong coupling is insufficient to achieve a flip in the A1 transition dipole, and longer chains are necessary.Finally, when the photon is completely detuned (rightmost column), there is negligible mixing between E1, E2, and the electromagnetic states, and therefore the excitations are analogous to the undressed ones.Therefore, we recover the resonant condition as a fundamental knob to observe local modifications of collective strong coupling since it determines the efficiency in the mixing of the undressed states.
Figure 24.LP and MP energies for three different photon frequencies as a function of the chain length.The avoided crossings between the polaritonic branches indicate a relevant mixing between the molecular states.In resonance with the perturbed excitation (central panel), the avoided crossing is less pronounced and already overcome for small N, suggesting a less efficient interaction.This result is consistent with the less pronounced local transition density modifications.Figure 25.Left, polaritonic absorption spectra for a 5 Å intermolecular separation H-aggregate (H2) N ≥6 (same as Figure 16).Right, spectra for a system with double dimension compared to left, (H2) N ≥12 , and coupling strength rescaled by 1/ √ 2. The LP and MP excitation energies are the same in these two realizations, consistent with the TC predictions in the thermodynamic limit.

A. TC-Kasha model for the molecular polaritons
The sign flip in the local transition density of A1 in Figure 19 can be modeled using a simplified picture that merges the TC and the Frenkel-Kasha models.Following the same procedure as the TC model, we consider the first excited state manifold of the N hydrogen dimers, each with the same excitation energy as the isolated molecules.Each H 2 is then coupled to the photon field via its isolated transition dipole moment as in Equation 34.Following the same assumptions as in the Kasha model in Equation 52, we introduce the intermolecular interactions via the transition dipole-dipole coupling of Equation 51.For simplicity, we couple only the nearest dimers, and the Hamiltonian matrix for the TC-Kasha model then reads Diagonalizing the Hamiltonian in Equation 53, we obtain the eigenvalues and eigenvectors corresponding to 3 polaritonic branches and N -2 (quasi) dark states.Compared to the ab initio simulations, which are limited by the steep computational costs of the CC methods, it is possible to perform simulations with large N .In The effects of rescaling the coupling strength λ for a fixed intermolecular separation by a factor α shifts the avoided crossing from length N roughly to N ′ ∼ N/α 2 , which is consistent to the thermodynamic limit requirement N p /V el = const.for N p , V el → ∞, where V el is the quantization volume λ ∝ V −1 el .Therefore, this local sign flip involving A1, which stems from the competition between longitudinal and transverse fields, is persistent even in the thermodynamic limit.In Figure 27, we show the coefficients of P (left), A1 (center), and A2 (right) for the LP obtained from the matrix Equation 53 where the bare coupling strength is λ 0 = 0.0 a.u., 0.01 a.u., 0.0125 a.u. and 0.015 a.u.further rescaled by 1/ √ N , that is, λ = λ 0 / √ N .In agreement with the ab initio simulations, outside the cavity (first row of Figure 27), the coefficients of P, A1, and A2 display alternate signs due to the intermolecular forces, and the contribution of A1 and A2 to the lowest excitation, mainly associated with P, increases for shorter distances.When the photon coupling is switched on, the coefficient of A1 still displays a sign flip, and the order of magnitude of the coefficients is the same in the three simulations.The position at which such a sign flip occurs depends on λ 0 , the intermolecular separation, the chain length N , and the excitation energy of the impurity.Therefore, the described effects are resistant to the thermodynamic limit for a single impurity P.
The simplified TC-Kasha model successfully reproduces the spectral behavior of the QED-CC results.However, in this framework, the transition and excited state densities can be reconstructed from Equation 53only from the isolated H 2 densities, which are available as an input for the model.Moreover, the ground state is assumed to be unchanged.For the H 2 system under consideration, the modifications of the ground and excited state densities are non-negligible for intermolecular distances below 5 Å, as can be seen from Figure 13, 14, and 15.Therefore, the mixing and rescaling of a few isolated densities coming from the eigenvectors of Equation 53could be unsuited for describing such local modifications, which can be responsible for alterations of chemical processes such as molecular reactivity.On the other hand, QED-CC includes all the quantum chemistry effects due to wavefunction overlap and Coulomb forces.Therefore, when the effect of the solvent is non-perturbative, although a simplified optical model can be used to reproduce the spectral polaritonic properties, an ab initio description of the solute-solvent will be necessary for a detailed investigation of the single-molecule (local) chemistry changes observed in several polaritonic chemistry experiments.While the current discussion has focused on electronic strong coupling, the same considerations hold for vibrational excitations.Quantum chemistry is then necessary to model the potential energy surfaces on which the nuclei move.An analogous interplay between longitudinal and transverse field effects is also expected in the framework of vibrational strong coupling.In particular, solute-solvent interactions, such as hydrogen bonds, will significantly modify the molecular vibrations.Moreover, the light-matter strong coupling can modify intermolecular forces, and it is suggested in particular that van der Waals forces can be particularly affected [22][23][24][25].While a few-level CI model might not capture such electronphoton correlation effects, they are naturally included in ab initio QED.
So far, we have focused on a single perturbed dimer P.However, the TC-Kasha model also allows us to study N rep replicas of the H 2 chain.In Figure 28, we show the TC-Kasha eigenvector coefficients for P, A1, and A2 (from left to right) for coupling strengths λ 0 = 0.004 a.u.(first three panes) and λ 0 = 0.002 a.u.(last three panels) for number of chain replicas N rep = 1, 4 and 16 (from top to bottom for each λ 0 ).The lightmatter coupling strength is then rescaled by 1/ N rep , that is, λ = λ 0 / N rep .For each λ 0 , the results show the same qualitative behavior for every N rep .We still predict a sign flip in the A1 coefficients, which means the qualitative behavior resists the thermodynamic limit.Nevertheless, notice that the coefficients in this case decrease of a factor ( N rep ) −1 and the coefficients will asymptotically approach zero.Similar results are obtained from the ab initio simulations, although the computational costs limit the number of replicas we can study using QED-CC.In Figure 29, we show the LP integrated local transition density matrix for a chain containing 2 perturbed dimers for 5 Å intermolecular distance, with coupling strength rescaled by the factor 1/ √ 2. The two perturbed dimers are separated by 6 unperturbed H 2 so that their direct influence is negligible.For comparison, the transition densities for a single perturbed dimer and unscaled light-matter coupling strength are shown in dotted lines.In agreement with the TC-Kasha model, the qualitative sign-inversion behavior of A1 is reproduced, the position N at which the sign-flip occurs is also the same, and the contribution of each dimer decreases by roughly the same scaling factor as the coupling strength.Nevertheless, we point out that the number of replicas is smaller than the number of coupled dimers, that is, Therefore, the presented results show a slower decay in the thermodynamic limit as 1/N rep instead of 1/N tot .In Figure 30, we plot the coefficients for P, A1, and A2 (left, center, and right panels, respectively) for up to 80 unperturbed dimers (upper panels) for a bare light-matter coupling strength λ of 0.01 a.u. and different intermolecular distances.In the second and third rows of Figure 30, we then duplicate the system, adding 5 (second row) and 10 (last row) identical replicas and consequently rescaling the fundamental coupling strength λ/ N × N rep .The figures show the same qualitative sign flip behavior as a function of N and the intermolecular separation.This confirms again that such a local change, stemming from the interplay between longitudinal (intermolecular forces) and transverse (photon) fields, is resistant to the thermodynamic limit.Figure 19, the coupling strength has been rescaled by ( √ 2) −1 , and the two perturbed H2 are separated by 6 dimers so that their mutual influence is negligible.The dotted lines show for reference the same results as in Figure 19 for a chain with a single perturbed dimer.The results for these two realizations are qualitatively the same and display the same order of magnitude.The differences are likely due to different intermolecular forces arising from the different number of dimers, dipole self-energy, and electron-photon correlation.Using a simplified TC-Kasha model also allows us to study the effect of disorder in the system.To this end, we consider P and its first solvation shell (A1) to be fixed while introducing disorder in the remaining dimers.Therefore, we model only long-range disorder while short-range disorder (in the vicinity of P) is suppressed for simplicity.A Gaussian distribution of the excitation energies, centered at the unperturbed dimers' excitation energies, introduces an energy broadening.Such an energy broadening can be thought of as a distribution of bond lengths in the dimers caused by thermal fluctuations.Moreover, the orientational disorder is introduced by assuming the dimers can form an angle θ k ∈ [0, π/2] with the X axis (k labels the different disordered dimers).At the same time, for simplicity, we keep them orthogonal to the Y axis.The restriction θ k ∈ [0, π/2] (instead of θ k ∈ [0, π]) ensures that the local transition dipole moments are still aligned in the same direction to avoid fictitious flipping of the wavefunction coefficients (see Equation 44).We also assume all the dimers to be kept at a fixed distance.Therefore, the Hamiltonian to be diagonalized here reads In Equation 54, the dimer labeled 1 is P, A1 is labeled 2, and the remaining N-2 dimers have energy disorder embedded in the excitation energy on the diagonal ω(k), and orientational disorder which causes a decrease of the light-matter coupling strength compared to the ordered case g0 (k) = g0 cos(θ k ).Moreover, the orientational disorder also affects the dipole-dipole couplings V (m, n), as seen in Equation 51.Refinement of such simplified disorder can be obtained by relaxing the rotation constraints and permitting a displacement of the H 2 .In Figure 31, we show the LP coefficients obtained from the matrix in Equation 54 for P, A1, and A2 (from left to right) for coupling strengths λ = 0.002 a.u., 0.003 a.u. and 0.004 a.u.
(from the top down to the bottom).These results are analogous to Figure 26.The order of magnitude of the wave function coefficients and the qualitative sign flip of the A1 coefficients are again obtained, although for a larger number of dimers N .Indeed, the collective coupling of the system is decreased because of disorder since the coupling to the cavity, compared to the perfectly ordered case, is scaled by cos θ k for each dimer.Therefore, a larger number of coupled molecules is generally necessary to reach the avoided crossing between the LP and the MP.Nevertheless, the fact that this effect resists the introduction of long-range disorder proves that it stems from the interplay between two distinct molecular environments: the local chemical environment (solvation shell) and the collective optical environment (polaritonic dressing).It also confirms the local nature of such changes.Therefore, we expect the presented results always to be relevant for systems characterized by short-range order, such as solute-solvent systems or aggregates.In Figure 32, we show the dressed absorption spectra for J-aggregate (see Figure 1) for 5 Å and 6 Å intermolecular separation.The cavity frequency is tuned to the undressed excitation of the isolated unperturbed dimers.Compared to the H-aggregates (see Figure 16), the middle polariton branch is less intense and loses oscillator strength with shorter intermolecular separations.This reflects the different effects the local arrangement can have on a solute due to the strong directionality of intermolecular forces.As the intermolecular forces become more relevant (shorter distances), the excitations redshift, in agreement with the predictions from the Kasha model (See section V A).In Figure 33, we plot the LP left and right transition densities for the Jaggregate with an intermolecular distance of 6 Å. Contrary to the H-aggregates, the local transition dipole moments are aligned for every N and show a smaller dependence on the chain length.In Figure 34, we show the corresponding middle polariton.As N increases, the transition dipole moment of A1, initially antialigned with P as in E2, changes sign.As for the H-aggregate, this proves again a local (single-molecule) change stemming from longitudinal and transverse fields.Moreover, the local changes here are more relevant for the middle polaritonic branches (instead of the lower branch as the H-aggregate), which provides a simple proof-of-concept example of how different chemical environments can influence chemical properties due to the directionality of the intermolecular forces.

Figure 1 .Y I − 2 Å
Figure 1.Structure of the (H 2 ) N aggregate and absorption spectra outside the cavity for 5 Å intermolecular separation.The local transition densities ρ j0 (x) (i.e., the transition density ϱ j0 (x, y, z) integrated perpendicularly to the dipole moment direction, that is, perpendicular to the H 2 bond ρ I j0 (x) = +∞ −∞ dz Y I +2 Å Y I −2 Å dy ϱ j0 (x, y, z) where I = P, A1, A2, see Sec.I in the SI) of the undressed excitations E1 (solid lines) and E2 (dotted lines) of (H 2 ) 7 are shown, showing specific transition moment alignment patterns caused by intermolecular interactions.Only the left densities are shown because the right densities display the same behavior (see Sec.I of the SI).

Figure 2 .
Figure 2. Polaritonic spectra for 5 Å intermolecular separation with coupling strength 0.005 a.u., polarization along the H 2 bonds, and photon frequency tuned to the undressed excitation of the isolated unperturbed dimers.The energies of the lower and middle polaritons as a function of the number of dimers N in the chain, depicted in the upper left panel, highlight an avoided crossing.The lower polariton LP integrated local transition densities for P, A1, and A2, are shown.The Coulomb forces compete with the transverse field (photon) to determine the alignment pattern of the transition moments, generating a sign flip in the transition density of A1.

Figure 3 .
Figure 3. Wave function coefficients c I where I = P, A1 from the TC-Kasha model as a function of the environment molecules N −1 and the intermolecular separation.a) The coupling strength is set to λ = 0.015 a.u.√ Ne .We qualitatively recover the same sign flip for A1 as Figure 2 for Nrep = 1 (upper panel) and 4 (lower panel), but the coefficient values approximately scale as 1/ Nrep.Notably, such a decay is much smaller than 1/ √ Ne. b) Wavefunction coefficients obtained by fixing λ and Nrep = 1 as N varies.This is consistent with experimental setups where the same cavity (i.e., λ fixed) is employed at different molecular concentrations (here, the solvent molecules N − 1).c) Same as (b) but for Nrep = 4. Compared to (b1), the overall collective coupling is maintained, but due to the larger number of solute molecules, the coefficients are rescaled by 1/ Nrep.The impurity concentration P then defines the decay of the wavefunction coefficients.Notice that P does not need to be in strong coupling.See Sec.V-A of the SI for more results.

Figure 4 .
Figure 4. LP integrated local transition densities for A1 for a cluster (H 2 ) 3 with intermolecular separation of 3.7 Å, together with N -3 dimers intermolecular at a distance of 20 Å. Due to the larger intermolecular forces between A1 and P, the excitation E1 is strongly blueshifted.Therefore, we use a slightly larger coupling and a larger number N compared to Figure 2 to reach the avoided crossing region between the LP and MP.The results show the same qualitative behavior as Figure 2.

Figure 5 .
Figure 5. Wavefunction coefficients obtained from a disordered TC-Kasha model where the first and second solvation shells around P are fixed while orientational disorder and inhomogeneous broadening are applied to the other dimers.The coupling strength is set to 0.005 a.u..We still recover the same qualitative behavior and coefficient magnitude of Figure 2 and 3(b), although for larger N because of the smaller collective coupling to the device due to disorder.Analogous results are obtained for different coupling strengths, see Sec.V-A of the SI for more results.

Figure 2 .
Figure 2.In the LP (upper panels), the perturbed and unperturbed H2 excitations are in phase (same sign of the CI coefficients) while in antiphase with the photon contribution.The coefficient of the impurity decreases with N , and the photon character becomes more and more relevant, while the unperturbed dimer contribution reaches a maximum around N = 9 and then decreases again.In the MP (lower panels), the perturbed dimers are in antiphase with the photon and the unperturbed H2 molecules.The polariton becomes more localized on the impurity while the photon and the impurity contribution approach zero.

Figure 3 .
Figure 3. Absorption spectra for an H-aggregate (H2) N ≥4 for 5 Å (upper panel), 6 Å (central panel), and 10 Å (lower panel) intermolecular separation.The excitation E2 depends on the number of dimers in the chain, and it blueshifts with stronger intermolecular forces and larger N .The excitation E1 is mainly localized on P and its neighbor dimers and is almost independent of N , while it redshifts with the longitudinal forces.N -2 quasi-dark states gain oscillator strength due to the intermolecular couplings, showing smaller shoulder peaks near E2.

Figure 4 .
Figure 4. Local left and right transition densities for E1 for an intermolecular separation of 5 Å in an H-aggregate configuration (parallel-bond arrangement, see Figure1).The excitation is independent of the chain length N and mainly localized on P, A1, and A2 with decreasing contributions.The local transition dipole moments are anti-aligned, producing an effective screening of the excitation compared to the isolated P.

Figure 5 .
Figure 5. Local left and right transition densities for E1 for an intermolecular separation of 6 Å in an H-aggregate configuration (parallel-bond arrangement, see Figure1).The excitation is almost independent of N and mainly localized on P, A1, and A2 with decreasing contributions.Compared to Figure4, the intermolecular forces are less relevant, and therefore, the excitation is more localized on P and less on A1 and A2.Still, the overall anti-alignment pattern of the molecular transition dipoles is maintained.

Figure 8 .
Figure 8. Absorption spectra for a head-tail aggregate (J-aggregate, see Figure 1) (H2) N ≥4 for 5 Å (upper left panel), 6 Å (lower left panel), 7 Å (upper right panel), and 8 Å (lower right panel) intermolecular separation.The excitation E2 depends on the number of dimers and redshifts with stronger intermolecular forces and larger N .The excitation E1 is mainly localized on P and its neighbor dimers with little dependence N and redshifts with the longitudinal forces.N -2 quasi-dark states gain oscillator strength due to the intermolecular couplings, showing peaks near E2.

Figure
Figure 9. Local left and right transition densities for E1 for an intermolecular separation of 6 Å in a J-aggregate configuration (head-tail arrangement, see Figure1).The excitation is independent of the chain length N and mainly localized on P, A1, and A2 with decreasing contributions.The local transition dipole moments are aligned, increasing the oscillator strength compared to the isolated P.

Figure 14 .
Figure 14.Variations in the first excited state (mainly associated with the impurity) density of two parallel H2 dimers (with bond length 0.76 Å and 0.78 Å placed on the XY plane, the perturbed dimer is placed at y = 0) compared to the isolated ground state of the unperturbed H2 and the excited state of the impurity, for different intermolecular distances.The differences are due to polarization and collective delocalization arising from the intermolecular coupling.

Figure 15 .Figure 16 .
Figure 15.Variations in the second excited state (mainly associated with the unperturbed dimer) density of two parallel H2 dimers (with bond length 0.76 Å and 0.78 Å placed on the XY plane, the perturbed dimer is placed at y = 0) compared to the isolated ground state of the impurity and the excited state of the unperturbed H2, for different intermolecular distances.The differences are due to polarization and collective delocalization arising from the intermolecular coupling.

Figure 17 .
Figure 17.LP and MP photon character as a function of the chain length N for 5 Å intermolecular distance.

Figure 18 .
Figure 18.LP and MP energy as a function of the chain length N for 5 Å intermolecular distance.

Figure 19 .
Figure 19.Local left and right transition densities for the lower polariton for an intermolecular separation of 5 Å in an Haggregate configuration (parallel-bond arrangement, see Figure1).As we increase the number of molecules coupled to the photon field, the local transition dipole moment of A1 changes sign.This change in the alignment patterns of the transition moments is a consequence of the competition between the longitudinal intermolecular forces and the collective coupling to the optical mode.This points out the local role of the concentration in polaritonic chemistry and hints at a modification of the solvation response under cooperative strong coupling.

Figure 20 .
Figure 20.Local left and right transition densities for the lower polariton for an intermolecular separation of 6 Å in an Haggregate configuration (parallel-bond arrangement, see Figure1).Compared to Figure19, the local transition dipole moments are already aligned even for short chains, so there is no sign-flip in the transition moment of A1.Nevertheless, we see that as N increases, the response of A1 becomes more and more relevant and more similar to A2.This behavior is still a consequence of the intermolecular forces between the dimers, but the collective photon effects here prevail since, at larger distances, the longitudinal coupling is smaller.

Figure 23 .
Figure 23.LP integrated local transition densities for P, A1 and A2 for different photon frequencies.When the photon is resonant with the impurity (either E1 or the undressed P excitation, first two columns), there are some modifications in the P transition density and smaller changes in A1 and A2.The excitations resemble E2, and no change in the alignment pattern of the transition moments is predicted.When the excitation is resonant with E2 (third column), we observe a trend similar to the cavity tuned to the undressed excitation of A1.Still, the collective strong coupling is insufficient to achieve a flip in the A1 transition dipole, and longer chains are necessary.Finally, when the photon is completely detuned (rightmost column), there is negligible mixing between E1, E2, and the electromagnetic states, and therefore the excitations are analogous to the undressed ones.Therefore, we recover the resonant condition as a fundamental knob to observe local modifications of collective strong coupling since it determines the efficiency in the mixing of the undressed states.

Fig- ure 26 ,
we show the coefficients of P (left), A1 (center), and A2 (right) for the LP from Equation 53 for coupling strengths λ = 0.002 a.u., 0.003 a.u., 0.004 a.u. and 0.005 a.u. as a function of the intermolecular separation and the chain length N .The qualitative behavior and the order of magnitude of the coefficients are the same in all cases and agree with the ab initio results in Figure19.

Figure 28 .Figure 29 .
Figure 28.LP wave function coefficients c of P, A1, and A2 (left to right).The first three panels refer to a bare coupling strength of λ0 = 0.004 a.u., while for the last three λ0 = 0.002 a.u..For each λ0, the calculations are repeated for 1, 4, and 16 duplicates of the chain (from top to bottom), and the coupling strength is rescaled by 1/ Nrep i.e. λ = λ 0 √ Nrep .

Figure 30 .
Figure 30.Lp wave function coefficients c of P, A1, and A2 (left to right) for 1, 6, and 11 duplicates (from the top down to the bottom) of the H-aggregate system with up to 40 bare molecules.The cavity strength is scaled as λ = 0.01/ √ Ntot to investigate the thermodynamic limit N, V el → ∞, N/V el = c.

Figure 31 .Figure 32 .
Figure 31.LP wave function coefficients c of P, A1, and A2 (left to right) with coupling strengths λ = 0.002 a.u., 0.003 a.u. and 0.004 a.u.(from the top down to the bottom) of the H-aggregate system with up to 250 bare molecules.Disorder is introduced in the system via a Gaussian energy broadening of the dimers (except for P and A1), with Gaussian width 0.02 eV.Moreover, the same dimers are given a random angle with the X axis (while keeping them perpendicular to Y) θ k ∈ [0, π/2].

Figure 33 .
Figure 33.Local left and right transition densities for the lower polariton for an intermolecular separation of 6 Å in a J-aggregate configuration head-tail arrangement, see Figure 1).The local transition dipole moments are all aligned for all N .

Figure 34 .
Figure 34.Local left and right transition densities for the middle polariton for an intermolecular separation of 6 Å in a Jaggregate configuration head-tail arrangement, see Figure 1).The local transition dipole moment of A1 changes sign with N .