# Solving the Tree Containment Problem Using Graph Neural Networks

Arkadiy Dushatskiy<EMAIL_ADDRESS>
Centrum Wiskunde & Informatica
Amsterdam, the Netherlands Esther Julien<EMAIL_ADDRESS>
Delft University of Technology
Delft, the Netherlands Leen Stougie<EMAIL_ADDRESS>
Centrum Wiskunde & Informatica
Vrije Universiteit Amsterdam
Amsterdam, the Netherlands Leo van Iersel<EMAIL_ADDRESS>
Delft University of Technology
Delft, the Netherlands

###### Abstract

Tree containment is a fundamental problem in phylogenetics useful for
verifying a proposed phylogenetic network, representing the evolutionary
history of certain species. Tree containment asks whether the given
phylogenetic tree (for instance, constructed from a DNA fragment showing tree-
like evolution) is contained in the given phylogenetic network. In the general
case, this is an NP-complete problem. We propose to solve it approximately
using Graph Neural Networks. In particular, we propose to combine the given
network and the tree and apply a Graph Neural Network to this network-tree
graph. This way, we achieve the capability of solving the tree containment
instances representing a larger number of species than the instances contained
in the training dataset (i.e., our algorithm has the inductive learning
ability). Our algorithm demonstrates an accuracy of over $95\%$ in solving the
tree containment problem on instances with up to 100 leaves.

##  1 Introduction

Evolutionary trajectories of biological species can be described using
directed acyclic graphs called _phylogenetic networks_ , which allow the
representation of common evolutionary events such as lateral gene transfer and
hybridization (Bapteste et al., 2013). Leaves of phylogenetic networks are
labeled with the names of studied species. In many cases, some parts of the
DNA structure evolve in graphs that are trees (Janssen & Murakami, 2021). A
valid phylogenetic network representing the evolution of certain species
should be consistent with the previously built trees, namely, it should
contain the trees which represent the evolutionary history of some pieces of
DNA of the considered species. Consequently, a mathematical problem arises
called tree containment: the task is to determine whether a given phylogenetic
tree is contained by a given phylogenetic network.

In general, the tree containment problem is _NP-complete_ (Kanj et al., 2008).
However, for some classes of networks (with certain topological restrictions),
polynomial-time algorithms (in some cases, even linear) were proposed, e.g.,
by Janssen & Murakami (2020), Gambette et al. (2018), and Janssen & Murakami
(2021). An algorithm that solves the tree containment problem should make use
of the labels of the leaves of both the given phylogenetic network and the
tree, as the containment property depends on them (the formal definition is
provided in Section 3.2). This naturally poses a challenge in using machine
learning approaches to solving it. It is especially difficult to handle
situations when the number of leaves in the test instances can be larger than
in the training ones. In this paper, we propose to tackle this challenge and
develop an approach to approximately solving tree containment, including
instances for which no polynomial-time algorithms are known. The main
contributions of this paper can be formulated as follows:

  1. 1.

We propose an approximate solution to the NP-complete tree containment problem
using Graph Neural Networks. To the best of our knowledge, this is the first
time a machine learning approach has been proposed to address this problem.
Our approach can be applied to evaluate the quality of a constructed
phylogenetic network: it can be used to efficiently and with high accuracy
estimate whether trees (e.g., gene trees) are contained in the built network.

  2. 2.

Our proposed approach demonstrates a generalization ability: when trained on
smaller instances (phylogenetic networks and trees with a smaller number of
studied species), it achieves high accuracy (on average, over $95\%$) on
larger instances not included in the training dataset. This inductive learning
ability is achieved through a specific technique: combining the given
phylogenetic network and tree into a single graph while respecting the labels
of the leaves, which is crucial for solving the tree containment problem.

  3. 3.

This work shows the potential of using GNNs for solving phylogenetic problems.
Potentially, ideas from the proposed approach to solving tree containment,
might be applied to solving other phylogenetic problems as well.

##  2 Related work

### Graph Neural Networks

Graph Neural Networks (GNNs) are a broad class of deep neural networks, which
can be applied to solve graph problems (e.g., node and edge classification,
graph classification) and can learn meaningful representations of nodes,
edges, and entire graphs. Message Passing Neural Networks (MPNN) (Gilmer et
al., 2017) is a general framework of neural networks for graphs. In each layer
of an MPNN, each node receives messages containing information about its
neighbors, then those messages are aggregated and combined with the node
information itself. Examples of popular algorithms which use the message-
passing paradigm are Graph Convolutional Networks (GCN) (Kipf & Welling,
2017), Graph Attention Networks (GAT) (Veličković et al., 2018), GraphSAGE
(Hamilton et al., 2017), and GIN (Xu et al., 2019). The comprehensive overview
of deep learning methods for graphs, including the recently proposed and non-
MPNN ones can be found in Khoshraftar & An (2024).

### GNNs for directed graphs

The phylogenetic networks and trees in the tree containment problem are
naturally directed graphs. Most of the GNNs work with undirected graphs.
However, many real-world graphs are directed, and therefore several approaches
have been proposed to adapt GNNs to directed graphs, e.g., (Zhang et al.,
2021b). Recently a more general approach to using GNNs for directed graphs was
proposed in Rossi et al. (2023). This approach called _Dir-GNN_ performs
message passing in two directions (basically transforming the graph to an
undirected one), but considers each direction of message passing separately
with its own learnable weights for each. This ultimately results in GNNs with
more expressive power (Rossi et al., 2023). Dir-GNN was shown to be a general
approach applicable to different GNNs such as GCN (Kipf & Welling, 2017), GAT
(Veličković et al., 2018), and GraphSAGE (Hamilton et al., 2017).

### The tree containment problem

Different approaches have been proposed to solve tree containment in
polynomial time for some specific classes of networks (that have certain
restrictions on the structure). Examples of such classes are binary near-
stable networks (Gambette et al., 2018) and binary tree-child networks
(Gunawan, 2018). In Janssen & Murakami (2020), a more general version of the
tree containment problem was considered (a network in network containment
instead of a tree in network containment), and a linear-time algorithm was
proposed to solve this problem for semi-binary tree-child networks. In the
follow-up work (Janssen & Murakami, 2021), a linear algorithm for tree
(network) containment was also proposed for the so-called cherry picking
networks and rigorous proofs of algorithm complexity were provided. A
different approach was proposed in van Iersel et al. (2023), and Huijsman
(2023): the main idea behind it is to design an algorithm that has an
exponential time complexity in terms of some input parameter, which results in
a polynomial complexity when this parameter is constant (fixed-parameter
tractability). Specifically, the treewidth (measuring the tree-likeness of a
network) was considered as such a parameter.

### Solving hard problems with Graph Neural Networks

GNNs have been applied to solving challenging, in many cases, NP-complete and
NP-hard combinatorial problems such as Traveling Salesman, Maximum Cut,
Minimum Vertex Cover (Barrett et al., 2020; Cappart et al., 2023; Li et al.,
2018; Kool et al., 2019). Usually, in such approaches, it is proposed to use
GNNs as part of the reinforcement learning algorithms: they can extract useful
features from a problem instance and a partial solution, and help in building
a high-quality solution further.

GNNs were also applied to solving a similar type of combinatorial problem to
the considered tree containment problem: graph matching and (sub-)Graph Edit
Distance (GED) calculation problems. There is a rather broad collection of
works in which GNNs are applied to (sub)graph matching problems. One common
approach (Bai et al., 2019; 2020; Lou et al., 2020) is to learn meaningful
embeddings of two graphs and then use them as inputs to the predictor, which
outputs the distance between the graphs. An alternative approach that is used,
for instance, in Doan et al. (2021); Bai et al. (2020), is to have an explicit
algorithm of matching the nodes of one graph with the corresponding nodes of
another (using the learned node-wise embeddings). In Zhang et al. (2021c) it
was proposed to use the hypergraph concept to further improve the algorithm
performance. Recently, Ranjan et al. (2022) proposed essentially a simpler
approach to (sub)graph edit distance learning that demonstrated state-of-the-
art performance outperforming previous works such as Lou et al. (2020). In
contrast to previous works, it does not construct an explicit node-to-node
similarity matrix, but instead just translates the distance between embeddings
of two graphs into their distance function. The embeddings of two graphs are
obtained using Siamese GNNs similar to Xu et al. (2019).

We note that the above-mentioned approaches to the graph matching and GED
problems are not directly applicable to tree containment as, in contrast to
the standard graph matching formulation, it has some specific rules on node
mapping, i.e., leaves mapping should respect their labels (more details are
provided in Section 3.2).

To the best of our knowledge, applying GNNs to solving phylogenetic problems
has not been widely explored yet. The existing works in this field (Zhang,
2022; Mimori & Hamada, 2023) mainly focus on the learning node representations
in the scope of phylogenetic inference, and we note that the phylogenetic
problem considered in this work is different from that.

##  3 Preliminaries

###  3.1 Phylogenetic networks

A rooted phylogenetic network constructed for a set of taxa $X$ is a directed
acyclic graph (DAG) with leaves labelled from the set of labels $x\in X$ and
with some specific restrictions on the in- and out-degrees of the nodes: the
network contains a single root (in-degree zero); all other nodes are either a
leaf (in-degree one, out-degree zero), a tree node (in-degree one, out-degree
greater than one), or a reticulation node (in-degree greater than one, out-
degree one). In this paper, we consider binary phylogenetic networks, meaning
that in- and out-degrees of all nodes are at most two. A phylogenetic tree is
a special case of a phylogenetic network: it is a phylogenetic network that
contains no reticulation nodes.

(a) A phylogenetic tree.

(b) A phylogenetic network that contains the given tree (a).

(c) A phylogenetic network that does not contain the given tree (a). Note: the
leaves’ labels have a different left-to-right order.

Figure 1: Examples of the tree containment problem instances: a tree (a) and
two networks, one (b) does contain the tree, and one (c) does not. Leaves are
labeled with letters (representing a set of taxa). The mapping between the
edges of the tree (a) and the paths of the network (b) (as explained in
Section 3) is depicted with the dashed curves with of corresponding colors.
For the second network (c) no such mapping exists.

###  3.2 Tree containment problem

The tree $T$ (in the general case, a network),
$T=(V^{\prime},E^{\prime},X^{\prime})$) is _contained_ (alternative terms
_embedded_ and _displayed_ can also be used) in network $N=(V,E,X)$ if there
exists an injective mapping of the nodes of $T$ to the nodes of $N$:
$V^{\prime}\mapsto V$ and edges of $T$ to node-disjoint paths in $N$, such
that the leaves’ labels are preserved and the mapping of the edges respects
the mapping of the nodes. In other words, $T$ is contained in $N$, if N has a
subgraph $T^{\prime}$, such that it can be obtained by replacing some edges of
$T$ by (directed) paths. Examples of cases where a tree is contained and is
not contained in a network are shown in Figure 1. In this paper, we consider
the problem of determining whether a phylogenetic tree on taxa $X$ is
contained in a binary phylogenetic network on the same taxa $X$: given a tree
$T$ and a network $N$, the task is to give an answer on whether $T$ is
contained in $N$ or not:

Network Containment
---
Input: binary phylogenetic network $N$ and tree $T$, both have the same set of
leaf labels $X$
Question: is $T$ displayed in $N$?

In the general case, when the network does not belong to one of the specific
classes, e.g., tree-child networks (Pons & Batle, 2021) (such networks that
each non-leaf node has a child that is either a tree vertex or a leaf), this
is a NP-complete problem (Kanj et al., 2008).

##  4 Solving the tree containment problem using a GNN

###  4.1 The summary of our approach

The primary challenge when applying a GNN to solve the tree containment
problem involves ensuring its awareness of the labels on the leaves while
retaining the ability to generalize to instances beyond the training set
(inductive learning ability). For example, employing the labels of one-hot
encoded leaves as node features, (as, e.g., in Zhang et al. (2021a)), does not
enable inductive learning ability. Completely ignoring the labels of the
leaves is a naive approach that, generally, can not be expected to work well,
because the answer to the tree containment problem might differ for two pairs
$(N,T)$ and $(N^{\prime},T^{\prime})$ where $N$ and $N^{\prime}$ have
identical topologies but different leaf labeling (similarly for $T$ and
$T^{\prime}$).

Our proposed approach to address the aforementioned challenge of applying GNNs
is as follows: we combine a network and a tree into a larger _network-tree
graph_ , ensuring shared corresponding leaves. This graph is also called the
_display graph_ in the phylogenetic literature (Bryant & Lagergren, 2006; van
Iersel et al., 2023). This process is illustrated in Figure 2 (left side).
With this technique, the labels of the leaves are taken into account in the
network-tree graph construction as they should be considered to correctly
solve the tree containment. However, we do not use the _specific_ leaves’
labels (for instance, as node features), therefore, this approach,
potentially, can work with the labels which have not been seen in the training
data. In other words, it holds the potential to enable inductive learning
ability.

Figure 2: The schematic illustration of our approach (_Combine-GNN_). To solve
a tree containment problem for the given network and tree, we start by
combining them in a single graph such that the leaves are shared, respecting
their labels. Then, a GNN is applied. We note that the input node features do
not contain the leaves’ labels. The GNN consists of multiple layers, then the
obtained embeddings are concatenated, and a node aggregation (graph readout)
is applied. Finally, an MLP is used to produce the final prediction.

After obtaining a network-tree graph, we apply a GNN to it. The architecture
of the used GNN is depicted in Figure 2 (right side). The goal of the GNN is
to extract the necessary information from the graph structure and produce the
answer to the tree containment problem. First, we use a multi-layer message-
passing GNN to obtain meaningful embeddings of the nodes (a comprehensive
survey of various GNNs can be found, for instance, in Hamilton (2020)).

As the considered phylogenetic networks and trees are directed graphs, simply
transforming them into undirected ones before passing to a GNN might lead to
the loss of some information. On the other hand, if the graphs are kept
directed, without special adaptations, message-passing GNNs cannot be
efficiently applied to our problem as, for instance, leaves, would not pass
any messages (their out-degree is zero) to other nodes, and this also leads to
the loss of information. To take into account the directed nature of the
graphs in the tree containment problem, we use the Dir-GNN approach (Rossi et
al., 2023). It allows bi-directional message-passing, but has separate
learnable weights of AGGREGATE operations in each direction. The message-
passing part of the GNN in our approach is therefore as follows:

$\displaystyle\textit{m}_{i,\leftarrow}^{k}\leftarrow\text{AGGREGATE}^{k}_{i,\leftarrow}\left(\\{(x_{j}^{k-1},x_{i}^{k-1}):(j,i)\in
E\\}\right)$
$\displaystyle\textit{m}_{i,\rightarrow}^{k}\leftarrow\text{AGGREGATE}^{k}_{i,\rightarrow}\left(\\{(x_{j}^{k-1},x_{i}^{k-1}):(i,j)\in
E\\}\right)$ $\displaystyle
x_{i}^{k}\leftarrow\text{COMBINE}^{k}(x_{i}^{k-1},\textit{m}_{i,\leftarrow}^{k},\textit{m}_{i,\rightarrow}^{k})$

where $x_{i}^{k}$ is a node embedding (associated feature vector) of node $i$
at layer $k$, $m_{i,\leftarrow}^{k}$ and $m_{i,\rightarrow}^{k}$ are passed
messages in two directions. AGGREGATE and COMBINE are crucial operations with
usually learnable parameters (weights), and their specific implementations
vary in different GNNs. For instance, in GCN (Kipf & Welling, 2017), they are
implemented as averaging with a consequent linear operation, usually followed
by applying a non-linear function (e.g., ReLU).

After obtaining node-wise embeddings, we gather and concatenate the embeddings
from different layers of the network, similar to the Jumping Knowledge GNN (Xu
et al., 2018) and Ranjan et al. (2022):

$\displaystyle
x_{i}^{final}\leftarrow\text{CONCATENATE}([x_{i}^{0},x_{i}^{1}\dots,x_{i}^{L}]),$

where $L$ is the number of layers, $x_{i}^{0}$ are the features passed to the
first layer. Here, we pass the normalized original node features of the graph,
followed by a linear projection (this projection might, in principle, be
omitted). With multiple GNN layers, we obtain a multi-scale representation of
the graph nodes. Subsequently, these concatenated embeddings are utilized to
obtain the whole graph representation using one of the readout operations
(e.g., an element-wise average pooling). Finally, a Multilayer Perceptron
(MLP) is employed to make the final prediction:

$\displaystyle
G\leftarrow\text{READOUT}(\\{x_{1}^{final},\dots,x_{N}^{final}\\})$
$\displaystyle\textit{output}\leftarrow\text{MLP}(G)$

We call the proposed approach to solving the tree containment problem
_Combine-GNN_. This name reflects the two main ideas behind it: combining a
network and a tree in a single graph, and applying a GNN to it. The code and
used tree containment instances are available at
https://github.com/ArkadiyD/PhyloGNN.

###  4.2 Node features

The considered phylogenetic graphs naturally do not have node features, except
for the leaves (they have the labels). However, we generate some features as a
preprocessing step in order to provide additional helpful information to a GNN
for solving the tree containment problem. We use two types of features for the
nodes of the obtained network-tree graphs: 1) In- and out-degrees of the
nodes; 2) Node origin: for each node of the network-tree graph, we specify
whether the node was a leaf, was part of the network or of the tree. This
information is provided in the form of one-hot encoding. Without explicitly
providing this information, the GNN can try to derive it implicitly, but this
might complicate the learning process.

###  4.3 Time complexity

As analyzed, for instance, in Kipf & Welling (2017), the time complexity of a
GNN consisting of multiple GCN layers is $\mathcal{O}(|E|)$ (assuming the
number of layers and node features are constants). In our case, this means
$\mathcal{O}(|E_{network}|+|E_{tree}|)$, which is equal to
$\mathcal{O}(|E_{network}|+|E_{network}|)=\mathcal{O}(|E_{network}|)$ because
if a tree is larger than the network than it is definitely not contained in
it, and we do not have to run the algorithm to determine the containment.
$\mathcal{O}(|E_{network}|)$ is, in turn, equal to $\mathcal{O}(N)$ ($N$ is
the number of nodes in the network) due to the sparse nature of binary
phylogenetic networks. The usage of Dir-GNN and data preprocessing (obtaining
the network-tree graph) do not change the time complexity of Combine-GNN. The
data processing step (including the network-tree graph creation) has
$\mathcal{O}(N)$ complexity: creating node features (specifically, indicators
of the node origin, see Section 4.2) has linear complexity as the features
need to be created for each node.

##  5 Experiments

###  5.1 Data

Our datasets consist of network-tree pairs (_N_ , _T_) along with
corresponding labels indicating tree containment (label "one" indicating the
network contains the tree, label "zero" indicating no containment).

#### Synthetic data

The data generation process is similar to the one used in Bernardini et al.
(2023). First, we generate a network using a Lateral Gene Transfer (LGT)
generator from Pons et al. (2019) and randomly select a tree subgraph of it.
This way, we obtain a positive example (the tree is contained in the network).
To obtain a negative example (a tree is not contained in the network), we
perform a number of so-called _tail moves_ (Janssen, 2021). A tail move is
defined on two edges $e=(u,v)$ and $f$: it simply moves the node $u$ to the
edge $f$, with respect to the corresponding edges. A valid tail move
guarantees that all necessary properties of the phylogenetic network (listed
in Section 3) are kept intact. Furthermore, a tail move does not change the
main graph properties (number of nodes, number of edges, etc.). After
performing a number of tail moves the tree containment is checked using the
exact algorithm called BOTCH (van Iersel et al., 2023; Huijsman, 2023).

Intuitively, determining whether a tree is contained within a network becomes
more challenging when the network under consideration resembles another
network that contains the specified tree. Thus, we can adjust the difficulty
of the problem by varying the number of tail moves performed. A higher number
of tail moves means greater modifications to the network, resulting in a
network more dissimilar to the original, which contains the given tree. We
create rather challenging instances by using the number of tail moves sampled
uniformly at random from the integers in the range $[1;5]$.

We vary the number of leaves in the generated datasets. We consider two
different scenarios: 1) _Transductive learning_ : the number of leaves in the
validation/test instances does not exceed the number of leaves in the training
data; 2) _Inductive learning_ : the number of leaves in the validation/test
instances is in the interval $(L_{train};2L_{train}]$, while the number of
leaves in the training set is in $[L_{min},L_{train}]$. The minimum number of
leaves $L_{min}$ in all instances is 5, and the maximum number considered is
100. We sample the parameter $\alpha$ of the LGT generator uniformly randomly
from $[0.1;0.2]$. The summary statistics of the generated graphs are provided
in Appendix, Figure 10.

In each setup (defined by the number of leaves in the training and
validation/test instances), we have 10000 training samples (network-tree
pairs), 1000 samples used as validation dataset, and 1000 for the test
dataset.

#### Real-world data

For the the real-world experiments, we use the dataset of gene trees from
bacterial and archaeal genomes (proposed in Beiko (2011), binarized in van
Iersel et al. (2022)). For each set of input trees, we construct a network
with the minimized number of reticulations such that it contains all given
trees using the TrivialRand algorithm from Bernardini et al. (2023). Then, to
obtain negative examples, we perform tail moves on the network similar to the
synthetic datasets construction. As obtaining the correct tree containment
labels using BOTCH is becoming extremely time consuming for large networks
(that have large number of leaves and/or high number of reticulations), we
limit the considered instances to maximum 50 leaves and 5 input trees. As we
have limited number of real-world instances (less than 300 network-tree pairs
for each value of $L\in\\{10,20,30,40,50\\}$, we use them only as validation
and test data while using the synthetic LGT data for training. We also do not
specifically tune any hyperparameters on the real-world data. For this
experiment, we sample the parameter $\alpha$ of the LGT generator uniformly
randomly from $[0.3;0.5]$ (to increase the number of reticulation nodes in the
networks). The number of leaves in the train instances is from $[5;50]$ for
all test instances.

###  5.2 Baselines

The first baseline we consider is a naive approach that disregards the labels
of the leaves. Instead, it attempts to predict the tree containment property
based on manually generated features from both the network and the tree, such
as the average distance from the root to the leaves. The complete list of used
features is listed in Appendix, Table 2. We train an XGBoost model (Chen &
Guestrin, 2016) using these features.

Next, we consider a Siamese GNN (similar to the GREED algorithm (Ranjan et
al., 2022)) applied to the tree containment problem. We make it aware of the
labels of the leaves by using one-hot encoding of the labels as the node
features. This approach does not enable the inductive learning ability, but
could potentially work well in the transductive learning setting. Moreover, in
theory, it can use the automatically extracted features from the network and
the tree (in contrast to the manually constructed features used in the first
considered baseline).

###  5.3 GNN architecture and training hyperparameters tuning

We train GNNs using the AdamW optimizer (Loshchilov & Hutter, 2017) with the
cosine annealing learning rate schedule, number of training steps is 5000;
batch size is 200 (similar parameters were used in Rossi et al. (2023)).

We performed a grid search to find the best performing architecture and the
corresponding training hyperparameters for Combine-GNN. The tuned
hyperparameters and their ranges are the following (the best values are
underlined): initial learning rate: {$10^{-4},\underline{10^{-3}},10^{-2}$};
weight decay: {$\underline{0},10^{-4},10^{-3}$}; dropout:
{$0.0,\underline{0.2}$}; number of GNN layers: {$3,\underline{5}$}, size of
node embeddings: {$32,64,\underline{128}$}; type of GNN operation: {GCN, GIN,
GAT}; readout operation: {mean, max, sum}. Grid search hyperparameter tuning
was also performed for the baselines. The best hyperparameters for them are
provided in Appendix, Table 1. During hyperparameter tuning, we used
validation subsets and random seeds different from those used to report our
main results.

###  5.4 Performance evaluation

We use the balanced accuracy score (Brodersen et al., 2010) as our quality
metric. Although we aimed at creating as balanced datasets as possible (as
shown in Appendix, Figure 10), they are not perfectly balanced (multiple tail
moves could still result in a network that contains the tree, though this is
unlikely). Thus, we use balanced accuracy instead of standard accuracy. In all
experiments, we report the results of test datasets using the epoch at which
the best performance on the validation dataset was achieved.

We compare the performance of Combine-GNN with the aforementioned baselines.
Furthermore, we study the scalability of Combine-GNN in terms of the training
dataset size and increasing size of the test instances (in the inductive
learning setup). Finally, we specifically study how the different design
choices of the GNN used in our algorithm affect the performance. The results
of statistical significance testing are provided in Appendix, Table 3.

We also empirically study the time-wise complexity of Combine-GNN (inference
time) compared to the exact algorithm for the tree containment problem,
namely, BOTCH (van Iersel et al., 2023; Huijsman, 2023). For this analysis, we
specifically consider more challenging instances, excluding tree-child
phylogenetic networks, because for them the tree containment can be solved in
linear time. We note that, for a rigorous time analysis, we adopt a
conservative approach in estimating the time performance of Combine-GNN:
processing one sample at a time (without using batching) and assuming that
both data and the GNN need to be transferred to the GPU each time a sample is
processed. This approach differs from a more practically realistic scenario
where the model is created and transferred to the GPU only once. We use a
system with Intel(R) Xeon(R) Silver 4110 CPU and Nvidia GeForce RTX 2080 Ti
GPUs (for this experiment, only one GPU is utilized).

##  6 Results

Figure 3: Main results: performance of our algorithm and the baselines in
terms of (balanced) accuracy on different datasets; the left part of the graph
shows the results in the inductive setting, the right part of the graph shows
the results in the transductive setting. For each dataset and each algorithm,
we perform five runs with different seeds. Bar height denotes the average
values; error bars denote the $95\%$ confidence interval.

The main results are shown in Figure 3. We see that Combine-GNN shows good
performance in both inductive and transductive learning settings: the average
balanced accuracy is $0.958$ and $0.970$ respectively, in most runs across
datasets of different sizes the accuracy is $>0.95$. Notably, its performance
does not substantially deteriorate on larger instances.

As anticipated, the baseline GNN approach, utilizing one-hot encoding of the
leaves, exhibits better performance in the transductive setting (albeit
substantially lower than Combine-GNN) than in the inductive setting. The
baseline which uses XGBoost on the manually engineered features shows
reasonably good performance, but its underperformance compared to Combine-GNN
indicates that relying solely on general network and tree characteristics does
not suffice to accurately determine tree containment.

#### Performance on much larger instances then used for training

The performance of Combine-GNN on test datasets with instances of increasing
size (including those which are much larger than the instances contained in
the corresponding training dataset) is presented in Figure 4. We observe that
while the performance deteriorates with the difference between size of the
training and test instances increasing, Combine-GNN demonstrates rather robust
behaviour, for instance, when trained on instances with up to 20 leaves, the
balanced accuracy score on datasets with instances with up to 100 leaves,
stays above 0.9. This highlights the generalization ability of the proposed
approach.

Figure 4: Performance (balanced accuracy) for test datasets containing
instances with increasing number of leaves. The leftmost point in each line
shows the performance in the transductive setup, the remaining ones show the
performance in the inductive setup. The number of leaves in the test instances
are sampled randomly uniformly from the specified interval. The results are
averaged over five runs with different seeds.

#### Impact of the training dataset size

Figure 5: Performance (balanced accuracy) for training datasets of different
sizes. The results are averaged over five inductive and five transductive
learning setups with a different number of leaves and five runs with different
seeds for each of them.

The performance of Combine-GNN using training datasets of varying sizes is
presented in Figure 5. For a smaller number of samples, the number of epochs
has been adjusted accordingly, so that the total number of processed samples
during training is equal in all cases ($10^{6}$). We can conclude that a
larger number of samples is clearly beneficial for the performance. It is
plausible that using more than 10000 samples in our training datasets can
result in higher accuracy scores. However, we did not generate larger datasets
due to the computational costs of generating true labels (scaling
exponentially with the instance size).

#### Impact of the number of reticulation nodes

Figure 6: Performance (balanced accuracy) for training and testing datasets
with varying percentage of reticulation nodes. Here, the transductive setup
denotes the test instances with $L$ from $[5;10]$ and inductive denotes the
instances with $L$ from $[11;20]$. The training instances in both setups have
the $L$ from with $[5;10]$. The results are averaged over five runs with
different seeds for each setup and each percentage of reticulation nodes.

The performance of Combine-GNN with training and testing datasets of varying
number of reticulation nodes is presented in Figure 6. We use as training
datasets the datasets with instances with $[5;10]$ leaves and varying number
of reticulations (by controlling the parameter $\alpha$ in the LGT generator).
The test datasets have been generated using the same $\alpha$ correspondingly.
We can conclude that a larger percentage (number) of reticulations makes the
task much more difficult for Combine-GNN to solve. However, we note that in
the transductive setup, the accuracy stays above 0.9 even for the dataset with
>20% reticulation nodes.

#### Real-world data results

Figure 7: Performance (balanced accuracy) for real-world data. The training
data was generated using LGT generator. The results are averaged over five
runs with different seeds for different instance sizes.

The results for real-world data experiments are shown in Figure 7. We note the
challenges posed in this experiment: the training data is a synthetic dataset
generated using LGT generator (real-world data is used only as validation/test
subsets); furthermore, the used hyperparameters were tuned on synthetic
datasets only; the real-world validation subsets are rather small causing
potential overfitting on them (we choose the best epoch based on the
validation performance).

Nevertheless, under all above-mentioned challenges, Combine-GNN demonstrates
reasonably good accuracy, in particular, 0.798 (averaged over five seeds) when
all considered real-world data ($L\in\\{10,20,30,40,50\\}$), is used for
validation and test subsets. It also outperforms the baselines by a
substantial margin. Thus, we can conclude that Combine-GNN can generalize to
the real-world data as well (in contrast to the baselines). Potentially, its
performance might likely be further improved by, for instance, specific
hyperparameter tuning, or using more and/or different training data.

One of the reasons for lower performance of Combine-GNN on the real-world data
is the higher percentage of reticulation nodes ($20\%$ for $L=10$, $17.6\%$
for $L=50$, slighly decreasing with the increase of $L$) than in the datasets
used in our main experiments (Figure 3) and hyperparameter tuning (the
percentage of reticulation nodes there was less than $8\%$). Moreover, it was
shown in Bernardini et al. (2023) that the LGT generator does not generate the
data that has exactly the same properties as this real-world dataset.

###  6.1 Runtime performance

The time performance (inference) results are shown in Figure 8. First, we
observe that the exact algorithm (BOTCH), as anticipated, takes exponential
time to solve tree containment (on instances that are not tree-child
networks). Our algorithm shows polynomial time complexity, namely, the best
linear (in the log-log axes) fit indicates the complexity of
$\mathcal{O}(N^{0.21})$. In particular, we observe that for the considered
graph sizes, the inference time itself scales extremely well (sublinearly).
The data processing step (including network-tree graph creation) has
approximately linear complexity, as expected. The model creation time
demonstrates approximately constant time complexity (the model size does not
depend on the graph size).

Figure 8: Left: scalability of Combine-GNN (along with the best linear fit in
the log-log axes), and the exact algorithm (BOTCH) in terms of the required
time to solve tree containment instances. Right: the breakdown of time
consumption by Combine-GNN. We average all measurements over 10 runs, and use
averaging on the x-axis over 10 bins (the points denote the centers of the
bins).

###  6.2 Ablation studies

The results of all ablation studies are shown in Figure 9.

Figure 9: Effects of different design choices of Combine-GNN to the
performance (balanced accuracy). Performance is averaged over five datasets in
the inductive learning setup and five random seeds for each dataset (25 values
in total for each shown bar). Bar height denotes the average values; error
bars denote the $95\%$ confidence interval.

#### Taking into account the direction of the edges

The use of Dir-GNN demonstrates its advantage over a simpler approach that
treats graphs as undirected. This is a reasonable behavior since retaining
information about the edge directions in the given network and tree is
expected to be helpful in solving the tree containment problem.

#### Node features

We observe that adding the origin of the node (whether it belongs to the
network, the tree, or is a leaf node) to the nodes of the network-tree graph
as features leads to performance improvement. This is an anticipated result,
since these supplementary features explicitly convey to the GNN whether a node
pertains to the network or the tree. We note that adding the node in- and out-
degrees as additional features does not increase the performance further.

#### GNN architecture design choices

We study how different GNN operations (GCN, GAT, GIN) affect the performance.
We can conclude that, while the GCN operation works better than the considered
alternatives, both GAT and GIN operations also demonstrate reasonable
performance. This highlights that our algorithm is not strongly reliant on the
GCN operation and suggests the potential for achieving even better performance
by exploring different GNN operations.

Similar conclusions can be made about the graph-level aggregation operation:
the used maximum pooling operation performs better than the considered
alternatives (specifically, we note that the summation aggregation
underperforms compared to the alternatives).

The embedding size appears to be a rather important hyperparameter: while
there is no substantial difference between performance of embedding size 128
and 64, the performance substantially deteriorates with embedding size 32, and
especially, 16.

##  7 Discussion

In this work, we demonstrated that GNNs can be used as approximate solvers of
the tree containment problem, i.e., they can predict whether a phylogenetic
network contains the given phylogenetic tree or not with high accuracy. While
this capability is valuable, for practical applications, an algorithm that not
only provides containment predictions but also offers a solution
representation might be advantageous. For instance, such an algorithm could
generate a mapping between the nodes of the tree and the network, enabling
solution verification. This verification process might help in filtering out
false positives in determining tree containment: ensuring that a positive
answer consistently reflects the true containment, although not all true
positive cases may be correctly identified as such. We consider this direction
as practically important direction of future research.

Another important direction of future work is exploring whether the proposed
approach (or adaptations of it) can be applied to solving more general
formulations of the tree containment problem, e.g., non-binary phylogenetic
networks and the network (instead of a tree) containment. Such an adaptation
might ease the application of the proposed approach to real-world phylogenetic
datasets.

There are other problems in phylogenetics related to tree containment, for
instance, hybridization. For that problem, the input consists of multiple
trees and the aim is to construct a phylogenetic network with at most $k$
reticulation nodes such that it contains all input trees. Potentially, at
least some of the ideas proposed in this work can be used to solve other
phylogenetic problems such as hybridization, though we acknowledge that this a
non-trivial task and requires additional research.

##  8 Conclusion

In this work, we demonstrated how GNNs can be used to approximately solve an
important NP-complete problem from the phylogenetics field: tree containment.
We proposed an approach named _Combine-GNN_ that allows inductive learning: it
consists of combining the given phylogenetic network and tree in a single
graph and applying a GNN to it. This way, we achieve the inductive learning
ability: it can work for instances with a larger number of leaves (i.e.,
studied species) than the number of leaves of the training instances. Combine-
GNN demonstrates high accuracy in both transductive and inductive learning
setups with instances having up to 100 leaves, clearly outperforming more
simple baseline approaches. Furthermore, we carefully studied how different
design choices of our algorithm affect its performance. Overall, our results
demonstrate that GNNs have grate potential in phylogenetics.

#### Acknowledgements

The work in this paper is supported by the Dutch Research Council (NWO)
through project OCENW.GROOT.2019.015 "Optimization for and with Machine
Learning (OPTIMAL)".

## References

  * Bai et al. (2019) Yunsheng Bai, Hao Ding, Song Bian, Ting Chen, Yizhou Sun, and Wei Wang.  Simgnn: A neural network approach to fast graph similarity computation.  In _Proceedings of the twelfth ACM international Conference on Web Search and Data Mining_ , pp. 384–392, 2019.
  * Bai et al. (2020) Yunsheng Bai, Hao Ding, Ken Gu, Yizhou Sun, and Wei Wang.  Learning-based efficient graph similarity computation via multi-scale convolutional set matching.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34:4, pp. 3219–3226, 2020.
  * Bapteste et al. (2013) Eric Bapteste, Leo van Iersel, Axel Janke, Scot Kelchner, Steven Kelk, James O McInerney, David A Morrison, Luay Nakhleh, Mike Steel, Leen Stougie, and James Whitfield.  Networks: expanding evolutionary thinking.  _Trends in Genetics_ , 29(8):439–441, 2013.
  * Barrett et al. (2020) Thomas Barrett, William Clements, Jakob Foerster, and Alex Lvovsky.  Exploratory combinatorial optimization with reinforcement learning.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34:4, pp. 3243–3250, 2020.
  * Beiko (2011) Robert G Beiko.  Telling the whole story in a 10,000-genome world.  _Biology Direct_ , 6:1–36, 2011.
  * Bernardini et al. (2023) Giulia Bernardini, Leo van Iersel, Esther Julien, and Leen Stougie.  Constructing phylogenetic networks via cherry picking and machine learning.  _Algorithms for Molecular Biology_ , 18(1):13, 2023.
  * Brodersen et al. (2010) Kay Henning Brodersen, Cheng Soon Ong, Klaas Enno Stephan, and Joachim M Buhmann.  The balanced accuracy and its posterior distribution.  In _2010 20th International Conference on Pattern Recognition_ , pp. 3121–3124. IEEE, 2010.
  * Bryant & Lagergren (2006) David Bryant and Jens Lagergren.  Compatibility of unrooted phylogenetic trees is FPT.  _Theoretical Computer Science_ , 351(3):296–302, 2006.
  * Cappart et al. (2023) Quentin Cappart, Didier Chételat, Elias B. Khalil, Andrea Lodi, Christopher Morris, and Petar Veličković.  Combinatorial optimization and reasoning with graph neural networks.  _Journal of Machine Learning Research_ , 24(130):1–61, 2023.
  * Chen & Guestrin (2016) Tianqi Chen and Carlos Guestrin.  XGBoost: A scalable tree boosting system.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pp. 785–794, 2016.
  * Doan et al. (2021) Khoa D. Doan, Saurav Manchanda, Suchismit Mahapatra, and Chandan K. Reddy.  Interpretable graph similarity computation via differentiable optimal alignment of node embeddings.  In _Proceedings of the 44th International ACM SIGIR Conference on Research and Development in Information Retrieval_ , pp. 665–674, 2021.
  * Gambette et al. (2018) Philippe Gambette, Andreas D.M. Gunawan, Anthony Labarre, Stéphane Vialette, and Louxin Zhang.  Solving the tree containment problem in linear time for nearly stable phylogenetic networks.  _Discrete Applied Mathematics_ , 246:62–79, 2018.
  * Gilmer et al. (2017) Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl.  Neural message passing for quantum chemistry.  In _International Conference on Machine Learning_ , pp. 1263–1272. PMLR, 2017.
  * Gunawan (2018) Andreas DM Gunawan.  Solving the tree containment problem for reticulation-visible networks in linear time.  In _Algorithms for Computational Biology: 5th International Conference, AlCoB 2018, Hong Kong, China, June 25–26, 2018, Proceedings 5_ , pp. 24–36. Springer, 2018.
  * Hamilton et al. (2017) Will Hamilton, Zhitao Ying, and Jure Leskovec.  Inductive representation learning on large graphs.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Hamilton (2020) William L Hamilton.  _Graph representation learning_.  Morgan & Claypool Publishers, 2020.
  * Huijsman (2023) Robbert Huijsman.  Treewidth based algorithms for tree containment in phylogenetics.  2023\.  URL http://resolver.tudelft.nl/uuid:3906ebda-d667-4d3e-8bee-3f1f4df78387.
  * Janssen (2021) Remie Janssen.  Heading in the right direction? Using head moves to traverse phylogenetic network space.  _Journal of Graph Algorithms and Applications_ , 25:263–310, 01 2021.
  * Janssen & Murakami (2020) Remie Janssen and Yukihiro Murakami.  Linear time algorithm for tree-child network containment.  In _Algorithms for Computational Biology_ , pp. 93–107, Cham, 2020. Springer International Publishing.
  * Janssen & Murakami (2021) Remie Janssen and Yukihiro Murakami.  On cherry-picking and network containment.  _Theoretical Computer Science_ , 856:121–150, 2021.
  * Kanj et al. (2008) Iyad A Kanj, Luay Nakhleh, Cuong Than, and Ge Xia.  Seeing the trees and their branches in the network is hard.  _Theoretical Computer Science_ , 401(1-3):153–164, 2008.
  * Khoshraftar & An (2024) Shima Khoshraftar and Aijun An.  A survey on graph representation learning methods.  _ACM Transactions on Intelligent Systems and Technology_ , 15(1):1–55, 2024.
  * Kipf & Welling (2017) Thomas N. Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  In _International Conference on Learning Representations_ , 2017.
  * Kool et al. (2019) Wouter Kool, Herke van Hoof, and Max Welling.  Attention, learn to solve routing problems!  In _International Conference on Learning Representations_ , 2019.
  * Li et al. (2018) Zhuwen Li, Qifeng Chen, and Vladlen Koltun.  Combinatorial optimization with graph convolutional networks and guided tree search.  _Advances in Neural Information Processing Systems_ , 31, 2018.
  * Loshchilov & Hutter (2017) Ilya Loshchilov and Frank Hutter.  Decoupled weight decay regularization.  In _International Conference on Learning Representations_ , 2017.
  * Lou et al. (2020) Zhaoyu Lou, Jiaxuan You, Chengtao Wen, Arquimedes Canedo, Jure Leskovec, et al.  Neural subgraph matching.  _arXiv preprint arXiv:2007.03092_ , 2020.
  * Mimori & Hamada (2023) Takahiro Mimori and Michiaki Hamada.  Geophy: Differentiable phylogenetic inference via geometric gradients of tree topologies.  In _ICML 2023 Workshop on Differentiable Almost Everything: Differentiable Relaxations, Algorithms, Operators, and Simulators_ , 2023.
  * Pons et al. (2019) Joan Carles Pons, Celine Scornavacca, and Gabriel Cardona.  Generation of level-$k$ LGT networks.  _IEEE/ACM Transactions on Computational Biology and Bioinformatics_ , 17(1):158–164, 2019.
  * Pons & Batle (2021) Miquel Pons and Josep Batle.  Combinatorial characterization of a certain class of words and a conjectured connection with general subclasses of phylogenetic tree-child networks.  _Scientific reports_ , 11(1):21875, 2021.
  * Ranjan et al. (2022) Rishabh Ranjan, Siddharth Grover, Sourav Medya, Venkatesan Chakaravarthy, Yogish Sabharwal, and Sayan Ranu.  Greed: A neural framework for learning graph distance functions.  _Advances in Neural Information Processing Systems_ , 35:22518–22530, 2022.
  * Rossi et al. (2023) Emanuele Rossi, Bertrand Charpentier, Francesco Di Giovanni, Fabrizio Frasca, Stephan Günnemann, and Michael Bronstein.  Edge directionality improves learning on heterophilic graphs.  _arXiv preprint arXiv:2305.10498_ , 2023.
  * van Iersel et al. (2022) Leo van Iersel, Remie Janssen, Mark Jones, Yukihiro Murakami, and Norbert Zeh.  A practical fixed-parameter algorithm for constructing tree-child networks from multiple binary trees.  _Algorithmica_ , 84(4):917–960, 2022.
  * van Iersel et al. (2023) Leo van Iersel, Mark Jones, and Mathias Weller.  Embedding phylogenetic trees in networks of low treewidth.  _Discrete Mathematics & Theoretical Computer Science_, 25(2), 2023.
  * Veličković et al. (2018) Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio.  Graph attention networks.  In _International Conference on Learning Representations_ , 2018.
  * Xu et al. (2018) Keyulu Xu, Chengtao Li, Yonglong Tian, Tomohiro Sonobe, Ken-ichi Kawarabayashi, and Stefanie Jegelka.  Representation learning on graphs with jumping knowledge networks.  In _International Conference on Machine Learning_ , pp. 5453–5462. PMLR, 2018.
  * Xu et al. (2019) Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka.  How powerful are graph neural networks?  In _International Conference on Learning Representations_ , 2019.
  * Zhang (2022) Cheng Zhang.  Learnable topological features for phylogenetic inference via graph neural networks.  In _International Conference on Learning Representations_ , 2022.
  * Zhang et al. (2021a) Muhan Zhang, Pan Li, Yinglong Xia, Kai Wang, and Long Jin.  Labeling trick: A theory of using graph neural networks for multi-node representation learning.  _Advances in Neural Information Processing Systems_ , 34:9061–9073, 2021a.
  * Zhang et al. (2021b) Xitong Zhang, Yixuan He, Nathan Brugnone, Michael Perlmutter, and Matthew Hirn.  Magnet: A neural network for directed graphs.  _Advances in Neural Information Processing Systems_ , 34:27003–27015, 2021b.
  * Zhang et al. (2021c) Zhen Zhang, Jiajun Bu, Martin Ester, Zhao Li, Chengwei Yao, Zhi Yu, and Can Wang.  H2mn: Graph similarity learning with hierarchical hypergraph matching networks.  In _Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining_, pp. 2274–2284, 2021c.

## Appendix

Figure 10: Number of nodes (upper left), number of edges (upper right), number
of reticulation nodes (bottom left) and label (whether the tree is contained
or not) statistics of the graphs in the used datasets.

GNN-based baseline
---
Hyperparameter | Considered values
initial learning rate |  {$10^{-4},\underline{10^{-3}},10^{-2}$}
weight decay |  {$0,\underline{10^{-4}},10^{-3}$}
dropout |  {$\underline{0.0},0.2$}
number of GNN layers |  {$3,\underline{5}$}
size of node embeddings |  {$32,64,\underline{128}$}
type of GNN operation |  {GCN, GIN, GAT}
readout operation |  {mean, max, sum}

XGBoost-based baseline
---
Hyperparameter | Considered values
estimators |  {$\underline{50},300$}
max depth |  {$3,5,10,\underline{\text{None}}$}
learning rate |  {$10^{-4},10^{-3},10^{-2},\underline{10^{-1}}$}
max number of leaves |  {$0,10,\underline{100}$}
subsample |  {$0.5,\underline{0.9}$}

Table 1: Search spaces of hyperparameters for the considered baselines. The best values are underlined. basic calculated values |  calculated features
---|---
$\\#N,\\#T$: number of nodes in the network and the tree respectively |  $\\#T/\\#N$
$N_{depth}$ \- network depth, $T_{depth}$ \- tree depth |  $N_{depth}/\\#N$
$T_{depth}/\\#T$
$d_{network},d_{tree}:$ arrays of distances to the corresponding leaves in network and tree respectively |  $\min(d_{network})/N_{depth}$
$\text{mean}(d_{network})/N_{depth}$
$\min(d_{tree})/T_{depth}$
$\text{mean}(d_{tree})/T_{depth}$
$\min(d_{network})-\max(d_{tree})$
$(\min(d_{network})-\max(d_{tree}))/\\#N$
$\min(d_{network}-d_{tree})$
$\max(d_{network}-d_{tree})$
$\text{mean}(d_{network}-d_{tree})$
Table 2: The features used in the XGBoost-based baseline.

Experiment |  Hypothesis |  p-value
---|---|---
Main results |  Combine-GNN performs better than GNN baseline |  0.000
Combine-GNN performs better than XGBoost baseline |  0.000
Ablation: GNN operation |  GCN performs better than GIN |  0.011
GCN performs better than GAT |  0.01
Ablation: readout operation |  max performs better than sum |  0.001
max performs better than mean |  0.000
Ablation: embedding size |  embedding size 128 performs better than embedding size 64 |  0.063
embedding size 128 performs better than embedding size 32 |  0.000
embedding size 128 performs better than embedding size 16 |  0.000
Ablation: using graph directionality |  Dir-GNN performs better standard GNN (on undirected graphs) |  0.000
Ablation: using different features |  Using node type features and node degree features performs better than using node type features |  0.336
Using node type features and node degree features performs better than using degree features |  0.000
Using node type features and node degree features performs better than using no features |  0.000
Table 3: The results of statistical significance tests. We use the one-sided Wilcoxon test, testing differences between the balanced accuracy scores. Statistically significant results at $\alpha=0.05$ and applied Bonferroni correction (in each experiment, we apply the correction with $n$ set to the number of tested hypotheses) are highlighted with bold font. The p-values are rounded to three decimal places. Setup | Balanced accuracy | Precision | Recall
---|---|---|---
transductive | 0.970 | 0.958 | 0.986
inductive | 0.957 | 0.942 | 0.976
Table 4: Balanced accuracy along with precision and recall metrics for
Combine-GNN. These numbers are averaged over all inductive and transductive
setups from Figure 3 (main results).


# Concerning the stability of biexcitons in hybrid HJ aggregates of
$\pi$-conjugated polymers

Eric R. Bittner<EMAIL_ADDRESS>Department of Chemistry, University
of Houston, Houston, Texas 77204, United States  Carlos Silva  School of
Chemistry and Biochemistry, Georgia Institute of Technology, 901 Atlantic
Drive, Atlanta, GA 30332, United States  School of Physics, Georgia Institute
of Technology, 837 State Street, Atlanta, GA 30332, United States  School of
Materials Science and Engineering, Georgia Institute of Technology, North
Avenue, Atlanta, GA 30332, United States

###### Abstract

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

##  I Introduction

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

$\displaystyle\ce{S0+2$\hbar\omega$->[2S1]$^{\ddagger}$->2e^{-}+2h^{+}}$ (1)

bimolecular annihilation [4],

$\displaystyle\ce{S1+S1->[2S1]$^{\ddagger}$->S0+S0}$ (2)

and singlet fission producing triplet (T1) states [5]

$\displaystyle\ce{S0+2$\hbar\omega$->[2S1]$^{\ddagger}$->T1+T1}.$ (3)

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

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

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

(a)(b) Figure 1: (a) Chemical structure of the PBTTT polymer. (b) Local
packing of excitonic units (ellipsoids) with superimposed transition dipoles
(yellow arrows) with J-like interactions along the horizontal axis and H-like
interactions along the vertical axis.

##  II Theoretical model

###  II.1 Homogeneous lattice model in 1- and 2-dimensions

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

$\displaystyle
a_{\textbf{n}}a_{\textbf{m}}^{\dagger}-(-1)^{\delta_{\textbf{nm}}}a_{\textbf{m}}^{\dagger}a_{\textbf{n}}=\delta_{\textbf{nm}}.$
(4)

These operators are Paulion operators that create and remove single
excitations on to a site labeled by n. On a given site, they obey the Fermion
relation $a_{\textbf{n}}a_{\textbf{n}}^{\dagger}+a_{n}^{\dagger}a_{n}=1$ which
enforces Pauli exclusion since $a_{n}a_{n}=a_{n}^{\dagger}a_{n}^{\dagger}=1$,
but commute “off-site” with $[a_{\textbf{n}},a_{\textbf{m}}^{\dagger}]=0$ when
$\textbf{n}\neq\textbf{m}$. This is different than the usual Fermion algebra
where the anti-commutation rule is applies over all sites, ultimately giving
rise to the exchange interaction. Further, we can write a generic multi-
exciton Hamiltonian as

$\displaystyle
H=\sum_{\textbf{nm}}h_{\textbf{nm}}a_{\textbf{n}}^{\dagger}a_{\textbf{m}}+\frac{1}{2}\sum_{\textbf{nm}}U_{\textbf{nm}}a^{\dagger}_{\textbf{n}}a^{\dagger}_{\textbf{m}}a_{\textbf{n}}a_{\textbf{m}}.$
(5)

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

For a one-dimensional chain with lattice spacing $a$, $n$ is simply an index
along the chain such that the site location is given by
$r_{\textbf{n}}=\textbf{n}a$. However, for 2- and 3-dimensional systems, we
shall take it as an n-tuple index specifying the site location. For the single
particle term, $h_{\textbf{nn}}$ is the excitation energy for single site
($\varepsilon_{\textbf{n}}$) and $h_{\textbf{nm}}$ (for
$\textbf{n}\neq\textbf{m}$) corresponds to the hopping integral between sites.
Upon transforming into the reciprocal space,

$\displaystyle\tilde{a}_{\textbf{k}}=\frac{1}{\sqrt{N}}\sum_{\textbf{n}}e^{-\textbf{k}\cdot\textbf{n}}a_{\textbf{n}}$
(6)

one finds the single particle energy dispersion as

$\displaystyle(\varepsilon(\textbf{k})-E)e^{i\textbf{k}\cdot\textbf{n}}=\sum_{\textbf{m}}h_{\textbf{nm}}e^{i\textbf{k}\cdot\textbf{m}}.$
(7)

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

$\displaystyle|\psi\rangle=\sum_{\textbf{kk}^{\prime}}c_{\textbf{kk}^{\prime}}\tilde{a}^{\dagger}_{\textbf{k}}\tilde{a}_{\textbf{k}^{\prime}}^{\dagger}.$
(8)

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

Figure 2: Sketch of 2D lattice model with superimposed transition dipole for
each site. In (A), we define a J-aggregate basis and create 2J biexciton
configurations along each row. In (B), we define an H-aggregate basis along
each column. As discussed in the text, these are equivalent representations of
the full 2D problem and are useful for reducing the 2D problem into two
separate (but formally identical) 1D problems.

###  II.2 Local exciton approximation

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

$\displaystyle\hat{J}_{m,k_{x}}=\frac{1}{\sqrt{N}}\sum_{n}e^{ik_{x}n}a_{(n,m)}$
(9)
$\displaystyle\hat{H}_{n,k_{y}}=\frac{1}{\sqrt{N}}\sum_{m}e^{ik_{y}m}a_{(n,m)}.$
(10)

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

Thus, we can write

$\displaystyle\sum_{\textbf{nm}}h_{\textbf{nm}}a^{\dagger}_{\bf n}a_{\bf m}$
$\displaystyle=\sum_{k_{x}mm^{\prime}}(\varepsilon_{J}(k_{x})\delta_{mm^{\prime}}+t_{H}(mm^{\prime})\hat{J}^{\dagger}_{m,k_{x}}\hat{J}_{m^{\prime},k_{x}}$
$\displaystyle=\sum_{k_{y}nn^{\prime}}(\varepsilon_{H}(k_{y})\delta_{nn^{\prime}}+t_{J}(nn^{\prime}))\hat{H}^{\dagger}_{n,k_{y}}\hat{H}_{n^{\prime}k_{y}}$
(11)

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

$\displaystyle\varepsilon_{J}(k_{x})$
$\displaystyle=\epsilon_{0}+2t_{J}\cos(k_{x})$ (12)
$\displaystyle\varepsilon_{H}(k_{y})$
$\displaystyle=\epsilon_{0}+2t_{H}\cos(k_{y})$ (13)

where $t_{J}<0$ is the nearest-neighbor coupling along the $x$-direction and
$t_{H}>0$ is the nearest-neighbor coupling taken in the $y$-direction. To
remind, throughout we are taking the wavevector $k\in[-\pi,\pi]$. Since the
problem is formally separable between $x$ and $y$ directions, the single-
particle terms do not mix wavevector components since $k_{x}$ and $k_{y}$ are
“good” quantum numbers for this system. Using the $\hat{J}$ and $\hat{H}$
operators, we can reduce the 2-particle/2-dimensional problem into a one for
pair of particles within a 1 dimensional frame without loss of generality.

Figure 3:  (a) Energy bands for attractive (2H) and repulsively (2J) bound
biexcitons. Dashed lines correspond to the non-interacting $U\to 0$ limit of
the model. Solid curves are plotted for the case of $U/2t>1$ as indicated. The
dispersion curves need to be taken at the corresponding points in $k$-space
with a corresponding shift in the energy origins for each exciton species.
(b)Variation in the biexciton binding energy at $k=0$ with increasing
exciton/exciton interactions.

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

If we use the $J_{km}$ (or $H_{kn}$) states as a basis for a given value of
$k$, we can use the Greens function approach to find the biexciton energies as

$\displaystyle E_{XX}(k)=2\varepsilon_{x}+\frac{4t^{2}}{U}\cos^{2}(k)+U$ (14)

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

###  II.3 Optical Transitions

It is important to keep in mind that two types of excitons sketched in Fig.1
occur at different wave-vectors. In Fig.1(a) the J-like exciton has its
longest wavelength in the $x$-direction (corresponding to $k_{x}=0$) while
having the shortest possible wavelength in the $y$-direction (i.e.
$k_{y}=\pm\pi/a$). Transitions to this state from the exciton vacuum (ground
states) are allowed under the dipole approximation for the transition moment.
Consequently, double excitations can create either pairs of free J-like
excitons or bound 2J excitons. Similarly, optically allowed H-like single
excitons are delocalized along the $y$-direction but localized along the
$x$-direction (i.e. ${\bf k}=\\{\pm\pi/a,0\\}$) implying that the double
excitations must also occur at these wave-vectors. Bearing this in mind, the
dispersion curves plotted in Fig. 3 need to be taken at the corresponding
points in $k$-space with a corresponding shift in the energy origins for each
exciton species.

###  II.4 Biexciton stability

According to our model, Frenkel biexcitons mix J-like and H-like character in
terms of their collective quantum behavior with the requirement that the ratio
of the exciton/exciton interaction and the perpendicular hopping term be
$U/t>0$ which gives rise to localized biexciton states in the perpendicular
direction. For the 1D $\delta$-function potential, any attractive interaction
with $U/t>0$ produces a localized state with localization length
$\lambda=\kappa^{-1}=2t/U$.111Note, that here $t=-\hbar^{2}/2\mu$ carries
units of $[{\rm Energy}]\cdot[{\rm length}]^{2}$ and $U$ carries units of
$[{\rm Energy}]\cdot[{\rm length}]^{-1}$ so that $\lambda$ carries the
appropriate units of [length]. For the 1D lattice, bound biexciton states
occur outside the band for free biexcitons. To gain further insight into the
stability of these states, we turn to a continuum model and work in a relative
coordinate reference frame where $x=|r_{1}-r_{2}|$ is the separation between
two localized excitons. Thus, the biexciton Schrödinger equation can be
approximated as

$\displaystyle t\psi^{\prime\prime}+U\delta(x)\psi=E\psi$ (15)

where $U$ is the contact interaction between the two excitons. For bound
states, $\psi(x)$ must vanish as $x\to\pm\infty$, giving that

$\displaystyle\psi(x)=\left\\{\begin{array}[]{cc}\sqrt{\kappa}e^{-\kappa
x}&x>0\\\ \sqrt{\kappa}e^{+\kappa x}&x<0\end{array}\right.$ (18)

where $\kappa=U/2t$ is a positive constant and $E=t\kappa^{2}$. In general, we
take $t=-\hbar^{2}/2\mu_{eff}$ and $U<0$ for an attractive potential giving
rise to a bound state energetically below the continuum for the scattering
states.

####  II.4.1 Lattice reorganization in the impurity model

Generally speaking, one cannot discount the role of lattice reorganization and
relaxation when discussing excitons and biexcitons in organic polymer
semiconducting systems. To study this, we append to the 1D impurity model a
continuum model for the medium a term coupling the biexciton to the lattice as
per the Davydov model. [21, 22, 23, 24, 25, 26, 27, 28, 29, 30] The resulting
equations of motion read

$\displaystyle i\hbar\dot{\psi}(x)$
$\displaystyle=(t\nabla^{2}+U\delta(x)+(E_{o}+2\chi\nabla u(x)))\psi(x)$
$\displaystyle\ddot{u}-\frac{k}{m}\nabla^{2}u$
$\displaystyle=2\frac{\chi}{m}\nabla|\psi|^{2}$ (19)

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

$\displaystyle\nabla u=-\frac{2\chi}{k(1-(v/c)^{2})}|\psi|^{2}$ (20)

that gives us a non-linear Schrödinger equation

$\displaystyle i\hbar\dot{\psi}=(t\nabla^{2}+g|\psi|^{2}+U\delta(x))\psi$ (21)

where $g=-4\chi^{2}/(k(1-(v/c)^{2}))$ and $c$ is the sound velocity. Note that
$E_{o}$ is a constant given by

$\displaystyle
E_{o}=E-2t+\frac{1}{2}\int_{-\infty}^{\infty}\left(m\ddot{u}^{2}+ku^{\prime\prime}\right)dx.$
(22)

that we can ignore for purposes of this analysis. The $\delta$-function
potential implies the wave function should have the form in Eq. 18. Taking
$\kappa$ as a variational parameter and minimizing the total energy, one
obtains

$\displaystyle\kappa=\frac{U}{2t}+\frac{g}{8t}.$ (23)

$\kappa>0$ is necessary to produce a localized state and from above $U/t>0$
and $g<0$ from its definition above, we have a stability requirement that if
$U>0$ and $t>0$, then $-g<4U$. Solving for the binding energy,

$\displaystyle E_{B}=\frac{(4U+g)^{2}}{64t}$ (24)

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

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

Figure 4: Biexciton binding energy (in reduced units) for the 2J and 2H cases
versus increasing exciton/lattice coupling. For the attractively bound 2H,
lattice reorganization is expected to stabilize the biexciton state while
destabilizing the 2J. The dashed line indicates the limit of stability for the
2J case.

##  III Comparison to experimental measurement

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

Figure 5: Spectral cuts along the two-quantum coherence energy axis, relative
to the two-quantum diagonal energy axis ($\hbar\omega_{2Q}=2\hbar\omega_{1Q}$,
$E_{XX}$) at fixed $\hbar\omega_{43}=2.06$ (H-aggregate resonance, blue),
$1.99$ (J-aggregate resonance, purple), and $1.94$ eV (J-aggregate resonance,
red). Figure from Ref.19

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

##  IV Perspective

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

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

###### Acknowledgements.

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

Author Declaration: Both authors contributed equally to the design and
analysis of this study. Theoretical methodology was developed by ERB and
experiments were designed and analyzed by CS. The manuscript was written and
edited by ERB. Both authors have read and agreed to the published version of
the manuscript. The authors have no conflicts to disclose.

Data Availability:The data that supports the findings of this study are
available within the article.

## References

  * Frenkel [1931] J. Frenkel, “On the transformation of light into heat in solids. i,” Phys. Rev. 37, 17 (1931).
  * Agranovich and Toshich [1968] V. Agranovich and B. Toshich, “Collective properties of frenkel excitons,” Sov. Phys. JETP 26, 104–112 (1968).
  * Silva  _et al._ [2001] C. Silva, A. S. Dhoot, D. M. Russell, M. A. Stevens, A. C. Arias, J. D. MacKenzie, N. C. Greenham, R. H. Friend, S. Setayesh, and K. Müllen, “Efficient exciton dissociation via two-step photoexcitation in polymeric semiconductors,” Phys. Rev. B 64, 125211 (2001).
  * Stevens  _et al._ [2001] M. A. Stevens, C. Silva, D. M. Russell, and R. H. Friend, “Exciton dissociation mechanisms in the polymeric semiconductors poly (9, 9-dioctylfluorene) and poly (9, 9-dioctylfluorene-co-benzothiadiazole),” Phys. Rev. B 63, 165213 (2001).
  * Silva  _et al._ [2002] C. Silva, D. M. Russell, A. S. Dhoot, L. M. Herz, C. Daniel, N. C. Greenham, A. C. Arias, S. Setayesh, K. Müllen, and R. H. Friend, “Exciton and polaron dynamics in a step-ladder polymeric semiconductor: the influence of interchain order,” J. Phys.: Condens. Matter 14, 9803 (2002).
  * Scholes  _et al._ [2011] G. D. Scholes, G. R. Fleming, A. Olaya-Castro, and R. Van Grondelle, “Lessons from nature about solar light harvesting,” Nat. Chem. 3, 763–774 (2011).
  * Avanaki and Schatz [2021] K. N. Avanaki and G. C. Schatz, “Mechanistic understanding of entanglement and heralding in cascade emitters,” J. Chem. Phys. 154, 024304 (2021).
  * Spano, Agranovich, and Mukamel [1991] F. C. Spano, V. Agranovich, and S. Mukamel, “Biexciton states and two-photon absorption in molecular monolayers,” J. Chem. Phys. 95, 1400–1409 (1991).
  * Vektaris [1994] G. Vektaris, “A new approach to the molecular biexciton theory,” The Journal of Chemical Physics 101, 3031–3040 (1994), https://doi.org/10.1063/1.467616 .
  * Guo, Chandross, and Mazumdar [1995] F. Guo, M. Chandross, and S. Mazumdar, “Stable biexcitons in conjugated polymers,” Phys. Rev. Lett. 74, 2086 (1995).
  * Gallagher and Spano [1996] F. B. Gallagher and F. C. Spano, “Theory of biexcitons in one-dimensional polymers,” Phys. Rev. B 53, 3790 (1996).
  * Mazumdar  _et al._ [1996] S. Mazumdar, F. Guo, K. Meissner, B. Fluegel, and N. Peyghambarian, “Exciton-to-biexciton transition in quasi-one-dimensional organics,” J. Chem. Phys. 104, 9292–9296 (1996).
  * Agranovich  _et al._ [2000] V. Agranovich, O. Dubovsky, D. Basko, G. La Rocca, and F. Bassani, “Kinematic frenkel biexcitons,” J. Lumin. 85, 221–232 (2000).
  * Kun  _et al._ [2009] G. Kun, X. Shi-Jie, L. Yuan, Y. Sun, L. De-Sheng, and Z. Xian, “Effect of interchain coupling on a biexciton in organic polymers,” Chin. Phys. B 18, 2961 (2009).
  * Xiang, Litinskaya, and Krems [2012] P. Xiang, M. Litinskaya, and R. V. Krems, “Tunable exciton interactions in optical lattices with polar molecules,” Phys. Rev. A 85, 061401 (2012).
  * Chakrabarti  _et al._ [1998] A. Chakrabarti, A. Schmidt, V. Valencia, B. Fluegel, S. Mazumdar, N. Armstrong, and N. Peyghambarian, “Evidence for exciton-exciton binding in a molecular aggregate,” Phys. Rev. B 57, R4206 (1998).
  * Klimov  _et al._ [1998] V. I. Klimov, D. McBranch, N. Barashkov, and J. Ferraris, “Biexcitons in $\pi$-conjugated oligomers: Intensity-dependent femtosecond transient-absorption study,” Phys. Rev. B 58, 7654 (1998).
  * McCulloch  _et al._ [2006] I. McCulloch, M. Heeney, C. Bailey, K. Genevicius, I. MacDonald, M. Shkunov, D. Sparrowe, S. Tierney, R. Wagner, W. Zhang, _et al._ , “Liquid-crystalline semiconducting polymers with high charge-carrier mobility,” Nat. Mater. 5, 328–333 (2006).
  * Gutiérrez-Meza  _et al._ [2021] E. Gutiérrez-Meza, R. Malatesta, H. Li, I. Bargigia, A. R. S. Kandada, D. A. Valverde-Chávez, S.-M. Kim, H. Li, N. Stingelin, S. Tretiak, E. R. Bittner, and C. Silva-Acuña, “Frenkel biexcitons in hybrid hj photophysical aggregates,” Science Advances 7, eabi5197 (2021), https://www.science.org/doi/pdf/10.1126/sciadv.abi5197 .
  * Note [1] Note, that here $t=-\hbar^{2}/2\mu$ carries units of $[{\rm Energy}]\cdot[{\rm length}]^{2}$ and $U$ carries units of $[{\rm Energy}]\cdot[{\rm length}]^{-1}$ so that $\lambda$ carries the appropriate units of [length].
  * Scott [1992] A. Scott, “Davydov’s soliton,” Physics Reports 217, 1–67 (1992).
  * Edler  _et al._ [2004] J. Edler, R. Pfister, V. Pouthier, C. Falvo, and P. Hamm, “Direct observation of self-trapped vibrational states in $\alpha$-helices,” Physical Review Letters 93, 106405 (2004).
  * Sun, Luo, and Zhao [2010] J. Sun, B. Luo, and Y. Zhao, “Dynamics of a one-dimensional holstein polaron with the davydov ansätze,” Phys. Rev. B 82, 014305 (2010).
  * Goj and Bittner [2011] A. Goj and E. R. Bittner, “Mixed quantum classical simulations of excitons in peptide helices,” The Journal of Chemical Physics 134, 205103 (2011), https://doi.org/10.1063/1.3592155 .
  * Goj and Bittner [2010] A. Goj and E. R. Bittner, “Mixed quantum classical simulations of vibrational excitations in peptide helices,” in  _Frontiers in Optics 2010/Laser Science XXVI_ (Optical Society of America, 2010) p. LTuC1.
  * Georgiev and Glazebrook [2019a] D. D. Georgiev and J. F. Glazebrook, “On the quantum dynamics of davydov solitons in protein $\alpha$-helices,” Physica A: Statistical Mechanics and its Applications 517, 257–269 (2019a).
  * Georgiev and Glazebrook [2019b] D. D. Georgiev and J. F. Glazebrook, “Quantum tunneling of davydov solitons through massive barriers,” Chaos, Solitons and Fractals 123, 275–293 (2019b).
  * Kerr and Lomdahl [1987] W. C. Kerr and P. S. Lomdahl, “Quantum-mechanical derivation of the equations of motion for davydov solitons,” Phys. Rev. B 35, 3629–3632 (1987).
  * Davydov [1977] A. Davydov, “Solitons and energy transfer along protein molecules,” Journal of Theoretical Biology 66, 379–387 (1977).
  * Davydov [1973] A. Davydov, “The theory of contraction of proteins under their excitation,” Journal of Theoretical Biology 38, 559–569 (1973).