Frenkel excitons are the primary photoexcitations in organic semiconductors and are ultimately responsible for the optical properties of such materials. They are also predicted to form *bound* exciton pairs, termed biexcitons, which are consequential intermediates in a wide range of photophysical processes. Generally, we think of bound states as arising from an attractive interaction. However, here, we report on our recent theoretical analysis, predicting the formation of stable biexciton states in a conjugated polymer material arising from both attractive and repulsive interactions. We show that in J-aggregate systems, 2J-biexcitons can arise from repulsive dipolar interactions with energies *E*_{2J} > 2*E*_{J}, while in H-aggregates, 2H-biexciton states with energies *E*_{2H} < 2*E*_{H} can arise corresponding to attractive dipole exciton/exciton interactions. These predictions are corroborated by using ultrafast double-quantum coherence spectroscopy on a [poly(2,5-bis(3-hexadecylthiophene-2-yl)thieno[3,2-b]thiophene)] material that exhibits both J- and H-like excitonic behavior.

## I. INTRODUCTION

It is generally understood that the primary photoexcitations in organic semiconducting materials are molecular *π*–*π** electronic singlet states (S_{1}) termed Frenkel excitons.^{1} While local in nature, at sufficiently high packing densities, excitons can be delocalized over several molecular units and sufficiently higher excitation densities, exciton–exciton interactions begin to dominate the optical properties of such materials.^{2} Biexcitons, bound pairs of excitons, are consequential intermediates in a wide range of photophysical processes, such as exciton dissociation into electrons (e^{−}) and holes (h^{+}),^{3}

bimolecular annihilation,^{4}

and singlet fission producing triplet (T_{1}) states,^{5}

In Ref. 4, we noted that bimolecular annihilation may be mediated both by resonance energy transfer and by diffusion-limited exciton–exciton scattering, but in the either case, we invoke the key intermediate [2S_{1}]^{‡}. Examples of this occur in biological light harvesting complexes where multi-exciton interactions may play important roles^{6} in the excitonic transport process, and biexcitons can be crucial in cascade quantum emitters as a source of entangled photons.^{7} While ample theoretical work points toward the existence of biexcitons in organic solids^{8–14} and in optical lattices,^{15} there has been only indirect evidence of the dynamic formation of two-quantum exciton bound states in polymeric semiconductors by incoherent, sequential ultrafast excitation.^{3–5,16,17}

Recently, we reported on the *direct* spectroscopic observation of *bound* Frenkel biexcitons, i.e., bound two-exciton quasiparticles ([2S_{1}]^{‡}), in a model polymeric semiconductor [poly(2,5-bis(3-hexadecylthiophene-2-yl)thieno[3,2-b]thiophene)]^{18} (PBTTT) using coherent two-dimensional ultrafast spectroscopy.^{19} The chemical structure of PBTTT is given in Fig. 1(a). PBTTT is unique in that depending on processing conditions, it can support the formation of both H and J aggregate single exciton states, suggesting an arrangement as sketched in Fig. 1(b) in which intra-chain J-like excitons can form along the chains spanning over several PBTTT subunits, while inter-chain H-like excitons can form due to the parallel stacking of several chains within the aggregate. The experimental observations revealed a correlation between peaks in the single and double quantum spectra that correspond to the formation of 2H and 2J biexciton species. This conclusion was supported by both a computational model and theoretical analysis based on a quasi-one-dimensional continuum model.

Here, we present an overview of the theoretical model we developed for biexcitons and use it to discuss biexcitons in related organic polymer materials. First, we show how one can reduce the two-dimensional lattice problem into two separate one-dimensional problems and use Green’s function approach to account for the contact interaction between excitons. This gives the criterion for the overall stability of the bound biexciton states in terms of the exciton bandwidth and contact interaction. Formally, the lattice model reduces to a one-dimensional continuum model with *δ*-function potential. We append to this “text-book” model a deformable classical medium to examine the contribution of the lattice reorganization about the bound-biexciton states and find that this stabilizes the attractively bound biexciton state but destabilizes the repulsively bound states. We conclude by discussing the experimental observations and the need for a better theoretical understanding of bound biexciton states.

## II. THEORETICAL MODEL

### A. Homogeneous lattice model in one- and two-dimensions

To explore the possibility of having multiple species of bound, biexcitons in the same system, we begin by writing a generic lattice model for the system by defining exciton operators *a*_{n} and $an\u2020$,

These operators are Paulion operators that create and remove single excitations on a site labeled **n**. On a given site, they obey the Fermion relation ${an,an\u2020}=1$, which enforces Pauli exclusion but commutes “off-site” with $[an,am\u2020]=0$ when **n** ≠ **m**. This is different from the usual Fermion algebra where the anti-commutation rule applies over all sites, ultimately giving rise to the exchange interaction. Furthermore, we can write a generic multi-exciton Hamiltonian as

where *h*_{nm} describes the single-exciton dynamics and *U*_{nm} represents the exciton/exciton interaction. In principle, the parameters entering into the Hamiltonian in Eq. (5) are defined by the system of interest. For the case of excitons, the diagonal elements of the single particle term define the energy to place an exciton in site **n**, and we write *h*_{nn} = *ϵ*_{n}. For a homogeneous lattice, all site energies are the same and *ϵ*_{n} = *ϵ*_{o}. Similarly, the off-diagonal elements of *h*_{nm} correspond to the matrix elements for transferring an excitation from site **n** → **m**. To a good approximation, the single-exciton transfer interaction can be described within the dipole–dipole approximation as described above. This model differs from the Hubbard model commonly studied in condensed matter physics in that we explicitly exclude double occupancy of each lattice and the exciton/exciton interaction is taken to be between occupied neighbors. Formally, a Frenkel exciton corresponds to a single electron/hole excitation on a given site. However, molecules are not point particles and excitons may acquire some intramolecular charge-transfer character. Therefore, we anticipate that *U*_{nm} is also dipole–dipole-like and reflects the relative orientation of the *static* exciton dipole moments.

For a one-dimensional chain with lattice spacing *a*, *n* is simply an index along the chain such that the site location is given by *r*_{n} = **n***a*. However, for two- and three-dimensional systems, we shall take it as an n-tuple index specifying the site location. For the single particle term, *h*_{nn} is the excitation energy for single site (*ɛ*_{n}) and *h*_{nm} (for **n** ≠ **m**) corresponds to the hopping integral between sites. Upon transforming into the reciprocal space,

one finds the single particle energy dispersion as

To determine the two-exciton states, we begin by writing

where *c*_{kk′} are the expansion coefficients for this state. At this point, there are various approaches one can take to find the general solutions for the Schrödinger equation for the two-exciton system. Indeed, for a small enough lattice, one can simply directly diagonalize the Hamiltonian in Eq. (5) for a finite sized grid. However, we are not interested in the full solution of this problem. Rather, we are only focused on solutions corresponding to *bound* exciton pairs, especially those bound pairs that retain their J- or H-like excitonic character. With this in mind, we develop an analytical solution that naturally extends to the full model.

### B. Local exciton approximation

In Ref. 19, both direct diagonalization and a lattice Greens function approach developed by Vektaris^{9} were employed to study the properties of the biexciton, including its dispersion and the effects of local disorder. A key assumption in our model is that we can define two equivalent quasi-one-dimensional representations for H-like or J-like excitons. For this, let us define a new set of exciton operators, $J\u0302kx\u2020(n)$ and $J\u0302kx(n)$, which creates or removes an exciton with wavevector *k*_{x} in the *x*-direction localized on the *n*th row of sites. Similarly, we define operators $H\u0302ky\u2020(m)$ and $H\u0302ky(m)$, which create and remove excitons with wavevector *k*_{y} in the *y*-direction, but are localized to the m-th column, as sketched in Fig. 2. These can be written in terms of the original lattice operators

Both are equivalent representations, and we can choose to use either (but not both) to rewrite the original problem in this new representation.

Thus, we can write

where in the first line, we diagonalized in the J-direction, and in the second, we diagonalized in the H-direction. This implies that we can think of a J-exciton state as moving in the H-direction with hopping integral *t*_{H} and H-exciton states as moving in the J-direction with hopping integrals *t*_{J}. The dispersion relations are, then, as usual given as

where *t*_{J} < 0 is the nearest-neighbor coupling along the *x*-direction and *t*_{H} > 0 is the nearest-neighbor coupling taken in the *y*-direction. To remind, throughout, we are taking the (dimensionless) wavevector *k* ∈ [−*π*, *π*]. Since the problem is formally separable between *x* and *y* directions, the single-particle terms do not mix wavevector components since *k*_{x} and *k*_{y} are “good” quantum numbers for this system. Using the $J\u0302$ and $H\u0302$ operators, we can reduce the two-particle/two-dimensional problem into one for a pair of particles within a one-dimensional frame without loss of generality. In both cases, the allowed optical transition occurs at *k*_{x} = *k*_{y} = 0.

This decomposition suggests that motion in one direction can be very different from the motion in the perpendicular direction due to the orientations of the transition dipoles between neighboring units. In Fig. 1(b), we suggest how this can be accomplished within the context of molecular aggregates with *π*-stacking. Here, excitonic sites are denoted as blue ellipsoids along with their respective local transition dipole moments (yellow arrows). Along the *x*-direction of the two-dimensional lattice sketched here, the transition moments are oriented more or less in a “head-to-tail” arrangement, producing a hopping matrix element in *h*_{nm} < 0. In the *y*-direction, however, the transition moments are aligned co-facially, producing a single-particle hopping matrix element *h*_{nm} > 0. In the former case, the optically bright state occurs at the bottom of the energy band (J-aggregate), while in the latter case, the optically bright state occurs at the top of the energy band (H-aggregate). The PBTTT material is unique in that both J- and H-aggregate states can be readily observed depending on the sample preparation.

If we use the *J*_{km} (or *H*_{kn}) states as a basis for a given value of *k*, we can use Green’s function approach to find the biexciton energies as

where *ϵ*_{x} is the excitation energy of either the *J* or *H* single exciton and *U* is the contact interaction.^{9} This expression is predicated on |*U*/2*t*| > 1 in order for the biexciton wavefunction to decay exponentially with exciton–exciton separation. These dispersions are plotted in Fig. 3(a) for both 2H and 2J biexcitons. In each case, the energy origin should be shifted to twice the J or H exciton energy. At *k* = 0, the difference between the interacting and localized excitons is the contact energy *U*, which defines the binding energy for the exciton pair. The binding energy must be at least greater in magnitude than 4*t*; else, the lowest energy interacting state will be still within the band for the freely dissociated pairs. These bands will split from the freely dissociated bands once *U*/2*t* > 1.

### C. Biexciton stability

According to our model, Frenkel biexcitons mix J-like and H-like character in terms of their collective quantum behavior with the requirement that the ratio of the exciton/exciton interaction and the perpendicular hopping term be *U*/*t* > 0, which gives rise to localized biexciton states in the perpendicular direction. For the 1D *δ*-function potential, any attractive interaction with *U*/*t* > 0 produces a localized state with localization length *λ* = *κ*^{−1} = 2*t*/*U*.^{20} For the 1D lattice, bound biexciton states occur outside the band for free biexcitons. To gain further insights into the stability of these states, we turn to a continuum model and work in a relative coordinate reference frame where *x* = |*r*_{1} − *r*_{2}| is the separation between two localized excitons. Thus, the biexciton Schrödinger equation can be approximated as

where *U* is the contact interaction between the two excitons. For bound states, *ψ*(*x*) must vanish as *x* → ±*∞*, giving that

where *κ* = *U*/2*t* is a positive constant and *E* = *tκ*^{2}. In general, we take *t* = −*ℏ*^{2}/2*μ*_{eff} and *U* < 0 for an attractive potential, giving rise to a bound state energetically *below* the continuum for the scattering states.

#### 1. Lattice reorganization in the impurity model

Generally speaking, one cannot discount the role of lattice reorganization and relaxation when discussing excitons and biexcitons in organic polymer semiconducting systems. To study this, we append to the 1D impurity model a continuum model for the medium, a term coupling the biexciton to the lattice as per the Davydov model.^{21–30} The resulting equations of motion read as

where *u*(*x*) is the lattice deformation, *χ* describes the linear coupling between the biexciton and the lattice, *m* is the mass of the lattice “atoms,” and *k* is the elastic modulus. If we seek traveling wave solutions, *u*(*x*, *t*) = *u*(*x* − *vt*), where *v* is the group velocity, we find a closure relation

that gives us a non-linear Schrödinger equation

where *g* = −4*χ*^{2}/(*k*(1 − (*v*/*c*)^{2})) and *c* is the sound velocity. Note that *E*_{o} is a constant, given by

which we can ignore for purposes of this analysis. The *δ*-function potential implies that the wavefunction should have the form in Eq. (16). Taking *κ* as a variational parameter and minimizing the total energy, one obtains

Here, *κ* > 0 is necessary to produce a localized state, and from the above equation, *U*/*t* > 0, and *g* < 0 from its definition; we have a stability requirement that if *U* > 0 and *t* > 0, then −*g* < 4*U*. Solving for the binding energy,

we obtain a straightforward estimate of the contribution of both the lattice and the exciton/exciton coupling to the biexciton binding.

In Fig. 4, we plot the biexciton binding energy vs the non-linearity parameter, *g*. For the attractively bound 2H case, lattice reorganization is expected to stabilize the biexciton state by further localizing the state (*κ* increases as *g* increases in magnitude). On the other hand, for the 2J state, increasing the magnitude of *g* decreases *κ* and destabilizes the otherwise bound 2J state. When −*g* = 2*U*, the state is fully delocalized and further increases in the lattice coupling lead to unbound solutions.

## III. COMPARISON TO EXPERIMENTAL MEASUREMENT

We have examined the theoretical concepts by means of two-dimensional coherent excitation spectroscopy on PBTTT, with the structure depicted in Fig. 1(a), which we reported extensively in Ref. 19. In that work, we identified spectral features associated with the 0–0 excitation peak of both the H- and J-aggregate components of the hybrid aggregate, with cross peaks reflecting spectral correlations due to their shared ground state. The origin of the H-aggregate vibronic progression was centered at 2.06 eV, while a weaker peak at 1.99 eV was assigned to the J-aggregate vibronic origin. By performing two-quantum coherence measurements, we found spectral signatures of both 2H and 2J biexcitons. A cross section of the 2D spectral data along the two-quantum energy axis (*ℏω*_{2Q} − *ℏω*_{diag}) relative to the two-quantum diagonal (*ω*_{diag} ≡ *ω*_{2Q} = 2*ω*_{1Q}) is shown in Fig. 5. We found that 2H biexcitons displayed *attractive* biexciton binding with energy −64 ± 6 meV, whereas 2J biexcitons displayed *repulsive* correlations with binding energy +106 ± 6 meV, which is that the energy of the 2H-biexciton resonance is lower than twice the H-aggregate resonance energy, while the corresponding energy for the 2J-biexciton resonance is higher.

We rationalized this observation, as depicted in Fig. 1: two quantum interactions for excitons dispersed along the polymer backbone (J aggregates) are with *J* < 0, while *J* > 0 for those between excitons dispersed on several chains (H aggregates). Considering physically reasonable parameters, we concluded that biexcitons in PBTTT are stable by the arguments depicted in Fig. 5.

## IV. PERSPECTIVE

We presented here theoretical and experimental evidence supporting the formation of bound Frenkel biexcitons in a molecular aggregate material. In our theoretical analysis, we solved the full 2D interacting model and gave the conditions necessary for the formation of stable, stationary states corresponding to bound exciton pairs. This model provides a road-map for developing a bi-exciton material genome in terms of the properties of the J- and H-excitons. We showed that for bound biexcitons, both the exciton–exciton interaction *U* and the single-particle hopping integral *t* must have the same sign. Furthermore, *U*/*t* > 2 so that the exciton/exciton potential interaction is greater than their total kinetic energy. Curiously, we found that while H-like excitons form bound biexcitons with attractive interactions, J-like exciton pairs form bound states arising from a repulsive interaction. Curiously, we found that while the 2H biexciton is stabilized by interactions with the lattice phonons, the 2J biexciton is destabilized to the extent that strong lattice/exciton interactions will only produce unbound biexciton states, which certainly complicates the observation of these states in systems with strong exciton/phonon coupling. Nonetheless, the 2J biexciton is clearly apparent in the 2D two-quantum experiments in Ref. 19.

An open question, however, is the nature of the exciton–exciton interaction itself. Here, we introduced it as a parameter into our model with a dipole–dipole-like form in that the exciton–exciton coupling in one direction is different from that in the perpendicular direction. It is also important to point out that while this interaction seems to have a dipole-like form, it necessarily must reflect the *static* dipole of the local Frenkel excitons rather than their transition dipole moments, which are responsible for the single exciton transfer between local sites. Computing these interactions from a first-principle *ab initio* theory remains a formidable challenge since it necessitates the accurate calculation of doubly excited states with some degree of charge-separation.

## ACKNOWLEDGMENTS

The work at the University of Houston was funded, in part, by the National Science Foundation (Grant No. CHE-2102506) and the Robert A. Welch Foundation (Grant No. E-1337). The work at Georgia Tech was funded by the National Science Foundation (Grant No. DMR-1904293).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

## REFERENCES

Note, that here *t* = −ℏ^{2}/2*μ* carries units of [Energy] · [length]^{2} and *U* carries units of [Energy] · [length]^{−1} so that *λ* carries the appropriate units of [length].