# Anomalous Phonon in Charge-Density-Wave Phase of Kagome Metal CsV3Sb5

Han-Yu Wang  Key Laboratory of Materials Physics,Institute of Solid State
Physics, HFIPS, Chinese Academy of Sciences,Hefei 230031, China Science Island
Branch of Graduate School, University of Science and Technology of China,
Hefei 230026, China Xiao-Cheng Bai  Key Laboratory of Materials
Physics,Institute of Solid State Physics, HFIPS, Chinese Academy of
Sciences,Hefei 230031, China Science Island Branch of Graduate School,
University of Science and Technology of China, Hefei 230026, China Wen-Feng Wu
Key Laboratory of Materials Physics,Institute of Solid State Physics, HFIPS,
Chinese Academy of Sciences,Hefei 230031, China Science Island Branch of
Graduate School, University of Science and Technology of China, Hefei 230026,
China Zhi Zeng  Key Laboratory of Materials Physics,Institute of Solid State
Physics, HFIPS, Chinese Academy of Sciences,Hefei 230031, China Science Island
Branch of Graduate School, University of Science and Technology of China,
Hefei 230026, China Da-Yong Liu  Department of Physics, School of Sciences,
Nantong University, Nantong 226019, China Liang-Jian Zou
<EMAIL_ADDRESS>Key Laboratory of Materials Physics,Institute of Solid
State Physics, HFIPS, Chinese Academy of Sciences,Hefei 230031, China Science
Island Branch of Graduate School, University of Science and Technology of
China, Hefei 230026, China

###### Abstract

CsV3Sb5, a notable compound within the kagome family, is renowned for its
topological and superconducting properties, as well as its detection of local
magnetic field and anomalous Hall effect in experiments. However, the origin
of this local magnetic field is still veiled. In this study, we employ the
first-principles calculations to investigate the atomic vibration in both the
pristine and the charge-density-wave phases of CsV3Sb5. Our analysis reveals
the presence of “anomalous phonons” in these structures, these phonon induce
the circular vibration of atoms, contributing to the phonon magnetic moments
and subsequently to the observed the local magnetic fields. Additionally, we
observe that lattice distortion in the charge-density-wave phase amplifies
these circular vibrations, resulting in a stronger local magnetic field,
particularly from the vanadium atoms. This investigation not only reveals the
potential relation between lattice distortion and atomic polarization but also
offers a novel idea to understand the origin of local magnetic moment in
CsV3Sb5.

††preprint: APS/123-QED

_Introduction.-_ The kagome lattice, which consists of planar triangles and
hexagons structures, was first introduced to statistical physics by Syozi [1].
In past decades, many unusual properties and behaviours, such as geometrical
frustration [2], quantum spin liquid states [3], charge-density-wave (CDW) [4,
5], topological insulator [5], and unconventional superconductivity (SC) [6],
etc., were proposed for this lattice, suggesting kagome compounds may display
numerous quantum phases. Very recently, some new kagome metals AV3Sb5 (A=K,
Rb, and Cs) have been discovered [7]. Universal chiral CDW order coexisting
with unconventional SC at low temperatures was observed [8, 9, 10] in these
compounds. The CDW phases show $2\times 2$ charge order, and the compounds
display three SC phases with increasing hydrostatic pressure; these kagome
metals demonstrate $Z_{2}$ topological surface state was observed [11], and so
on. These compounds integrate many important concepts in condensed matter
physics, raising great attention and debates on the microscopic origins of the
CDW order and SC order.

Among these unusual properties, two phenomena are especially attractive. On
the one hand, both the neutron scattering experiment [12] and the muon spin
relaxation ($\mu$SR) experiment [13, 14] demonstrated the absence of
measurable local spin moments or magnetic correlations, implying a perfect
nonmagnetic phase in kagome metals AV3Sb5. However, a giant anomalous Hall
effect (AHE) was observed, concomitant with the charge order at low
temperatures [15, 16, 17], indicating the presence of intrinsic magnetic
field. On the other hand, the time-reversal symmetry breaking (TRSB) also
occurs with the chiral CDW ordering in the scanning tunneling microscopy
experiment [17, 18, 19]. Whereafter, chiral orbital flux or loop current was
put forward by Feng et al. [20] and Denner et al. [21] for addressing the CDW
ordering and TRSB; on the contrary, several experiments denied the presence of
the chiral flux current in bulk AV3Sb5 [22, 23]. The topical issue thus arises
what causes the AHE and the TRSB in the CDW phase of these nonmagnetic metals.

It is widely recognized that when charged ions undergo circularly or
elliptically vibrations, they generate extra orbital angular momenta and extra
orbital magnetic moments. Recently, Zhang et al. proposed that coherent
motions of two components of charged ions could result in chirality of phonon
[24]. This specialized phonon, known as chiral phonon, emerge from the global
spatial inversion or the time-reversal symmetry breaking within specific
system [24, 25], including hexagonal lattice [24, 26, 27] and kagome lattice
[28]. The circular or elliptical vibrations of chiral phonons thus may
contribute extra orbital angular momenta and orbital magnetic moment, leading
to the phonon magnetic moment [29, 30, 31]. Conversely, in certain
conventional materials preserving the spatial inversion and time-reversal
symmetries, charged ions perhaps still exhibit circular or elliptical
vibration in particular phonons. However, due to the presence of opposite
polarizations of atomic vibrations, the total polarizations of these phonons
are zero, thus they do not have chirality [24, 28, 26, 27]. To distinguish
this kind of phonon from the conventional one, here, we term it “anomalous
phonon”, which has local orbital magnetic moment like chiral phonon. In this
context, the relationship between the chiral phonon and the anomalous one
resembles that between ferromagnet and antiferromagnet.

In this Letter, we investigate the atomic polarization in the pristine and the
CDW phases of CsV3Sb5, utilizing an extension from chiral phonon to anomalous
phonon for topologically kagome metals. We reveal the existence of anomalous
phonons in these systems and show the vanadium atoms contribute extra phonon
magnetic moment in the CDW phase. Such an orbital magnetic moment provide a
novel idea to explain the emergence of measurable local magnetic field [13,
14], AHE [15, 16, 17] and TRSB [17, 18, 19] in low-T phase of AV3Sb5.

_Crystal structures.-_ To investigate the phonon of CsV3Sb5, we initially
optimize the structures by the first-principles calculations. Fig.1 (a) and
(b) show the optimized pristine phase of CsV3Sb5. The structure showcases a
layered arrangement of V-Sb atomic sheets interspersed with Cs atomic layers,
and the V atoms construct the kagome plane in $ab$-plane. The crystal
structure is classified to the $P6/mmm$ space group with $C_{6}$ symmetry. Cs
and V atoms each have an equivalent site, while Sb atoms have two inequivalent
sites, which are labeled as Sb1 and Sb2 in Fig.1 (a) and (b). The lattice
constants of the $ab$-plane and the $c$-axis are 5.45 Å($a$,$b$) and 9.30
Å($c$), respectively, in agreement with experimental observation at high
temperature [19].

Figure 1:  Crystal structures of CsV3Sb5. (a) The oblique view of the pristine
phase. The red and blue layers indicate the kagome plane and Cs plane,
respectively. (b) The top view of the pristine phase with the $2\times 2$ unit
cells. (c) The top view of the ISD $2\times 2$ CDW phase. The red arrows mark
the distortion direction of V atoms in $ab$-plane.

As the temperature decreases to the critical value, CsV3Sb5 transits from
high-temperature homogenous pristine phase to the low-temperature CDW phase.
In this CDW phase, the V atoms undergo a shrinkage resembling an inverse Star
of David (ISD) phase with the $2\times 2$ supercell [32, 8], as illustrated in
Fig.1 (c). Upon transitioning to the ISD $2\times 2$ CDW phase, the six V
atoms in hexagonal configuration move toward to the center, while the three V
atoms in the triangle configuration shift outward, corresponds to a breathing-
in phonon within the kagome plane. This optimized crystal structure of the CDW
phase agrees with the prior findings [32, 8], maintaining the $P6/mmm$ space
group with $C_{6}$ symmetry. However, due to such lattice distortion, there
emerge two inequivalent sites (V1 and V2) for the V atom, as well as two
inequivalent sites (Sb’1 and Sb’2) for Sb2 atom, which are labeled in Fig.1
(c).

_Anomalous phonon.-_ Chen et al. [28] reported the presence of chiral and
anomalous phonons in the kagome lattice. To ascertain whether anomalous
phonons also exist in CsV3Sb5, we investigate the phonon dispersion and atomic
polarization $S^{z}_{\alpha}$ of both the pristine and CDW phases within the
first Brillouin zone (BZ), the magnitude of $S^{z}_{\alpha}$ indicates the
extent of atomic circular vibration. For detailed calculations and discussions
on the polarization $S^{z}_{\alpha}$, please refer to the Supplemental
Materials Sec.S3.

Fig.2 shows the phonon dispersions with atomic polarization $S^{z}_{\alpha}$
distribution of the pristine phase in the high-symmetry path of the BZ. The
unit cell comprises eight atoms, resulting in 24 phonon branches in the
dispersion. As shown in Fig.2 (a), all Cs, V and Sb1 atoms exhibit zero
polarization $S^{z}_{\alpha}$, indicating the absence of circular vibration
for these atoms. Conversely, Fig.2 (b) illustrates the polarization of an Sb2
atom, this atom displays nonzero polarization $S^{z}_{\alpha}$ in certain
phonons, signifying its engagement in circular vibration and the presence of
anomalous phonons in the pristine phase. These anomalous phonons predominantly
occur at low frequencies ranging from 1 THz to 4 THz. It is noteworthy that
Cs, V and Sb1 atoms reside within the symmetric $ab$-plane, representing
spatial inversion symmetry points. Therefore, the local atomic circular
vibrations in these phonons may originate from the inherent asymmetry of these
sites.

Figure 2:  Phonon dispersions of the pristine phase with atomic polarization
$S^{z}_{\alpha}$ along the high-symmetry path. (a) Phonon dispersion with the
atomic polarization of Cs, V and Sb1 atoms. (b) Phonon dispersion with the
atomic polarization of Sb2 atom. The red and blue parts represent the opposite
direction of circular vibration.

In the CDW phase, the Cs atoms occupy two inequivalent sites and the atomic
polarizations of these atoms remain zero. Similarly, the sites of Sb1 atoms,
which were equivalent in the pristine phase, evolute into two inequivalent
sites in the CDW phase, with both exhibiting zero atomic polarization.
Moreover, the sites of Sb2 atoms in the pristine phase also evolute to two
inequivalent sites (Sb’1 and Sb’2) in the CDW phase, and the nonzero
polarization distribution resembles that of the pristine phase. Distinctively,
the two inequivalent V1 and V2 atoms also display circular vibration in
certain phonon modes. Fig.3 (a) and (b) show the phonon dispersions and atomic
polarizations of V atoms in the CDW phase. Due to the distortion in the
2$\times$2 supercell, this dispersion comprises 96 phonon branches, with those
involving circular vibrations of V atoms exhibiting higher frequencies
compared to those involving circular vibrations of Sb atoms. Compared to the
pristine phase, V atoms experience the $ab$-plane distortion in the CDW phase.
Consequently, the emergence of V atomic polarization stems from the distortion
of the kagome plane.

Figure 3:  Phonon dispersions of the CDW phase with atomic polarization
$S^{z}_{\alpha}$ along the high-symmetry path. (a) Phonon dispersion with the
atomic polarization of V1 atom. (b) Phonon dispersion with the atomic
polarization of V2 atom. The red and blue parts represent the opposite
direction of circular vibration.

These nonzero atomic polarizations indicate the presence of anomalous phonons
in both the pristine and CDW phases of CsV3Sb5. However, according to our
analysis, the total atomic polarizations of these anomalous phonon is zero,
thus these anomalous phonons are not chiral phonons. We will delve into this
discussion in detail in the subsequent section.

_Atomic polarization and circular vibration.-_ Atomic polarization serves as a
measure of the magnitude of circular vibration. To investigate the
polarization values and the nature of atomic vibration, we numerically label
these phonon modes based on their frequency magnitude ranging from low to high
and focus on the high-symmetry $K$ point as an example, to elucidate the
anomalous phonons. Table.1 lists the nonzero polarizations of atoms in select
anomalous phonons. Notably, Mode 8 in the pristine phase and Mode 27 in the
CDW phase exhibit similar frequencies, suggesting an evolutionary relationship
between them. Besides, Sb’2 atoms evolve from Sb2 atoms following structural
distortion. In these phonon modes, both Sb2 and Sb’2 atoms exhibit have
nonzero atomic polarization. Considering that the phonon eigenvector is
normalized, it is appropriate to compare polarizations after expanding the
unit cell of the pristine phase to a $2\times 2$ supercell. Please refer to
Supplemental Materials Sec.S3 for discussions on this matter. Following this
expansion, the polarization value of each atom will be 1/4, i.e., 0.0625,
which is half of the polarization value of Sb’2 atoms in the CDW phase. This
indicates that the Sb’2 atoms in the CDW phase exhibit a greater degree of
circular polarization compared to their counterparts. Regarding the V2 atoms
in the CDW phase, although their polarization is relatively smaller compared
to the Sb’2 atoms, they still exhibit significant polarization properties,
surpassing the polarization of Sb2 atoms in the pristine phase. The
corresponding energy of this anomalous phonon is higher, approximately 5.90
THz.

Table 1:  The magnitude of atomic polarization in specific anomalous phonon modes at the high-symmetry $K$ point. The sites of Sb2, Sb’2, and V2 atoms are labeled in Fig.1 and 4. Phase |  Pristine Phase |  CDW Phase
---|---|---
Atom | Sb2 |  Sb’2 | V2
Mode |  8 | 16 |  27 | 61 | 77
Frequency (THz) | 2.27 | 3.72 | 2.35 | 3.98 | 5.90
Polarization ($\hbar$) | 0.250 | 0.250 | 0.127 | 0.134 | 0.086
Figure 4:  Atomic vibrations in specific anomalous phonons at the high-
symmetry $K$ point. (a) The atomic vibrations of Mode 8 and 14 in the pristine
phase. (b), (c), (d) The atomic vibration of Mode 27, 61 and 77 in the CDW
phase, respectively. The radius of the circular trajectory around atoms are
proportional to the amplitude of atomic circular vibration. The red and blue
parts represent the opposite directions of atomic vibration.

To facilitate a more intuitive discussion, we depict the vibration trajectory
of each atom in these phonon modes, as shown in Fig.4. The trajectories of
these atoms maintain $C_{6}$ symmetry within these crystal structures,
consistent with the symmetry of the crystal lattice. In these phonon modes,
for Sb2 and Sb’2 atoms, the vibration trajectories form complete circles
within the $ab$-plane. Additionally, the circularly polarized Sb atoms are
paired in a unit cell, with their vibration directions opposing each other.
Thus, the total phonon polarization is zero, and there is no chirality present
in this phonon mode, owing to the spatial inversion symmetry inherent in the
system. However, as discussed above, for the other Sb atoms located at spatial
inversion symmetry points, they do not exhibit circular vibration. For V
atoms, their atomic vibrations are shown in Fig.4 (d). In this phonon mode,
except for V2 atoms, the other atoms remain relatively stationary. This
indicates that these V atoms predominantly contribute to the total energy of
harmonic vibration in this mode. Additionally, three adjacent V atoms vibrate
in opposite directions to the other three adjacent V atoms, resulting in a
total phonon polarization of zero. Comparing the vibrations of Sb atoms in
Fig.4 (a-c), V atoms exhibit a larger amplitude in one direction, although the
polarization trajectory forms an ellipse. Thus, as shown in Table.1, its
polarization value of V atom are smaller than those of Sb atoms in Fig.4 (c).
Moreover, we find that lager structural distortions induce greater
polarization of V2 atoms. For a detailed discussion, please refer to
Supplemental Materials Sec.S5.

_Local magnetic field.-_ Circular atomic vibration can result in phonon
magnetic moment, which, typically, shares a similar order of magnitude with
the nuclear magnetic moment [30]. Moreover, the $\mu$SR experiments have
detected local magnetic fields in AV${3}$Sb${5}$, with magnitudes also
comparable to the nuclear magnetic moment [13]. In this section, we explore
the potential attribution of the observed local magnetic field in experiments
to the phonon magnetic moment.

Figure 5:  The simulation of local magnetic field distribution by vibrations
of three neighbor V2 atoms at the Cs plane in Mode 77. (a) Vibrations of atoms
in V2 site in Mode 77. The point A is the projection of the center of the
triangle constructed by three V atoms along $z$-axis in Cs plane. (b),(c),(d)
The intensity distribution of the local magnetic field along the $z$, $x$ and
$y$ directions around A point, respectively.

The magnitude of the phonon moment correlate to the gyromagnetic ratio
$\gamma_{\alpha}=e\mathbf{Z}_{\alpha}/(2M_{\alpha})$. However, Born effective
charges $\mathbf{Z}_{\alpha}$ are ill-definition in metal, like CsV3Sb5 [33].
Nonetheless, Born effective charges reflect the amount of carried charges of
ions, and estimate them by calculating the charge density distribution. In the
Supplemental Materials Sec.S4, we provide support for the charge density
distribution of CsV3Sb5. Compared to Sb atoms, V atoms have a smaller
effective mass $M_{\alpha}$. Moreover, based on the charge density
distribution, V atoms harbor numerous localized electrons, resulting in a
larger magnitude of Born effective charges $\mathbf{Z}_{\alpha}$. Thus, V atom
has a greater gyromagnetic ratio $\gamma_{\alpha}$, leading a larger local
phonon magnetic moment in circular vibration.

Yu et al. [14] reported the detection of local magnetic fields along different
directions in the Cs plane following the CDW phase transition, as observed
through the $\mu$SR experiments. Here, we take Mode 77 as an example and
analyze the local magnetic fields generated by the polarizations of V atoms in
the Cs plane. Figure.5 (a) shows the vibrations of V atoms in this mode. As
mentioned in the previous section, the vibration trajectory of V atoms forms
an ellipse located in the $ab$-plane. To investigate the local field
distribution generated by the three neighboring V atoms in a triangular
configuration, we select a point A that is vertically aligned with the center
of the three V atoms in the Cs plane and simulate the nearby local field
distribution. We designate point A as the center of the $xy$-plane and
simulate the local magnetic field distribution in three different directions:
$B_{z},B_{x},B_{y}$, as shown in Figure.5 (b-d), respectively. The local
magnetic field along $z$ direction ($B_{z}$) exhibits the largest magnitude.
Furthermore, compared to the local magnetic field along $x$ direction
($B_{x}$), the local magnetic field in others directions are more localized
around A point.

According to the analysis above, when these anomalous phonons are excited, the
phonon moment of V atoms will induce a local magnetic field in Cs plane. In
the $\mu$SR experiment, muons decayed into positrons. These positrons carry
energy and specific momentum will scatter with phonons, potentially exciting
anomalous phonons. This phenomenon is akin to the excitation of chiral phonon
by the photons [34, 35] and also occurs in the AHE. Furthermore, we analyze
the magnitude of magnetic field intensity and magnetic moment induced by this
mode. We estimate that the magnetic field intensity approximately 10-5 T at
point A, which falls in the range of magnitude in order detected by the
$\mu$SR experiment. Meanwhile, the magnetic moment approaches
10${}^{-4}\mu_{B}$ per V atom. This value is comparable with the magnitude of
the magnetic moment reported in antiferromagnetic materials with considerable
AHE [36]. For a detailed discussion, please refer to the Supplemental
Materials Sec.S6. Thus, the phonon magnetic moment suggests an explanation on
the origin of local magnetic field observed in the $\mu$SR and AHE
experiments, and it is worth verifying experimentally.

_Discussions.-_ We have extended the theory of chiral phonon to the anomalous
phonon in kagome metals, the global chiral phonon in noncentrosymmetric
crystal structures is replaced by the local anomalous phonon in
centrosymmetric crystal structures, just as the global magnetization in
ferromagnet is replaced by sublattice magnetization in antiferromagnets. In
the CDW phase, we do not observe anomalous acoustic phonons for V atoms;
instead, we predict that the significant experimental characteristic will be
anomalous optical phonons. The possible experimental characteristics should be
either the Raman spectroscopy or the infrared circular dichroism spectra. As
we have shown in the proceeding section, for the former, the energy range of
anomalous phonons of V atom fall in 5 THz to 9 THz, and for the latter, the
presence of the local phonon magnetic moment will demonstrate in infrared
circular dichroism. We anticipate further experiments could verify the
prediction of the present theory.

There are two additional points that should be mentioned. (1) The vibrations
of V atoms in the neighboring kagome plane contribute to an opposite local
field in the Cs plane, resulting in the overall local field being zero along
the $x$ and $y$ directions. However, if the CDW phase adopts a $2\times
2\times 2$ configuration with $\pi$-phase shift rather than the $2\times
2\times 1$ configuration mentioned above, the entire local field in the Cs
plane along the $x$ and $y$ directions would be nonzero. For further details,
please refer to the Supplemental Materials Sec.S7. (2) Because of the time-
reversal symmetry, at +q and -q points, the vibration directions are opposite
for the same atom, resulting in a zero the macroscopical phonon moment.
Nevertheless, in the $\mu$SR experiment, muons decay into positrons, creating
an effective scattering cross section of positrons and leading to a
nonequilibrium state of the system. Consequently, it results in a different
number of phonons at $\pm$q points and a nonzero macroscopical phonon magnetic
moment in real space, which further causes the TRSB. Please see the
Supplemental Materials Sec.S8 for detailed discussion.

Yu et al. [16] and other authors [37, 38] have reported giant AHE in KV3Sb5
and RbV3Sb5, in their CDW phase, which has sparked significant interest. The
absence of local magnetic moment [16, 13] in these kagome metals ruled out the
skew scattering mechanism, while the intrinsic Berry phase alone is
insufficient to account for the observed anomalous conductivity in
experiments. This the necessity of a new mechanism. The presence of anomalous
phonons and their associated phonon magnetic moment could provide an
explanation for the microscopic origin of the AHE and the TRSB observed in
CsV3Sb5 and KV3Sb5. As demonstrated above, the polarized motions of V ions in
the CDW phase may contribute orbital magnetic moment of phonons, resulting in
extraordinary scattering to the Hall transport.

Moreover, as we address in the beginning, the debate surrounding the
microscopic origin of the TRSB could be spontaneously resolved within the
framework of present anomalous phonon and phonon magnetic moment. The orbital
magnetic moment arising from anomalous phonons contributes to the local
magnetic field, thereby breaking time-reversal symmetry. Hence, the concept of
anomalous phonons and phonon magnetic moments could naturally explain
microscopic origin of the puzzled AHE and the TRSB effects in the CDW phase of
the kagome metals AV3Sb5, elucidate other hidden magnetic properties in
distorted materials.

_Summary.-_ We perform the first-principles calculations to investigate the
atomic vibrations of the pristine and the CDW phases of CsV3Sb5. Both phases
maintain spatial inversion symmetry, thus chiral phonon do not emerge, but
anomalous phonons are present in certain high-symmetry paths. In the CDW
phase, the lattice distortion of kagome plane induces polarizations of V
atoms, with the magnitude of polarization dependent on the degree of
distortion. These polarizations manifest as elliptical vibrations of V atoms
within kagome plane, leading to the emergence of orbital magnetic moments and
contributing to the local magnetic field. The presence of orbital magnetic
moments in phonons offers a novel explanation for observations in the $\mu$SR
experiment, even the AHE in CsV3Sb5, and deserved to further investigate in
experimental study.

The authors thank the supports from the NSFC of China under Grant Nos.
11974354, Program of Chinese Academy of Sciences, Science Challenge Project
No. TZ2016001. Numerical calculations were partly performed at the Center for
Computational Science of CASHIPS, the ScGrid of the Supercomputing Center, the
Computer Network Information Center of CAS, the CSRC computing facility, and
Hefei Advanced Computing Center.

## References

  * Syozi [1951] I. Syozi, Prog. Theor. Phys. 6 (1951).
  * Ramirez [1994] A. P. Ramirez, Annu. Rev. Mater. Sci. 24, 453 (1994).
  * Yan  _et al._ [2011] S. M. Yan, A. H. David, and R. W. Steven, Science 332, 1173 (2011).
  * van Heumen [2021] E. van Heumen, Nat. Mater. 20, 1308 (2021).
  * Guo and Franz [2009] H. M. Guo and M. Franz, Phys. Rev. B 80, 113102 (2009).
  * Ko  _et al._ [2009] W. H. Ko, P. A. Lee, and X. G. Wen, Phys. Rev. B 79, 214502 (2009).
  * Ortiz  _et al._ [2019] B. R. Ortiz, L. C. Gomes, J. R. Morey, M. Winiarski, M. Bordelon, J. S. Mangum, I. W. H. Oswald, J. A. Rodriguez-Rivera, J. R. Neilson, S. D. Wilson, E. Ertekin, T. M. McQueen, and E. S. Toberer, Phys. Rev. Mater. 3, 094407 (2019).
  * Ortiz  _et al._ [2020] B. R. Ortiz, S. M. L. Teicher, Y. Hu, J. L. Zuo, P. M. Sarte, E. C. Schueller, A. M. M. Abeykoon, M. J. Krogstad, S. Rosenkranz, R. Osborn, R. Seshadri, L. Balents, J. He, and S. D. Wilson, Phys. Rev. Lett. 125, 247002 (2020).
  * Chen  _et al._ [2021a] H. Chen, H. Yang, B. Hu, Z. Zhao, J. Yuan, Y. Xing, G. Qian, Z. Huang, G. Li, Y. Ye, S. Ma, S. Ni, H. Zhang, Q. Yin, C. Gong, Z. Tu, H. Lei, H. Tan, S. Zhou, C. Shen, X. Dong, B. Yan, Z. Wang, and H. J. Gao, Nature 599, 222 (2021a).
  * Chen  _et al._ [2021b] K. Y. Chen, N. N. Wang, Q. W. Yin, Y. H. Gu, K. Jiang, Z. J. Tu, C. S. Gong, Y. Uwatoko, J. P. Sun, H. C. Lei, J. P. Hu, and J. G. Cheng, Phys. Rev. Lett. 126, 247001 (2021b).
  * Hu  _et al._ [2022] Y. Hu, S. M. Teicher, B. R. Ortiz, Y. Luo, S. Peng, L. Huai, J. Ma, N. C. Plumb, S. D. Wilson, J. He, and M. Shi, Sci. Bull. 67, 495 (2022).
  * Xie  _et al._ [2022] Y. F. Xie, Y. K. Li, P. Bourges, A. Ivanov, Z. J. Ye, J. X. Yin, M. Z. Hasan, A. Luo, Y. G. Yao, Z. W. Wang, G. Xu, and P. C. Dai, Phys. Rev. B 105, L140501 (2022).
  * Kenney  _et al._ [2021] E. M. Kenney, B. R. Ortiz, C. Wang, S. D. Wilson, and M. J. Graf, J. Phys.: Condens. Matter 33, 235801 (2021).
  * Yu  _et al._ [2021a] L. Yu, C. Wang, Y. Zhang, M. Sander, S. Ni, Z. Lu, S. Ma, Z. Wang, Z. Zhao, H. Chen, K. Jiang, Y. Zhang, H. Yang, F. Zhou, X. Dong, S. L. Johnson, M. J. Graf, J. Hu, H.-J. Gao, and Z. Zhao, (2021a), arXiv:2107.10714 .
  * Mielke Iii  _et al._ [2022] C. Mielke Iii, D. Das, J. X. Yin, H. Liu, R. Gupta, Y. X. Jiang, M. Medarde, X. Wu, H. C. Lei, J. Chang, P. Dai, Q. Si, H. Miao, R. Thomale, T. Neupert, Y. Shi, R. Khasanov, M. Z. Hasan, H. Luetkens, and Z. Guguchia, Nature 602, 245 (2022).
  * Yu  _et al._ [2021b] F. H. Yu, T. Wu, Z. Y. Wang, B. Lei, W. Z. Zhuo, J. J. Ying, and X. H. Chen, Phys. Rev. B 104, L041103 (2021b).
  * Shumiya  _et al._ [2021] N. Shumiya, M. S. Hossain, J.-X. Yin, Y.-X. Jiang, B. R. Ortiz, H. Liu, Y. Shi, Q. Yin, H. Lei, S. S. Zhang, G. Chang, Q. Zhang, T. A. Cochran, D. Multer, M. Litskevich, Z.-J. Cheng, X. P. Yang, Z. Guguchia, S. D. Wilson, and M. Z. Hasan, Phys. Rev. B 104, 035131 (2021).
  * Jiang  _et al._ [2021] Y.-X. Jiang, J.-X. Yin, M. M. Denner, N. Shumiya, B. R. Ortiz, G. Xu, Z. Guguchia, J. He, M. S. Hossain, X. Liu, J. Ruff, L. Kautzsch, S. S. Zhang, G. Chang, I. Belopolski, Q. Zhang, T. A. Cochran, D. Multer, M. Litskevich, Z.-J. Cheng, X. P. Yang, Z. Wang, R. Thomale, T. Neupert, S. D. Wilson, and M. Z. Hasan, Nat. Mater. 20, 1353 (2021).
  * Zhao  _et al._ [2021] H. Zhao, H. Li, B. R. Ortiz, S. M. L. Teicher, T. Park, M. Ye, Z. Wang, L. Balents, S. D. Wilson, and I. Zeljkovic, Nature 599, 216 (2021).
  * Feng  _et al._ [2021] X. Feng, K. Jiang, Z. Wang, and J. Hu, Sci. Bull. 66, 1384 (2021).
  * Denner  _et al._ [2021] M. M. Denner, R. Thomale, and T. Neupert, Phys. Rev. Lett. 127, 217601 (2021).
  * Li  _et al._ [2022a] H. Li, H. Zhao, B. R. Ortiz, T. Park, M. Ye, L. Balents, Z. Wang, S. D. Wilson, and I. Zeljkovic, Nat. Phys. 18, 265 (2022a).
  * Li  _et al._ [2022b] H. Li, S. Wan, H. Li, Q. Li, Q. Gu, H. Yang, Y. Li, Z. Wang, Y. Yao, and H.-H. Wen, Phys. Rev. B 105, 045102 (2022b).
  * Zhang and Niu [2015] L. Zhang and Q. Niu, Phys. Rev. Lett. 115, 115502 (2015).
  * Zhu  _et al._ [2018] H. Zhu, J. Yi, M. Li, J. Xiao, L. Zhang, C. Yang, R. A. Kaindl, L. J. Li, Y. Wang, and X. Zhang, Science 359, 579 (2018).
  * Ptok  _et al._ [2021] A. Ptok, A. Kobialka, M. Sternik, J. Lazewski, P. T. Jochym, A. M. Oles, S. Stankov, and P. Piekarz, Phys. Rev. B 104, 054305 (2021).
  * Wang  _et al._ [2021] H. Wang, Z. C. Zhou, H. Chen, C. Q. Xia, L. F. Zhang, and X. Li, J. Phys.: Condens. Matter 33, 285704 (2021).
  * Chen  _et al._ [2019] H. Chen, W. Wu, S. Yang, X. Li, and L. Zhang, Phys. Rev. B 100, 094303 (2019).
  * Luo  _et al._ [2023] J. Luo, T. Lin, J. Zhang, X. Chen, E. R. Blackert, R. Xu, B. I. Yakobson, and H. Zhu, Science 382, 698 (2023).
  * Juraschek and Spaldin [2019] D. M. Juraschek and N. A. Spaldin, Phys. Rev. Mater. 3, 064405 (2019).
  * Xiong  _et al._ [2022] G. Xiong, H. Chen, D. Ma, and L. Zhang, Phys. Rev. B 106, 144302 (2022).
  * Tan  _et al._ [2021] H. Tan, Y. Liu, Z. Wang, and B. Yan, Phys. Rev. Lett. 127, 046401 (2021).
  * Gonze and Lee [1997] X. Gonze and C. Lee, Phys. Rev. B 55, 10355 (1997).
  * Cheng  _et al._ [2020] B. Cheng, T. Schumann, Y. Wang, X. Zhang, D. Barbalas, S. Stemmer, and N. P. Armitage, Nano Letters 20, 5991 (2020), pMID: 32633978\.
  * Nova  _et al._ [2017] T. F. Nova, A. Cartella, A. Cantaluppi, M. Först, D. Bossini, R. V. Mikhaylovskiy, A. V. Kimel, R. Merlin, and A. Cavalleri, Nature Physics 13, 132 (2017).
  * Nakatsuji  _et al._ [2015] S. Nakatsuji, N. Kiyohara, and T. Higo, Nature 527, 212+ (2015).
  * Yang  _et al._ [2020] S.-Y. Yang, Y. Wang, B. R. Ortiz, D. Liu, J. Gayles, E. Derunova, R. Gonzalez-Hernandez, L. Šmejkal, Y. Chen, S. S. P. Parkin, S. D. Wilson, E. S. Toberer, T. McQueen, and M. N. Ali, Sci. Adv. 6, eabb6003 (2020).
  * Yin  _et al._ [2021] Q. Yin, Z. Tu, C. Gong, Y. Fu, S. Yan, and H. Lei, Chin. Phys. Lett. 38, 037403 (2021).


# Singlet exciton dynamics of perylene diimide and tetracene based
hetero/homogeneous substrates via an ab initio kinetic Monte Carlo model

A. Manian<EMAIL_ADDRESS>These authors contributed equally.  ARC
Centre of Excellence in Exciton Science, School of Science, RMIT Univeristy,
Melbourne, 3000, Australia.  F. Campaioli<EMAIL_ADDRESS>These authors contributed equally.  ARC Centre of Excellence in Exciton
Science, School of Science, RMIT Univeristy, Melbourne, 3000, Australia.  I.
Lyskov  ARC Centre of Excellence in Exciton Science, School of Science, RMIT
Univeristy, Melbourne, 3000, Australia.  J. H. Cole  ARC Centre of Excellence
in Exciton Science, School of Science, RMIT Univeristy, Melbourne, 3000,
Australia.  S. P. Russo<EMAIL_ADDRESS>ARC Centre of Excellence in
Exciton Science, School of Science, RMIT Univeristy, Melbourne, 3000,
Australia.

###### Abstract

Luminescent solar concentrators (LSCs) are devices that trap a portion of the
solar spectrum and funnel it towards photon harvesting devices. The modelling
of LSCs at a quantum chemical level however, remains a challenge due to the
complexity of exciton and photon dynamic modelling. This study examines
singlet exciton dynamics occurring within a typical LSC device. To do this, we
use a rejection-free kinetic Monte Carlo method to predict diffusion lengths,
diffusion coefficients, substrate anisotropy, and average exciton lifetimes of
perylene diimide (PDI) and tetracene based substrates in the low concentration
scheme. Ab initio rate constants are computed using time-dependant density
functional theory based methods. PDI type substrates are observed to display
enhanced singlet exciton transport properties when compared to tetracene.
Simulations show that substrates with dipole-aligned chromophores are
characterised by anisotropic exciton diffusion, with slightly improved
transport properties. Finally, a PDI-tetracene substrate is simulated for both
disordered and dipole-aligned chromophore configurations. In this multi-dopant
substrate transport is predominantly mediated by PDI due to the asymmetry in
the transport rates between the two dyes considered. We conclude discussing
the properties of multi-dopant substrates and how they can impact the design
of next generation LSCs.

††preprint: AIP/123-QED

##  I Introduction

Luminescent Solar Concentrators (LSCs) are devices which use the concept of
exciton hopping and total internal reflection to trap light from a large
surface area, and guide photons to photovoltaic (PV) cells of a smaller
surface area[1, 2, 3, 4]. Here, total internal reflection allows for transport
of energy across macroscopic scales, while exciton dynamics mediate the
efficient transfer of energy between molecular bodies, in addition to the
recycling of otherwise unusable excitions via multi-exciton interactions[5, 6,
7, 8, 9]. LSCs are obtained by embedding chromophores in a transparent matrix
substrate with a large refractive index. When a photon enters the LSC, it is
absorbed by a chromophore. Then, upon photon emission, such as fluorescence,
the photon is guided within the matrix due to a condition of total internal
reflection[10, 11]. Light is thus concentrated within the matrix, and is
funnelled towards the edges, where PV cells are fixed to absorb each
photon[1].

While the idea at the core of this strategy is simple, a few hurdles
complicate the design of effective LSCs. On one hand, a large concentration of
chromophores within the substrate is required such that photons will always
interact with the device. However within such an amorphous matrix, the
chromophores tend to aggregate and clump together, creating undesired sources
of luminescent quenching[12, 13], where quenching refers to processes
competing with radiative pathways. On the other hand, if the concentration of
chromophores is too low, then the majority of photons are lost as they pass
through the LSC[10]. Other sources of loss include photons escaping the
device[11] upon fluorescence, and the intrinsic photoluminescence quantum
yield (PLQY) loss, i.e., the probability of fluorescence of the host
chromophore. Previously, we have investigated the parameters behind excitonic
quenching[12, 13, 13, 14], providing us a strong foundation with which to
consider the processes within a given device.

Modelling of an LSC device at a quantum level can be split into two
categories. The first concerns modelling of exciton dynamics, i.e., how long
the exciton takes to dissociate and return to a photon. The second concerns
modelling of the photon dynamics, i.e., when an exciton has dissociated into a
photon, how that energy packet propagates through the substrate, whether it is
reabsorbed by another chromophore, by the PV, or escapes. The former requires
an understanding of the quantum chemistry and corresponding rate constants of
the numerous de-excitation pathways, while the latter requires tracking of
photons through the substrate, or ray tracing[15, 16].

The dynamics of both interacting and non-interacting excitons have been
simulated using kinetic Monte Carlo (KMC) methods in a variety of studies,
providing insights on the fundamental features of exciton transport as well as
practical guidelines for the design of optoelectronics and PV applications
[17, 18, 19, 20, 5, 21, 22, 23, 24]. Engels in particular has studied exciton
diffusion and charge-carrier mobilities using Marcus theory[25, 26, 27].
However, it was noted that diffusion properties calculated using Marcus theory
were drastically underestimated with respect to experiment and other
methods[28], despite being qualitatively correct with respect to trends.
Further, the Monte Carlo approach allows simulation of transport properties
within a time dependent framework. As such, a method involving computation of
the full spectral density within a time dependent framework may allow for
accurate results comparable with experiment.

While modelling of the photons’ dynamics in the substrate has been studied
rigorously[4, 29, 30, 31], a treatment at the atomic scale is very complex; we
will therefore address this problem in a future study. In this paper we focus
on exciton migration through a medium such as an LSC. To do so, we use a KMC
method based on transition rates that we obtain from quantum chemical
calculations [14]. KMC methods have been used in surface based
simulations[32], and exciton diffusion/hopping studies[33], making it an ideal
tool for our purpose. Herein, we simulate exciton diffusion and dissociation
in homogeneous and heterogeneous Perylene Diimide (PDI) and Tetracene
chromophores in a Toluene medium. PDIs[13, 10] are important in the field of
organic light emitting diodes (OLEDs) and LSCs, typically displaying large
PLQYs, while Tetracene[34, 35, 36] has been well studied as a host for singlet
fission type applications. We investigate transport properties such as the
average diffusion lengths, diffusion coefficient, and the anisotropy profile
of the exciton migration across disordered chromophore configurations with
both random and forcibly aligned transition dipoles. For convenience, we refer
to these two types of systems as disordered configurations and dipole-aligned
configurations respectively. We also study the frequency of different
dissociation pathways within favoured dissociation sites, and the average
exciton lifetimes, which can provide valuable guidelines for LSCs design.

##  II Theory & Implementation

###  II.1 Transport Model

Exciton transport, often referred to as exciton diffusion, is one of the key
processes behind the operation of organic optoelectronic devices[37, 38, 33].
In addition to the exciton transport, we also need to model the processes that
result in loss of the electron-hole quasiparticle, such as charge
recombination, internal conversion (IC) and inter-system crossing (ISC). These
latter two are relaxation processes mediated by heat, in which heat is
exchanged with the environment when spin is and is not conserved [14, 13]
respectively.

In essence, exciton transport corresponds to the the migration of an energy
packet from some molecule A to another molecule B, simply shown as:

${}^{1}A^{\star}+{}^{1}B\rightarrow{}^{1}A+{}^{1}B^{\star}.$ (1)

Here, an exciton residing on chromophore A is transferred, or hops, to
chromophore B. Exciton hopping can be modelled within both Coulombic and
Exchange picture, using either Förster theory[39] or Dexter theory[40]
respectively. In this work, we only consider Förster type hopping, which is
modelled using rigorous quantum chemical calculations, as discussed in the
next section.

Figure 1: Energy diagram illustrating how the system excites and relaxes.
Here, an incident photon excited the ground state chromophore to an excited
state, be it the first or second singlet excited state. Following relaxation
to the ground vibrational state of first singlet electronic excited state, the
energy can transition between nearby electronic excited states.

Two chromophores that undergo energy transfer are assumed to be in resonance,
where molecule A donates an exciton to an accepting molecule B, which is known
as Förster Resonance Energy Transfer (FRET). The transition rate of a FRET
mechanism can be modelled using Fermi’s Golden Rule, where the matrix element
can be shown using the simple dipole approximation[41], given as:

$\left|\left<f\left|H^{\left(1\right)}\right|i\right>\right|=\frac{1}{4\pi\epsilon_{0}}\frac{\kappa\left|\mu_{D}\right|\left|\mu_{A}\right|}{n^{2}r^{3}},$
(2)

where $\mu_{D}$ ($\mu_{A}$) is the transition dipole moment of a given donor
(acceptor) chromophore, $\epsilon_{0}$ is the permittivity of free space,
$\kappa$ is the dipole angular orientation factor, given in Eq. (4), $n$ is
the refractive index of the medium, and $r$ is the interchromophore
separation. The density of states across the energy continuum $\rho$ can be
expressed as the integral of the normalised absorption band of the acceptor
molecule and emission band of the donor, $\textrm{Abs}\left(\omega\right)$ and
$\textrm{Emi}\left(\omega\right)$ respectively, weighted by the PLQY of the
donor in the absence of the acceptor, $Q_{D}$:

$\rho=Q_{D}\int\textrm{Abs}\left(\omega\right)\textrm{Emi}\left(\omega\right)d\omega.$
(3)

The angular orientation factor $\kappa$ can be expressed as:

$\kappa=\hat{\mu}_{A}\cdot\hat{\mu}_{D}-3\left(\hat{\mu}_{D}\cdot\hat{R}\right)\left(\hat{\mu}_{A}\cdot\hat{R}\right),$
(4)

where $\hat{R}$ is the normalized inter-chromophore separation, and where all
three vectors are unit vectors. We then substitute both Equation 2 and
Equation 3 into Fermi’s Golden Rule to yield the following:

$\begin{split}k_{FRET}=&\frac{2\pi}{\hbar}\frac{1}{\left(4\pi\epsilon_{0}\right)^{2}}\frac{\kappa^{2}\left|\mu_{D}\right|^{2}\left|\mu_{A}\right|^{2}}{n^{4}r^{6}}\times\\\
&\times
Q_{D}\int\textrm{Abs}\left(\omega\right)\textrm{Emi}\left(\omega\right)d\omega.\end{split}$
(5)

All terms within the integrand are often collected together and termed the
spectral overlap.

When we consider excitonic quenching, we examine the probability of radiative
decay upon charge recombination. This comes back to the quantum yield $Q_{D}$
of the chromophore. As per Kasha’s Rule[42], fluorescence will occur from the
first singlet excited state. Therefore, the quantum yield of fluorescence is a
ratio of radiative decay rates against the combined rate of all possible
mechanisms, as shown in Figure 1. In the case of an isolated monomer, the
likely competing mechanisms are IC and ISC. Each process is dependant on
different photophysical coupling terms: fluorescence on the transition dipole
moment, IC on vibronic coupling between singlet states, and ISC the spin-orbit
coupling between singlet and triplet states. These rate constants can be very
difficult to compute, however we have recently developed a methodology to do
so at the Time-Dependant Density Functional Theory (TDDFT) level of
theory[14].

###  II.2 The KMC method

KMC is a numerical method to simulate the dynamics of a process, based on
known transition rates between the states of a considered system. KMC methods
are widely adopted in many areas of Physics and Chemistry due to their
computational efficiency, and can be used to simulated both equilibrium and
non-equilibrium processes [43, 44, 45, 46, 47, 48, 49, 50, 51, 52]. They
provide an ideal approach to simulate FRET and Dexter exciton transport
processes in disordered systems, which are known to correspond to a classical
random walk over an ensemble of sites [53, 38]. Moreover, when the occurrence
of multi-exciton interactions processes, such as triplet fusion and singlet
fission, can be neglected, each exciton can be treated independently, turning
the KMC simulation into an embarrassingly parallel computational problem.

In this work we adopt a rejection-free KMC method to simulate the dynamics of
an ensemble of non-interacting excitons that undergo FRET-mediated transport,
fluorescence, IC or ISC. The molecular geometries and rates associated with
these processes were computed using a combination of DFT and DFT based
multireference configuration interaction DFT/MRCI approaches[54]. The
considered systems are static arrangements of PDI and Tetracene chromophores
in a substrate. This is done to simulate the transport of excitons in a solid,
where chromophores’ position and dipole orientation do not vary during the
transport process. The photo-physical properties of the substrate are based on
Toluene — despite it being a liquid solvent —, due to its high refractive
index, ideal for the requirements of total internal reflection necessary for
the LSCs working mechanism. For each considered system, we simulate transport
within a cubic portion of the LSCs, whose size is such that the probability of
transport outside the simulation box is negligible. In all cases we consider
concentrations corresponding to a low concentration solution, such that
chromophores are close enough to maximise exciton transfer, but not so close
such that exciton quenching pathways facilitated by molecular aggregation can
be considered negligible.

It should be noted that since we do not consider photon dynamics, we therefore
do not consider photon reabsorption, a process whereby exciton dissociation at
site $\vec{r}_{i}$ results in a photon propagating through the substrate, and
is reabsorbed by another site $\vec{r}_{j}$. Rather, simulations begin and
terminate at the birth and death of a given exciton, due to the foreseen
complexity of ray tracing at the quantum chemical level.

First, we generate the geometry of the considered systems. A number of sites,
corresponding to the chromophores locations, is randomly generated to match
the dyes concentration within the simulation box. The sites are uniformly
distributed in the considered volume. Then, the dipole orientation of each
chromophore is randomly generated, with different distributions being used for
disordered and dipole-aligned materials. Once the geometry of the systems is
fixed, exciton dynamics is simulated using the KMC method. A state
$\mathcal{S}_{t}=\\{t,\vec{r}_{i}\\}$ of the process at time $t$ is given by
the location (site) $\vec{r}_{i}$ of a singlet exciton. Exciton transport is
simulated propagating the state in time
$\mathcal{S}_{t}\to\mathcal{S^{\prime}}_{t+\Delta t}$, using the routine
outlined in Alg. 1, until a dissociation event occurs. Transport occurs via
FRET from the current site $\vec{r}_{i}$ to a site $\vec{r}_{j}$ within the
interaction range $r_{\text{max}}$, which depends on the dye and is determined
as specified in Sec. IV.1.

1

input : The current state $\mathcal{S}_{t}$.

output : The propagated state $\mathcal{S}^{\prime}_{t+\Delta t}$ if FRET
occurs, the dissociation pathway otherwise.

2

3Calculate the FRET transition rates $k_{FRET}(i\to j)$ from the current site
$\vec{r}_{i}$ to any site $\vec{r}_{j}$ within interaction range
$r_{\text{max}}$;

4Determine the dissociation rates $k_{r}$, $k_{IC}$, and $k_{ISC}$, which
depend on the chromophore species of the current site $\vec{r}_{i}$;

5Calculate the cumulative function $R_{m}=\sum_{m=1}^{M}k_{m}$ of all the $M$
rates involved, and the cumulant $Q=R_{M}$;

6Sample a uniform random number $u\in(0,1]$;

7Determine which event occurs, by choosing the rate $k_{m}$ such that
$R_{m-1}\leq uQ\leq R_{m}$;

8if  _$k_{m}=k_{FRET}(i\to j)$_ then

9 FRET occurs, the new site is $\vec{r}_{j}$;

10 Sample another uniform random number $u^{\prime}\in(0,1]$;

11 Determine the time-interval $\Delta t=Q^{-1}\ln(1/u^{\prime})$;

12 Return the new state $\mathcal{S}^{\prime}_{t+\Delta t}$.

13else

14 Dissociation occurs via $k_{r}$, $k_{IC}$ or $k_{ISC}$;

15 Return the dissociation pathway, end the routine.

16

Algorithm 1 Rejection-free KMC routine

##  III Computational Details

###  III.1 Systems generation and KMC implementation

Using the KMC method described in the previous section we simulated exciton
transport for single and multi-dopant substrates for two classes of 3D
systems: disordered chromophore configurations with both random and aligned
transition dipoles, labelled as disordered configurations and dipole-aligned
configurations respectively. The former case of disordered geometries were
chosen as they most closely represented exciton diffusion between nanofibres
in bulk[33]. The latter case corresponds to the scenario for which most
dipoles are aligned and parallel to a given axis. Despite of the fact that
dipole alignment methods can vary greatly[55], the resulting geometries allow
for possible circumvention of high dye concentration loss mechanisms[29, 29,
56, 9].

Table 1:  Calculated rates and transfer properties for Perylene Diimide and Tetracene monomers simulated in Toluene.  Species |  $\mu_{D}$ $(au)$ |  $\mu_{A}$ $(au)$ |  $J$ $(au^{-1})$ | $Q_{D}$ |  $k_{r}$ $(s^{-1})$ |  $k_{IC}$ $(s^{-1})$ |  $k_{ISC}$ $(s^{-1})$ |  $R_{C}$ $(nm)$ |  $R_{I}$ $(nm)$
---|---|---|---|---|---|---|---|---|---
PDI | 4.005 | 3.836 | 46.243 | 0.99 | $1.505\times 10^{8}$ | $1.316\times 10^{4}$ | $1.025\times 10^{6}$ | 1.6 | 15
Tetracene | 1.340 | 1.299 | 31.734 | 0.26 | $2.245\times 10^{7}$ | $2.707\times 10^{7}$ | $3.565\times 10^{7}$ | 1.4 | 8

To generate each arrangement, chromophore concentrations were set to
$6.02\times 10^{6}~{}\mu m^{-3}$ and $8.03\times 10^{6}~{}\mu m^{-3}$ for
homogeneous PDI and Tetracene systems, respectively (see Supplementary
Information for further details on heterogeneous chromophore concentrations).
Simulations were carried out in a cubic box with side length of $0.25~{}\mu
m$. Random sites were uniformly generated to match the desired concentration.
An iterative approach was used to ensure that inter-chromophore separations
were always larger than the culling sphere radius $r_{\text{min}}$, which
depends on the chromophore species, as discussed in Sec. IV. This was done by
repopulating the simulation box after removing the sites that were too close
too each other.

For disordered arrangements, dipole orientations
$\hat{\mu}=(\cos\theta\sin\phi,\sin\theta\sin\phi,\cos\phi)$ were randomly
generated by uniformly sampling angles $\theta\in(-\pi,\pi]$ and
$\phi\in(0,\pi]$, to obtain a uniform distribution on the sphere. For the case
of dipole-aligned arrangements, the two angles were sampled uniformly in
$\theta\in(-p\pi,p\pi]$ and $\phi\in((1-p)\pi/2,(1+p)\pi/2]$, with $p=0.01$,
thus returning dipoles mostly aligned with the $\hat{x}$ direction. Defects
were simulated by randomly changing $p\to 1$, with defect probability of
$1\%$.

For every KMC trajectory a singlet exciton was initialised at time $t_{0}=0$
ns, on a random site $\vec{r}_{i}$ within a sphere of radius $2r_{\text{max}}$
located in the centre of the simulation box. Initial states
$\mathcal{S}_{t_{0}}$ were then propagated with the KMC routine describe in
Alg. 1 until a dissociation event occurred. For each system, 5000 trajectories
were sampled.

Figure 2: (_Top_) Normalised absorption and emission spectra for homogeneous
A$\left)\right.$ Perylene Diimide and B$\left)\right.$ Tetracene systems
solvated in Toluene, in atomic units. (_Bottom_) Spectral overlap schemes for
heterogeneous 2 dye substrates when the exciton transport path is from
C$\left)\right.$ Perylene Diimide to Tetracene, and D$\left)\right.$ Tetracene
to Perylene Diimide. Paired overlap functions are computed upon combination of
absorption and emission spectra, as per Equation 2.  Figure 3: Excitonic
transfer rate computed using Förster theory as a function of interchromophore
separation for homogeneous PDI and Tetracene (TET) substrates. The horizontal
lines represent the combined dissociation rates
$k_{\text{diss}}=k_{r}+k_{IC}+k_{ISC}$ for the two dyes. Interaction ranges
for PDI and Tetracene are represented by the vertical dashed lines, while the
culling radii, equal to 16 Å and 14 Å for PDI and Tetracene, respectively, are
approximately represented by the solid black vertical line.

###  III.2 Quantum Chemistry

Quantum chemically optimised geometries of the relevant excited states, rate
constants, and photophysical properties of PDI and tetracene were computed as
per our previous study[14]. Here, electronic ground and excited state
geometries were computed at the DFT level of theory using the Becke
3-parameter Lee–Yang–Parr exchange-correlation hybrid functional B3LYP[57]
using the valence triple-zeta polarization basis set TZVP[58] as implemented
in the Gaussian16 software package[59]. A solvent of Toluene was simulated
using a polarizable continuum model (PCM). The same parameters were used to
compute the electronic Hessian. Single-point calculations were performed using
DFT/MRCI, using the def2-TZVP basis set. The one-particle basis was computed
using the Becke half-and-half Yang-Lee-Parr BHLYP exchange-correlation
functional[60] as implemented in the Turbomole software package[61]. Here, the
PCM was employed using the Cosmo module[62]. Spin-orbit matrix elements were
computed using the Spock.ci module[63, 64, 65] of the DFT/MRCI platform.
System temperature was set to 300 K, corresponding to standard room
temperature conditions.

The culling radius was derived from the Lennard-Jones potential. Here, each
atom is modelled as a sphere, with radius equal to it’s collision radius
dependant on it’s isotropic static polarizability. Following, a box is built
around these spheres, such that the volume of this structure in minimised,
hugging the spheres[66, 67]. The culling radius is taken as the largest
dimension of this box, and is computed using the Chemcraft visualisation
package[68].

##  IV Results & Discussion

###  IV.1 Quantum Chemistry

Expansion of the DFT/MRCI wavefunction for PDI shows a strong
$S_{0}\rightarrow S_{1}$ transition from the highest unoccupied molecular
orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO), agreeing
well with Yang & Jang[69]. Here, absorption and emission energies were
computed to be 2.36 eV and 2.09 eV respectively, with a corresponding
adiabatic energy of 2.23 eV. Absorption and emission spectra, shown in Figure
2A, compares well with experiment[10]. While typically at the Franck-Condon
point we expect a degree of symmetry with respect to the zero-phonon line of
absorption to emission spectra, here we see a very slight shift in the density
of states in the absorption spectra, with a slightly smaller secondary peak
with respect to the emission bandstructure. We attribute this to the switching
of intensities between two vibrational normal modes at the Franck-Condon
points of the two transitions, more commonly known as the Duschinsky
effect[70]. This change does not appear to affect the rate of either exciton
transfer or de-excitation pathways. Tetracene also yielded a strong
HOMO$\rightarrow$LUMO transition for the first singlet excited state upon
DFT/MRCI wavefunction expansion. Absorption and emission energies were found
to be 2.61 eV and 2.28 eV, with a resulting adiabatic energy of 2.44 eV.
Similar to PDI, the spectral bandshapes compare well to literature[36], while
also displaying some slight asymmetry about the 0-0 phonon line, as shown in
Figure 2B.

Rate constants and important quantum chemical qualities are reported in Table
1. For a homogeneous PDI exciton transfer pathway, the computed bandshapes
result is a large spectral overlap of 46 $au^{-1}$. Further, donor and
acceptor transition dipole moments of 4.00 au and 3.84 au were measured. For
the three monomer deexciatation pathways, rate constants of $1.505\times
10^{8}~{}s^{-1}$, $1.316\times 10^{4}~{}s^{-1}$, and $1.025\times
10^{6}~{}s^{-1}$ were used for fluorescence, IC and ISC respectively, as per
Ref. 14, resulting in a near-unity fluorescence quantum yield of 0.99.

Tetracene on the other hand displays significantly weaker fluorescent
properties. While displaying a sizeable homogeneous overlap of 32 $au^{-1}$,
the donor and acceptor transition dipole moments are almost four times smaller
than that of PDI, computed to be 1.34 $au$ and 1.30 $au$ respectively. The
computed rate constants of $2.245\times 10^{7}~{}s^{-1}$, $2.707\times
10^{7}~{}s^{-1}$, and $3.565\times 10^{7}~{}s^{-1}$ for fluorescence, IC, and
ISC respectively result in a much smaller quantum yield of 0.26. It is worth
noting here that while photophysical properties computed for Tetracene in
Toluene here are very similar to those computed in CycloHexane[14], smaller
derivative components for the second order corrections result in a
significantly smaller ISC rate constant, and therefore a larger quantum yield.

Upon comparison to experimental data reported by Burgdorff and coworkers[71],
we see that our predicted PLQY is slightly overestimated with respect to
experiment. The fast fluorescence lifetime noted as 4.8 ns infers the
stabilisation of higher lying excited states, increasing the effect of second
order terms. However, Burdett and coworkers[36] noted a range of lifetimes
spanning 20-100 ns. This suggests that aggregation in Tetracene has a negative
effect on the photophysics of the substrate, with strong aggregation slowing
down the fluorescence lifetime. As such, the rates computed here can be deemed
accurate with respect to a chromophore in isolation.

Table 2:  Computed spectral overlap for multi-dopant substrates. Donors are given in column, while acceptors are given in row. Overlaps in units of $au^{-1}$.  Transfer Path $\nearrow$ | PDI | Tetracene
---|---|---
PDI | 46.243 | 0.056
Tetracene | 59.875 | 31.734

If we examine how the exciton hopping rate of an example PDI dimer system
evolves as a function of the interchromophore separation, as shown in Figure
3, we can see that any separation beyond 50 nm yields rate constants smaller
than $10^{3}~{}s^{-1}$. Here, all terms used for computation of the rate
constant are fixed, as per Equation 5, assuming this test system is
anisotropic[72] and taking the standard value of $\kappa^{2}=\frac{2}{3}$ for
the orientation. With $r_{\text{min}}$ as a limiting factor, measured here as
16 Å (see Supplementary Information), we may then define a range of inter-
chromophore separations in which excitonic processes dominate. Here, we note
that the rate drops below $10^{5}~{}s^{-1}$ for separations larger than 15 nm.
In this regime, exciton dissociation is the dominating process. We will
therefore adopt this as the maximum to the interaction range $r_{\text{max}}$
for PDI chromophores. The same logic applies for a tetracene anisotropic
dimer, where the exciton hopping rate drops to less than $10^{5}~{}s^{-1}$
when dimer separation increases beyond $r_{\max}=8$ nm. Furthermore, the
culling sphere radius $r_{\text{min}}$ for Tetracene can be taken as 14 Å (see
Supplementary Information).

In this study we also consider a novel system whereby an even mixture of the
PDI and Tetracene dyes are in the same substrate. For this system, transition
probabilities will not be symmetric upon the exchange of donor and acceptor,
and must be calculated separately. Our results show that PDI-Tetracene
displays a drastically different overlap to Tetracene-PDI, as shown in Table 2
and Figure 2. More specifically, for the transfer of an exciton from PDI-
Tetracene, a very poor overlap of 0.056 $au^{-1}$ is observed, as shown in
Figure 2C. Here, the respective donor and acceptor spectra are not compatible,
suggesting that PDI-PDI exciton hopping is more favourable.

However, in the reverse case, where an exciton propagates via a Tetracene-PDI
pathway, we see an appreciable overlap of 60 $au^{-1}$, as shown in Figure 2D.
Here, the primary peaks of both respective donor and acceptor spectra overlap
considerably, resulting a less intense overlap across a larger portion of the
energetic continuum. This larger overlap suggests that an exciton pathway of
Tetracene-PDI is in fact more favourable than homogeneous hopping, with very
interesting implications on exciton transport properties for the case of mixed
dyes. In this particular case, Tetracene has significantly poorer
photophysical properties with respect to PDI. Therefore, in this 2-dye system,
an exciton is less likely to hop from a PDI onto a Tetracene, and if it does,
is likely to jump back. This asymmetry can be exploited, when trying to ensure
charge recombination via a particular pathway, in this case fluorescence.

###  IV.2 Exciton transport properties

Table 3: Exciton transport properties for disordered and dipole-aligned
configurations. Generalised diffusion coefficient $D$ and diffusion anomaly
$\alpha$ have been obtained fitting the data to the models of Eq. (7) and Eq.
(8), respectively. Mean exciton lifetime $\langle\tau\rangle$ and average
diffusion length $\langle r\rangle$ have been obtained averaging over the KMC
trajectories. Transport anisotropy is shown via the $D_{x}$, $D_{y}$, and
$D_{z}$ components of $D$. The ratio of each dissociation pathways is shown
for radiative (Rad), inter-conversion (IC), and inter-system crossing (ISC).
See sec. IV.2 for further details. Disordered Configurations:
---
| Species | $D\;(\text{nm}^{2}/\text{ns})$ | $\alpha$ |  $\langle\tau\rangle$ $(\text{ns})$ |  $\langle r\rangle$ $(\text{nm})$ | Anisotropy $(D_{x},D_{y},D_{z})$ | Rad (%) | IC (%) | ISC (%)
---|---|---|---|---|---|---|---|---
PDI | 2.82 | 0.94 | 6.51 | 4.76 | 0.33, 0.32, 0.34 | 99.34 | 0.64 | 0.02
TET | 0.02 | 0.84 | 11.89 | 0.46 | 0.33, 0.32, 0.34 | 26.12 | 42.82 | 31.06
PDI + TET | 0.89 | 0.94 | 7.65 | 2.78 | 0.33, 0.33, 0.33 | 83.14 | 9.80 | 7.06
Dipole-aligned Configurations:
| Species | $D\;(\text{nm}^{2}/\text{ns})$ | $\alpha$ |  $\langle\tau\rangle$ $(\text{ns})$ |  $\langle r\rangle$ $(\text{nm})$ | Anisotropy $(D_{x},D_{y},D_{z})$ | Rad (%) | IC (%) | ISC (%)
---|---|---|---|---|---|---|---|---
PDI | 3.42 | 0.95 | 6.66 | 5.23 | 0.51, 0.24, 0.25 | 99.34 | 0.64 | 0.02
TET | 0.03 | 0.81 | 11.58 | 0.5 | 0.50, 0.24, 0.26 | 26.04 | 41.58 | 32.38
PDI + TET | 1.10 | 0.97 | 7.56 | 3.15 | 0.49, 0.25, 0.26 | 85.76 | 8.49 | 5.74
Figure 4: Exciton transport in single-dopant substrates for disordered
(_left_) and dipole-aligned (_right_) configuration. Singlet transport is more
efficient in PDI (_top_) than in Tetracene (_bottom_), despite the former
being characterised by shorter exciton lifetimes. Dipole-aligned
configurations have enhanced transport properties and a slight anisotropy
towards the dipoles’ orientation, which coincides with the $x$-axis.
Generalised diffusion coefficient $D$ and diffusion anomaly $\alpha$ are
fitted to the data as discussed in Sec. IV.2. Figure 5: Exciton transport in
multi-dopant substrates for disordered (_left_) and dipole-aligned (_right_)
configuration. Singlet transport in PDI-Tetracene mixtures is vastly more
efficient than in Tetracene, as excitons spend the majority of their lifetime
on PDI chromophores. As for the case of single-dopant substrates, dipole-
aligned configurations have enhanced transport properties and display
anisotropy towards the dipoles’ orientation, given by the $x$-axis.
Generalised diffusion coefficient $D$ and diffusion anomaly $\alpha$ are
fitted to the data as discussed in Sec. IV.2.

The results of KMC simulation were analysed to obtain average transport
properties. The mean square displacement

$\langle
r^{2}(t)\rangle=\frac{1}{N}\sum_{i=1}^{N}\lVert\vec{r}_{i}(t)-\vec{r}_{i}(t_{0})\rVert^{2},$
(6)

was calculated for each ensemble of $N=5000$ trajectories, where
$\vec{r}_{i}(t)$ is the position of the exciton at time $t$ for the $i$-th
trajectory. The mean square displacement was then used to infer the diffusion
transport properties according to the following model:

$\langle
r^{2}(t)\rangle=2D\alpha^{2}k^{-\alpha}\gamma(\alpha,kt)\frac{\Gamma(\alpha)}{\Gamma(\alpha+1)},$
(7)

which represents the time-evolution of the mean square displacement for
anomalous diffusion with generalised diffusion coefficient $D$, diffusion
anomaly $\alpha$ (regular diffusion for $\alpha=1$, sub-diffusive transport
for $\alpha<1$), and dissociation rate $k$, while the special functions
$\Gamma$ and $\gamma$ in Eq. (7) are the gamma and the lower incomplete gamma
functions, respectively (see Supplementary Information for the derivation of
the model).

Exciton transport in disordered materials is well known to have diffusive to
sub-diffusive character [73]. To accurately infer the anomaly $\alpha$ of the
considered diffusion processes we also simulated dissociation-free exciton
transport. This was done by removing the dissociation rates from the KMC
routine, and sampling 5000 trajectories until $t>2\tau_{\text{diss}}$, with
$\tau_{\text{diss}}^{-1}=k_{r}+k_{IC}+k_{ISC}$ being the theoretical exciton
lifetime of the considered dye (or the largest dissociation lifetime for the
case of mixed dyes). The mean square displacement obtained from these
ensembles was used to fit the anomaly $\alpha$ using the following
dissociation-free model

$\langle r^{2}(t)\rangle\propto t^{\alpha}.$ (8)

The model of Eq. (7) was then used to fit the remaining parameters $D$ and
$k$, by fixing the anomaly $\alpha=\alpha_{\texttt{Fit}}$ gathered from the
dissociation-free trajectories. For single-dopant substrates, the dissociation
rate is expected to be equal to $k_{r}+k_{IC}+k_{ISC}$. This allows us to use
Eq. (7) as a consistency check of the transport model, since one has to expect
$k_{\texttt{Fit}}\approx k_{r}+k_{IC}+k_{ISC}$. The mean exciton lifetime
$\langle\tau\rangle$ was also calculated from the KMC trajectories,

$\langle\tau\rangle=\frac{1}{N}\sum_{i=1}^{N}\tau_{i},$ (9)

with $\tau_{i}$ being the time at which dissociation occurred for the $i$-th
trajectory. For single-dopant substrates $\langle\tau\rangle$,
$\tau_{\text{diss}}$ and $k^{-1}_{\texttt{Fit}}$ are all expected to be in
good agreement with each other, providing a consistency condition for the
transport simulation and model.

In order to gather information on the anisotropy of the transport process,
each cartesian component of the square displacement for the dissociation-free
trajectories was studied to obtain $D_{x}$, $D_{y}$ and $D_{z}$, where
$D=D_{x}+D_{y}+D_{z}$, $\langle r^{2}(t)\rangle=\langle
x^{2}(t)\rangle+\langle y^{2}(t)\rangle+\langle z^{2}(t)\rangle$ and where

$\langle
x^{2}(t)\rangle=\frac{1}{N}\sum_{i=1}^{N}|x_{i}(t)-x_{i}(t_{0})|^{2},$ (10)

and equivalently for $y$ and $z$. The anisotropy profile, shown in Table 3,
was then determined comparing the relative magnitude of each component with
respect to their sum, $D_{j}/D$ for $j=x,y,z$.

Dissociation pathways, sites and dyes were also analysed. For single-dopant
substrate the frequency of each dissociation pathway ($r$, $IC$ or $ISC$) is
expected to be in agreement with its respective dissociation rate ($k_{r}$,
$k_{IC}$ or $k_{ISC}$). For multi-dopant substrates instead, KMC simulation
results provide insight on the interplay between different chromophore
species, valuable to guide the design of efficient LSCs.

####  IV.2.1 Disordered configuration

The exciton transport properties for disordered configurations are presented
in Table 3 and Figures 4 and 5. In disordered substrates transport is
isotropic, as expected, due to the random relative orientation of
chromophores. Single-dopant PDI substrates display the best transport
properties, with an average diffusion length $\langle
r\rangle=4.76\;\text{nm}$, near-unit anomaly $\alpha=0.94$, and diffusion
coefficient $D=2.82\;\text{nm}^{2}/\text{ns}$, despite having the shortest
exciton lifetime $\langle\tau\rangle=6.51\;\text{ns}$. As expected from the
PLQY, PDI’s dissociation is predominately radiative (99.34%), with marginal
contributions from IC (0.65%) and ISC (0.02%) recombination pathways; a
desirable property for the realisation of LSCs. The mean exciton lifetime
$\langle\tau\rangle=6.51\;\text{ns}$ and that obtained by fitting to the model
of Eq. (7) $\tau_{\texttt{fit}}=k_{\texttt{fit}}^{-1}=6.72\>\text{ns}$ are in
good agreement with PDI’s theoretical lifetime
$\tau_{\text{diss}}=k_{\text{diss}}^{-1}=6.60\;\text{ns}$, confirming the
validity of the model.

Transport in Tetracene is sub-diffusive, with $\alpha=0.84$, and diffusion
coefficient $D=0.02\;\text{nm}^{2}/\text{ns}$. The average diffusion length
$\langle r\rangle=0.46\;\text{nm}$ is about one order of magnitude smaller
than PDI, despite the longer exciton lifetime
$\langle\tau\rangle=11.89\;\text{ns}$. For comparison, singlet-exciton
diffusion in Tetracene crystals is around 12 nm [21]. As opposed to PDI,
Tetracene dissociation is dominated by IC (42.82%) and ISC (31.06%), while
radiative recombination only accounts for 26.12% of the dissociation events,
establishing it as a poor choice for single-dopant substrate LSCs. The mean
exciton lifetime $\langle\tau\rangle=11.89\;\text{ns}$ and that obtained by
fitting to the model of Eq. (7)
$\tau_{\texttt{fit}}=k_{\texttt{fit}}^{-1}=11.62\>\text{ns}$ are in good
agreement with Tetracene’s theoretical lifetime
$\tau_{\text{diss}}=k_{\text{diss}}^{-1}=11.74\;\text{ns}$, confirming the
validity of the model.

Figure 6: Dissociation sites in single and multi-dopant substrates with
disordered configurations for samples of 5000 exciton trajectories. All
exciton trajectories start in the centre of the simulation box. Remarkably,
exciton transport properties of PDI-Tetracene mixtures are almost as good as
in PDI, and vastly better than in Tetracene. As a result, mixed-dye materials
provide a PDI-Tetracene singlet-bus mechanism which allows excitons to reach
Tetracene sites that would otherwise be out of range in single-dopant
substrates (see the Supplementary Information for the case of dipole-aligned
configurations).

The transport properties of the multi-dopant PDI+Tetracene substrate are
remarkably unbalanced towards PDI’s properties. This is due to the fact that
the excitons spend the majority of their lifetime ($\langle\tau\rangle\approx
7.65\;\text{ns}$) on PDI chromophores, as expected from the strong asymmetry
in the inter-species FRET rates shown in Table 2. Transport is near-diffusive
with $\alpha=0.94$, diffusion coefficient $D=0.89\;\text{nm}^{2}/\text{ns}$,
and average diffusion length $\langle r\rangle=2.78\;\text{nm}$. Radiative
dissociation (83.14%) is remarkably more frequent than IC (9.82%) and ISC
(7.06%), returning an effective PLQY of around 0.83. While the PDI+Tetracene
transport properties do not seem to favor a multi-dopant substrates over
single-dopant ones, they suggest that mixed-dyes LSCs could prove effective
for harnessing singlet excitons, from PDI, and up-converted triplets, from
Tetracene, while avoiding Tetracene-induced singlet-quenching. The mean
exciton lifetime $\langle\tau\rangle=7.65\;\text{ns}$ is in excellent
agreement with that obtained by fitting to the model of Eq. (7)
$\tau_{\texttt{fit}}=k_{\texttt{fit}}^{-1}=7.65\>\text{ns}$, confirming the
validity of the model.

####  IV.2.2 Dipole-aligned configuration

The exciton transport properties for dipole-aligned configurations are
summarised in Table 3. In dipole-aligned substrates transport is mildly
anisotropic. Perfect anisotropy, characteristic of purely crystalline
materials, is not to be expected, as the orientation factor $\kappa$ depends
on donor and acceptor dipoles and their orientation with respect to the inter-
chromophore separation. Here, instead, the combination of aligned dipoles and
randomly distributed sites returns an intermediate regime of anisotropy. For
both single-dopant and multi-dopant system the anisotropy is unbalanced
towards the $x$ direction, i.e., the typical dipole orientation, with no
preferential contribution between the $y$ and $z$ components.

Exciton diffusion is noticeably enhanced by the alignment of the dipoles, with
improved transport properties across all considered substrates. As for the
case of disordered configurations, single-dopant PDI substrates display the
best transport properties, with an average diffusion length $\langle
r\rangle=5.23\;\text{nm}$, near-unit anomaly $\alpha=0.95$, and diffusion
coefficient $D=3.43\;\text{nm}^{2}/\text{ns}$, despite having the shortest
exciton lifetime $\langle\tau\rangle=6.66\;\text{ns}$. PDI’s dissociation
pathways are unchanged with respect to the case of disordered configurations.
Numerical, fitted and theoretical exciton lifetimes are in good agreement with
each other, with $\langle\tau\rangle=6.66\;\text{ns}$,
$\tau_{\texttt{fit}}=k_{\texttt{fit}}^{-1}=6.69\>\text{ns}$ and
$\tau_{\text{diss}}=k_{\text{diss}}^{-1}=6.60\;\text{ns}$.

Transport in Tetracene is sub-diffusive, with $\alpha=0.81$, and diffusion
coefficient $D=0.03\;\text{nm}^{2}/\text{ns}$. The average diffusion length
$\langle r\rangle=0.5\;\text{nm}$ is only marginally improved by the alignment
of the dipoles, and the dissociation pathways are unaffected. Numerical,
fitted and theoretical exciton lifetimes are in good agreement with each
other, with $\langle\tau\rangle=11.89\;\text{ns}$,
$\tau_{\texttt{fit}}=k_{\texttt{fit}}^{-1}=11.88\>\text{ns}$ and
$\tau_{\text{diss}}=k_{\text{diss}}^{-1}=11.74\;\text{ns}$.

Similar improvement to the transport properties are also obtained for the case
of PDI+Tetracene substrates. Transport is near-diffusive with $\alpha=0.97$,
diffusion coefficient $D=1.10\;\text{nm}^{2}/\text{ns}$, and average diffusion
length $\langle r\rangle=3.15\;\text{nm}$. Radiative dissociation (85.76%)
remains dominant over IC (8.49%) and ISC (5.74%), returning an effective PLQY
of around 0.85. Numerical and fitted exciton lifetimes are in excellent
agreement with each other, with $\langle\tau\rangle=7.56\;\text{ns}$,
$\tau_{\texttt{fit}}=k_{\texttt{fit}}^{-1}=7.56\>\text{ns}$.

It should be noted that the direction of photon emission during a fluorescence
event is directly related to the direction of transition dipole alignment[74].
While in the case of a truly disordered system, this propagation is random, in
the case of a dipole-aligned system, emitted photons can all be considered to
propagate in the same direction. However, aligned-dipole systems have been
shown to increase the probability of molecular aggregation[74], resulting in
an increased chance of excitonic quenching pathways[13]. As such, some work
into the molecular dynamics of dipole-aligned heterogeneous systems may need
to be considered, in addition to bulkier derivative which mitigate aggregative
tendencies.

These results highlight the possible advantages of using multi-dopant
substrates. Here, we have shown that for heterogeneous PDI-Tetracene
substrates, transport is mostly mediated by a particular species of
chromophore (PDI). Although transport tends to terminate upon such dye, it is
possible for it to also end on Tetracene chromophores. More importantly, it
can end on Tetracene chromophores that are much further away from the original
excitation site than for single-dopant Tetracene materials, as shown in Figure
6. As a result, mixed-dye materials provide a PDI-Tetracene singlet-bus
mechanism which allows excitons to reach Tetracene sites that would otherwise
be out of range in single-dopant Tetracene substrates.

##  V Conclusion

Quantum chemical analysis revealed very strong emission properties for PDI, in
contrast to weaker emission properties for Tetracene. Simulated homogeneous
substrates were observed to display strong spectral overlaps, resulting in
diffusion properties which agreed well with the literature. Heterogeneous PDI-
Tetracene systems on the other hand, displayed a favouritism between diffusion
sites. Here, singlet excitons residing on Tetracene sites were observed to
have higher probabilities to hop to PDI sites.

The KMC simulations provided a fast, computationally inexpensive and reliable
approach to simulate the transport properties of singlet excitons in single
and multi-dopant substrates. While some exciton transport properties of
single-dopant substrates can be predicted analytically, multi-dye materials
often necessitate a numerical treatment. The results presented in this work
have been entirely obtained from ab-initio calculations and did not require
any experimental nor phenomenological parametrisation. This establishes our
approach as a powerful technique for LSCs design with existing and yet-to-be
synthesised dyes. Furthermore, the approach used in this work can be extended
to other exciton pathways, assuming that the corresponding rates are known or
can be calculated.

While here we have only considered singlet exciton transport, the inclusion of
triplet excitons dynamics highlights another benefit of multi-dopant
substrates, i.e., the harvesting of a broader portion of the available
spectrum via upconversion pathways. This can be achieved by engineering
substrates containing chromophores with efficient singlet transport ($A$), and
dyes with efficient triplet transport and up-conversion mechanism ($B$), such
that the singlet-bus can be utilised to maximize device efficiency. The dyes
would be required to favour FRET transport of singlets from $B$ to $A$, and
Dexter transport of triplets from $A$ to $B$. In such way, singlet excitons
would mostly undergo FRET transport in $A$, while triplet excitons can be
recycled into higher-energy singlets instead of being lost to non-radiative
pathways. The approach used in this work can be adapted to address the design
of these materials, and we expect it to be used in future studies. Our outlook
to build on this work will involve adding in ray tracing elements, and multi-
exciton interactions, so that we can push towards an even more accurate model.

##  VI Supplementary Information

Supplementary information contains the method used to ascertain a chromophore
concentration analogous to the low concentration scheme, pictorial
representation of the culling radius, and a derivation of the anomalous
diffusion with dissociation rates. We also present the dissociation sites for
dipole-aligned configurations. Finally, the optimised cartesian coordinates of
all relevant electronic states are given.

##  VII Authors’ Contributions

A.M. and F.C. contributed equally to this work.

##  VIII Competing Interests

This work was supported by the Australian Government through the Australian
Research Council (ARC) under the Centre of Excellence scheme (project number
CE170100026). This work was also supported by computational resources provided
by the Australian Government through the National Computational Infrastructure
National Facility and the Pawsay Supercomputer Centre.

##  IX Data Availability

The data that supports the findings of this study are available within the
article [and its Supplementary Information].

##  X References

## References

  * Banal  _et al._ [2016] J. L. Banal, B. Zhang, D. J. Jones, K. P. Ghiggino, and W. W. H. Wong, “Emissive molecular aggregates and energy migration in luminescent solar concentrators,” Acc Chem Res 50, 49–57 (2016).
  * Debije and Verbunt [2011] M. G. Debije and P. P. C. Verbunt, “Thirty years of luminescent solar concentrator research: Solar energy for the built environment,” Advanced Energy Materials 2, 12–35 (2011).
  * McKenna and Evans [2017] B. McKenna and R. C. Evans, “Towards efficient spectral converters through materials design for luminescent solar devices,” Advanced Materials 29, 1606491 (2017).
  * Zhang  _et al._ [2020a] B. Zhang, H. Yang, T. Warner, P. Mulvaney, G. Rosengarten, W. W. H. Wong, and K. P. Ghiggino, “A luminescent solar concentrator ray tracing simulator with a graphical user interface: features and applications,” Methods and Applications in Fluorescence 8, 037001 (2020a).
  * van Eersel, Bobbert, and Coehoorn [2015] H. van Eersel, P. A. Bobbert, and R. Coehoorn, “Kinetic Monte Carlo study of triplet-triplet annihilation in organic phosphorescent emitters,” Journal of Applied Physics 117, 115502 (2015).
  * Mou  _et al._ [2013] W. Mou, S. Hattori, P. Rajak, F. Shimojo, and A. Nakano, “Nanoscopic mechanisms of singlet fission in amorphous molecular solid,” Applied Physics Letters 102, 173301 (2013).
  * Nakano  _et al._ [2019] M. Nakano, K. Okada, T. Nagami, T. Tonami, R. Kishi, and Y. Kitagawa, “Monte carlo wavefunction approach to singlet fission dynamics of molecular aggregates,” Molecules 24, 541 (2019).
  * Ansari-Rad [2021] M. Ansari-Rad, “Transport and annihilation of the triplets in organic phosphorescent systems: Kinetic monte carlo simulation and modeling,” The Journal of Physical Chemistry C 125, 5760–5770 (2021).
  * Xie and Ma [2016a] X. Xie and H. Ma, “Opposite anisotropy effects of singlet and triplet exciton diffusion in tetracene crystal,” ChemistryOpen 5, 201–205 (2016a).
  * Zhang  _et al._ [2017a] B. Zhang, H. Soleimaninejad, D. J. Jones, J. M. White, K. P. Ghiggino, T. A. Smith, and W. W. H. Wong, “Highly fluorescent molecularly insulated perylene diimides: Effect of concentration on photophysical properties,” Chem Mater 29, 8395–8403 (2017a).
  * Tummeltshammer  _et al._ [2016] C. Tummeltshammer, A. Taylor, A. Kenyon, and I. Papakonstantinou, “Losses in luminescent solar concentrators unveiled,” Sol Energ Mat Sol 144, 40–47 (2016).
  * Zhang  _et al._ [2020b] B. Zhang, I. Lyskov, L. J. Wilson, R. P. Sabatini, A. Manian, H. Soleimaninejad, J. M. White, T. A. Smith, G. Lakhwani, D. J. Jones, K. P. Ghiggino, S. P. Russo, and W. W. H. Wong, “FRET-enhanced photoluminescence of perylene diimides by combining molecular aggregation and insulation,” J Mater Chem C 8, 8953–8961 (2020b).
  * Meftahi  _et al._ [2020] N. Meftahi, A. Manian, A. J. Christofferson, I. Lyskov, and S. P. Russo, “A computational exploration of aggregation-induced excitonic quenching mechanisms for perylene diimide chromophores,” J. Chem. Phys. 153, 064108 (2020).
  * Manian  _et al._ [2021] A. Manian, R. A. Shaw, I. Lyskov, W. Wong, and S. P. Russo, “Modeling radiative and non-radiative pathways at both the franck–condon and herzberg–teller approximation level,” The Journal of Chemical Physics 155, 054108 (2021).
  * Leow  _et al._ [2013] S. W. Leow, C. Corrado, M. Osborn, and S. A. Carter, “Monte carlo ray-tracing simulations of luminescent solar concentrators for building integrated photovoltaics,” in  _High and Low Concentrator Systems for Solar Electric Applications VIII_, edited by A. P. Plesniak (SPIE, 2013).
  * Zhang  _et al._ [2017b] F. Zhang, N.-N. Zhang, Y. Zhang, S. Yan, S. Song, B. Jun, and G. Chen, “Theoretical simulation and analysis of large size BMP-LSC by 3d monte carlo ray tracing model,” Chinese Physics B 26, 054201 (2017b).
  * Athanasopoulos  _et al._ [2007] S. Athanasopoulos, J. Kirkpatrick, D. Martínez, J. M. Frost, C. M. Foden, A. B. Walker, and J. Nelson, “Predictive study of charge transport in disordered semiconducting polymers,” Nano Letters 7, 1785–1788 (2007).
  * Kimber  _et al._ [2012] R. G. E. Kimber, E. N. Wright, S. E. J. O’Kane, A. B. Walker, and J. C. Blakesley, “Mesoscopic kinetic Monte Carlo modeling of organic photovoltaic device characteristics,” Physical Review B 86, 235206 (2012).
  * Yost  _et al._ [2012] S. R. Yost, E. Hontz, S. Yeganeh, and T. Van Voorhis, “Triplet vs singlet energy transfer in organic semiconductors: The tortoise and the hare,” Journal of Physical Chemistry C 116, 17369–17377 (2012).
  * Akselrod  _et al._ [2014a] G. M. Akselrod, F. Prins, L. V. Poulikakos, E. M. Lee, M. C. Weidman, A. J. Mork, A. P. Willard, V. Bulović, and W. A. Tisdale, “Subdiffusive exciton transport in quantum dot solids,” Nano Letters 14, 3556–3562 (2014a).
  * Xie and Ma [2016b] X. Xie and H. Ma, “Opposite Anisotropy Effects of Singlet and Triplet Exciton Diffusion in Tetracene Crystal,” ChemistryOpen 5, 201–205 (2016b).
  * Kaiser  _et al._ [2018] W. Kaiser, J. Popp, M. Rinderle, T. Albes, and A. Gagliardi, “Generalized Kinetic Monte Carlo Framework for Organic Electronics,” Algorithms 11, 37 (2018).
  * Ligthart  _et al._ [2018] A. Ligthart, X. de Vries, L. Zhang, M. C. W. M. Pols, P. A. Bobbert, H. van Eersel, and R. Coehoorn, “Effect of Triplet Confinement on Triplet–Triplet Annihilation in Organic Phosphorescent Host–Guest Systems,” Advanced Functional Materials 28, 1804618 (2018).
  * Saxena  _et al._ [2020] R. Saxena, T. Meier, S. Athanasopoulos, H. Bässler, and A. Köhler, “Kinetic Monte Carlo Study of Triplet-Triplet Annihilation in Conjugated Luminescent Materials,” Physical Review Applied 14, 034050 (2020).
  * Marcus [1956] R. A. Marcus, “On the theory of oxidation-reduction reactions involving electron transfer. i,” J Chem Phys 24, 966–978 (1956).
  * Stehr  _et al._ [2011] V. Stehr, J. Pfister, R. F. Fink, B. Engels, and C. Deibel, “First-principles calculations of anisotropic charge-carrier mobilities in organic semiconductor crystals,” Physical Review B 83, 155208 (2011).
  * Settels  _et al._ [2014] V. Settels, A. Schubert, M. Tafipolski, W. Liu, V. Stehr, A. K. Topczak, J. Pflaum, C. Deibel, R. F. Fink, V. Engel, and B. Engels, “Identification of ultrafast relaxation processes as a major reason for inefficient exciton diffusion in perylene-based organic semiconductors,” Journal of the American Chemical Society 136, 9327–9337 (2014).
  * Stehr  _et al._ [2014] V. Stehr, R. F. Fink, B. Engels, J. Pflaum, and C. Deibel, “Singlet exciton diffusion in organic crystals based on marcus transfer rates,” Journal of Chemical Theory and Computation 10, 1242–1255 (2014).
  * McDowall, Johnson, and Patrick [2010] S. McDowall, B. L. Johnson, and D. L. Patrick, “Simulations of luminescent solar concentrators: Effects of polarization and fluorophore alignment,” Journal of Applied Physics 108, 053508 (2010).
  * Şahin, Ilan, and Kelley [2011] D. Şahin, B. Ilan, and D. F. Kelley, “Monte-carlo simulations of light propagation in luminescent solar concentrators based on semiconductor nanoparticles,” Journal of Applied Physics 110, 033108 (2011).
  * Schüler  _et al._ [2008] A. Schüler, A. Kostro, C. Galande, M. V. del Olmo, E. de Chambrier, and B. Huriet, “Principles of monte-carlo ray-tracing simulations of quantum dot solar concentrators,” in  _Proceedings of ISES World Congress 2007 (Vol. I – Vol. V)_ (Springer Berlin Heidelberg, 2008) pp. 1033–1037.
  * Andersen, Panosetti, and Reuter [2019a] M. Andersen, C. Panosetti, and K. Reuter, “A practical guide to surface kinetic monte carlo simulations,” Front. Chem. 7 (2019a), 10.3389/fchem.2019.00202.
  * de Sousa  _et al._ [2018] L. E. de Sousa, D. A. da Silva Filho, R. T. de Sousa, and P. H. de Oliveira Neto, “Exciton diffusion in organic nanofibers: A monte carlo study on the effects of temperature and dimensionality,” Sci. Rep. 8 (2018), 10.1038/s41598-018-32232-5.
  * Suarez, Menger, and Faraji [2020] L. E. A. Suarez, M. F. S. J. Menger, and S. Faraji, “Singlet fission in tetracene: an excited state analysis,” Mol Phys 118, e1769870 (2020).
  * Burdett and Bardeen [2013] J. J. Burdett and C. J. Bardeen, “The dynamics of singlet fission in crystalline tetracene and covalent analogs,” Acc. Chem. Res. 46, 1312–1320 (2013).
  * Burdett  _et al._ [2010] J. J. Burdett, A. M. Müller, D. Gosztola, and C. J. Bardeen, “Excited state dynamics in solid and monomeric tetracene: The roles of superradiance and exciton fission,” The Journal of Chemical Physics 133, 144506 (2010).
  * Brédas  _et al._ [2009] J.-L. Brédas, J. E. Norton, J. Cornil, and V. Coropceanu, “Molecular understanding of organic solar cells: The challenges,” Accounts of Chemical Research 42, 1691–1699 (2009).
  * Balzer  _et al._ [2021] D. Balzer, T. J. A. M. Smolders, D. Blyth, S. N. Hood, and I. Kassal, “Delocalised kinetic monte carlo for simulating delocalisation-enhanced charge and exciton transport in disordered materials,” Chemical Science 12, 2276–2285 (2021).
  * Förster [1948] T. Förster, “Energy migration and fluorescence,” J Biomed Opt 437, 55–75 (1948).
  * Dexter [1953] D. L. Dexter, “A theory of sensitized luminescence in solids,” J Chem Phys 21, 836–850 (1953).
  * Banerjee  _et al._ [2016] S. Banerjee, A. Baiardi, J. Bloino, and V. Barone, “Vibronic effects on rates of excitation energy transfer and their temperature dependence,” J Chem Theory Comput 12, 2357–2365 (2016).
  * Kasha [1950] M. Kasha, “Characterization of electronic transitions in complex molecules,” Discuss Faraday Soc , 14–19 (1950).
  * Voter [2007] A. F. Voter, “Introduction To The Kinetic Monte Carlo Method,” Radiation Effects in Solids , 1–23 (2007).
  * Stamatakis and Vlachos [2011] M. Stamatakis and D. G. Vlachos, “A graph-theoretical kinetic Monte Carlo framework for on-lattice chemical kinetics,” The Journal of Chemical Physics 134, 214115 (2011).
  * Walker, Vogt, and Cole [2015] K. A. Walker, N. Vogt, and J. H. Cole, “Charge filling factors in clean and disordered arrays of tunnel junctions,” Scientific Reports 2015 5:1 5, 1–8 (2015).
  * Ostroverkhova [2016] O. Ostroverkhova, “Organic Optoelectronic Materials: Mechanisms and Applications,” Chemical Reviews 116, 13279–13412 (2016).
  * Oberhofer, Reuter, and Blumberger [2017] H. Oberhofer, K. Reuter, and J. Blumberger, “Charge Transport in Molecular Materials: An Assessment of Computational Methods,” Chemical Reviews 117, 10319–10357 (2017).
  * Lin, Huang, and Hin [2018] F. Lin, J. Huang, and C. Hin, “Electron Transport from Quantum Kinetic Monte Carlo Simulations,” Journal of Physical Chemistry C 122, 20550–20554 (2018), arXiv:1704.07545 .
  * Jørgensen and Grönbeck [2018] M. Jørgensen and H. Grönbeck, “MonteCoffee: A programmable kinetic Monte Carlo framework,” The Journal of Chemical Physics 149, 114101 (2018).
  * Andersen, Panosetti, and Reuter [2019b] M. Andersen, C. Panosetti, and K. Reuter, “A practical guide to surface kinetic Monte Carlo simulations,”  (2019b), arXiv:1904.02561 .
  * Rego and Brandao [2020] A. S. Rego and A. L. Brandao, “General method for speeding up kinetic monte carlo simulations,” Industrial and Engineering Chemistry Research 59, 9034–9042 (2020).
  * Vogt, Jeske, and Cole [2013] N. Vogt, J. Jeske, and J. H. Cole, “Stochastic bloch-redfield theory: Quantum jumps in a solid-state environment,” Physical Review B 88, 174514 (2013).
  * Kranz and Elstner [2016] J. J. Kranz and M. Elstner, “Simulation of Singlet Exciton Diffusion in Bulk Organic Materials,” Journal of Chemical Theory and Computation 12, 4209–4221 (2016).
  * Lyskov, Kleinschmidt, and Marian [2016] I. Lyskov, M. Kleinschmidt, and C. M. Marian, “Redesign of the DFT/MRCI hamiltonian,” J Chem Phys 144, 034104 (2016).
  * Grell and Bradley [1999] M. Grell and D. D. C. Bradley, “Polarized luminescence from oriented molecular materials,” Advanced Materials 11, 895–905 (1999).
  * Lunt  _et al._ [2009] R. R. Lunt, N. C. Giebink, A. A. Belak, J. B. Benziger, and S. R. Forrest, “Exciton diffusion lengths of organic semiconductor thin films measured by spectrally resolved photoluminescence quenching,” Journal of Applied Physics 105, 053711 (2009).
  * Becke [1993] A. D. Becke, “Density-functional thermochemistry. III. the role of exact exchange,” J Chem Phys 98, 5648–5652 (1993).
  * Schäfer, Huber, and Ahlrichs [1994] A. Schäfer, C. Huber, and R. Ahlrichs, “Fully optimized contracted gaussian basis sets of triple zeta valence quality for atoms li to kr,” J Chem Phys 100, 5829–5835 (1994).
  * Frisch  _et al._ [2016] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, and D. J. Fox, “Gaussian16 Revision B.01,” (2016), gaussian Inc. Wallingford CT.
  * Becke [1988] A. D. Becke, “Density-functional exchange-energy approximation with correct asymptotic behavior,” Phys Rev A 38, 3098–3100 (1988).
  * of Karlsruhe [2007] U. of Karlsruhe, “TURBOMOLE V7.3 2018, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, TURBOMOLE GmbH, since 2007;,”  (1989-2007).
  * Schäfer  _et al._ [2000] A. Schäfer, A. Klamt, D. Sattel, J. C. W. Lohrenz, and F. Eckert, “COSMO implementation in TURBOMOLE: Extension of an efficient quantum chemical code towards liquid systems,” Phys Chem Chem Phys 2, 2187–2193 (2000).
  * Kleinschmidt, Tatchen, and Marian [2002] M. Kleinschmidt, J. Tatchen, and C. M. Marian, “Spin-orbit coupling of DFT/MRCI wavefunctions: Method, test calculations, and application to thiophene,” J Comput Chem 23, 824–833 (2002).
  * Kleinschmidt and Marian [2005] M. Kleinschmidt and C. M. Marian, “Efficient generation of matrix elements for one-electron spin–orbit operators,” Chem Phys 311, 71–79 (2005).
  * Kleinschmidt, Tatchen, and Marian [2006] M. Kleinschmidt, J. Tatchen, and C. M. Marian, “SPOCK.CI: A multireference spin-orbit configuration interaction method for large molecules,” J Chem Phys 124, 124101 (2006).
  * Cambi  _et al._ [1991] R. Cambi, D. Cappelletti, G. Liuti, and F. Pirani, “Generalized correlations in terms of polarizability for van der waals interaction potential parameter calculations,” The Journal of Chemical Physics 95, 1852–1861 (1991).
  * Sharipov  _et al._ [2014] A. S. Sharipov, B. I. Loukhovitski, C.-J. Tsai, and A. M. Starik, “Theoretical evaluation of diffusion coefficients of (al2o3)n clusters in different bath gases,” The European Physical Journal D 68 (2014), 10.1140/epjd/e2014-40831-2.
  * Zhurko and Zhurko. [2018] G. Zhurko and D. Zhurko., “Chemcraft. version 1.7 (build 132).”  (2018).
  * Yang and Jang [2020] L. Yang and S. J. Jang, “Theoretical investigation of non-förster exciton transfer mechanisms in perylene diimide donor, phenylene bridge, and terrylene diimide acceptor systems,” J Chem Phys 153, 144305 (2020).
  * Duschinsky [1937] F. Duschinsky, “The importance of the electron spectrum in multi atomic molecules concerning the franck-condon principle,” Acta Physicochim. U.R.S.S. 7, 551–566 (1937).
  * Burgdorff, Ehrhardt, and Loehmannsroeben [1991] C. Burgdorff, S. Ehrhardt, and H. G. Loehmannsroeben, “Photophysical properties of tetracene derivatives in solution. 2. halogenated tetracene derivatives,” The Journal of Physical Chemistry 95, 4246–4249 (1991).
  * Lakowitz [2007] J. R. Lakowitz,  _Principles of Fluorescence Spectroscopy_, 3rd ed. (Springer US, 2007).
  * Akselrod  _et al._ [2014b] G. M. Akselrod, P. B. Deotare, N. J. Thompson, J. Lee, W. A. Tisdale, M. A. Baldo, V. M. Menon, and V. Bulovic, “Visualization of exciton transport in ordered and disordered molecular solids,” Nature Communications 5, 1–8 (2014b).
  * Senes  _et al._ [2016] A. Senes, S. C. J. Meskers, W. M. Dijkstra, J. J. van Franeker, S. Altazin, J. S. Wilson, and R. A. J. Janssen, “Transition dipole moment orientation in films of solution processed fluorescent oligomers: investigating the influence of molecular anisotropy,” Journal of Materials Chemistry C 4, 6302–6308 (2016).