# GT-CausIn: a novel causal-based
insight for traffic prediction

Ting Gao, Rodrigo Kappes Marques, Lei Yu This paper was produced by the IEEE
Publication Technology Group. They are in Piscataway, NJ.Manuscript received
April 19, 2021; revised August 16, 2021.

###### Abstract

Traffic forecasting is an important application of spatiotemporal series
prediction. Among different methods, graph neural networks have achieved so
far the most promising results, learning relations between graph nodes then
becomes a crucial task. However, improvement space is very limited when these
relations are learned in a node-to-node manner. The challenge stems from (1)
obscure temporal dependencies between different stations, (2) difficulties in
defining variables beyond the node level, and (3) no ready-made method to
validate the learned relations. To confront these challenges, we define
legitimate traffic causal variables to discover the causal relation inside the
traffic network, which is carefully checked with statistic tools and case
analysis. We then present a novel model named Graph Spatial-Temporal Network
Based on Causal Insight (GT-CausIn), where prior learned causal information is
integrated with graph diffusion layers and temporal convolutional network
(TCN) layers. Experiments are carried out on two real-world traffic datasets:
PEMS-BAY and METR-LA, which show that GT-CausIn significantly outperforms the
state-of-the-art models on mid-term and long-term prediction.

###### Index Terms:

Traffic forecasting, Causal Discovery, Spatiotemporal Network

††publicationid: pubid: 0000–0000/00$00.00 © 2021 IEEE

##  I Introduction

As an important aspect of smart cities[1], traffic forecasting is gaining more
and more popularity in this data era. Traffic forecasting aims at predicting
future traffic volumes with historical records. Different models have been
proposed and evolved. Knowledge-driven models focus on defining functional
characteristics of traffic elements and simulating various traffic service
demands [2], where analysis is generally theory-based and requires constrained
mathematical assumptions. Data-driven models contain three major categories:
time-series analysis methods (e.g., ARIMA [3], SARIMA [4]), traditional
machine learning methods (e.g., SVM [5], SVR [6]) and deep learning based
methods (e.g., CNN [7, 8]).

Deep learning methods outperform the others by a large margin and have become
the mainstream solution in the industry. Guo et al. [8] partition a city into
a grid-based map and use 3D convolution to capture the spatiotemporal
correlation of traffic data. Geng et al. [9] build road graphs from multiple
views based on pixel images and design a Contextual Gated Recurrent Neural
Network (CGRNN) for information integration. Li et al. [10] model traffic flow
as a diffusion process for directed graph, encoder-decoder structure and
scheduled sampling are further proposed to reduce error accumulation and
improve model performance. Zheng et al. [11] skip convolution layers and
benefit from multi-head attention layers to embed spatial and temporal
information.

The traffic road network provides a natural premium connection map for traffic
flow analysis, nevertheless, we argue that its inherent regime is not limited
to node-to-node impacts. For example, the inflow of a highway service station
is equal to its outflow within a certain period of time. Therefore, it is
interesting to analyze the relations between different traffic variables. In
recent years, causal discovery [12] opens up the door to the analysis of
causal relations behind statistical correlations and is widely explored in the
computational biology domain [13, 14]. However, to the best of our knowledge,
few works consider causal observation for spatiotemporal forecasting.

In this work, we focus on improving traffic speed prediction with insight
revealed by a causal discovery program. Speed change is thoroughly used for
causal discovery, as it is a good indication of traffic conditions and can
lead to changes in neighboring roads. In summary, our main contributions are:

  * •

We define legitimate traffic variables for causal discovery and conduct dense
experiments on the PEMS-BAY dataset to learn their causal relations. We also
find out that this knowledge can be generalized to the METR-LA dataset.

  * •

We present an effective and efficient framework to capture spatial-temporal
dependencies, which is named Graph Spatial-Temporal Network Based on Causal
Insight (GT-CausIn). The core idea is to assemble causal insight, spatial
dependency modeling, and temporal dependency modeling in a way that
information can flexibly flow between different perspectives at different
scales.

  * •

We validate our model on two real-world public datasets and achieve state-of-
the-art performance.

##  II Related work

Causal structural discovery The objective of causal structural discovery is to
combine statistical and logical reasoning to obtain causal relations between
variables, which are represented by a graph. Fast Causal Inference (FCI) was
proposed in [15], which estimates the Markov equivalence class of causal
Maximal Ancestral Graph (MAG). Recently, many papers came up with ideas to
improve the high computational cost of causal inference. Zečević et al. [16]
step towards this problem by learning interventional distributions with an
over-parameterized sum-product network. Akbari et al. [17] propose a recursive
constraint-based method that removes a specific type of variables at each
iteration. Rohekar et al. [18] further introduce Iterative Causal Discovery
(ICD) based on FCI. Each iteration of ICD refines a graph recovered from
previous iterations with a smaller set of conditions and higher statistical
reliability. As a result, ICD requires fewer conditional independence tests
and is more computational resources friendly. In this work, we benefit from
ICD to discover the causal structure of defined variables.

Spatial dependency Because of its intrinsic advantages of exploring
topological spaces, graph-based neural networks stand out from diverse deep
learning models and achieve so far the best results. Cao et al. [19] develop
the StemGNN layer to extract in-series and intra-series correlations with the
help of Graph Fourier Transformer (GFT) and Inverse Graph Fourier Transformer
(IGFT). Derrow-Pinion et al. [20] rethink the relations between road segments
and collects information from node level, segment level, and super-segment
level. Dynamic transition matrices are adopted in many works [21, 22, 23],
where nodes with similar embeddings are assumed to have higher connection
weights. However, there may be some problems with this intuitive hypothesis.
Matrices are fixed after training, and embedding similarity may yaw because of
sensor dysfunction and extreme external situations (e.g., hot/cold weather),
as a result, the adaptive matrix can get biased. In this work, we employ graph
diffusion layers designed in [10] with predefined roads as a transition
matrix, which provides reliable physical connections in the real world.

Temporal dependency The most popular way to discover temporal dependency is
through RNN-based networks. Long Short-Term Memory neural network (LSTM) is
adopted in the work [24, 25], while Li et al. [10] and Jin et al. [26] arrange
RNN block in their core model structure. Geng et al. [9] and Lv et al. [27]
apply Gated Recurrent Unit (GRU). Due to the time-consuming problem and
complex gate mechanism of RNN-based networks, we assign TCN to model temporal
dependency, which captures long sequences in a non-recursive manner [22, 23].

##  III Preparatory work

###  III-A Problem formulation

In a traffic forecasting scenario, sensors are scattered on the road to record
traffic speed at a given time, and thus can be represented as a directed
graph. In this work, a graph is denoted by
$\mathcal{G}(\mathcal{V},\mathcal{E},\mathbf{W})$ where $\mathcal{V}$ is the
node set ($|\mathcal{V}|=N$), $\mathcal{E}$ is the edge set and $\mathbf{W}$
is the adjacency matrix ($|\mathbf{W}|=N\times N$). Main notations used in the
paper are summarised in Appendix B. In the real world, the adjacency matrix is
a representation of proximity and connectivity between node sensors.

Suppose each node’s feature is of dimension $P$ (e.g., velocity, volume,
occupancy), at timestamps $t$, a graph signal can be defined as
$\mathbf{X}^{(t)}(\mathbf{X}^{t}\in\mathbb{R}^{N\times P})$. The predicted
signal at time $t$ is noted as
$\hat{\mathbf{X}}^{(t)}(\hat{\mathbf{X}}^{(t)}\in\mathbb{R}^{N\times Q})$,
where $Q$ is the number of the desired features. We consider the prediction
task as using $T^{\prime}$ prior timestamps to predict $T$ following, in other
words, given graph $\mathcal{G}$ and complementary information $\mathcal{I}$
(e.g., day of the week), we aim to find the best mapping function $f$ from
space $\mathbb{R}^{N\times P\times T^{\prime}}$ to $\mathbb{R}^{N\times
Q\times T}$:

$(\mathcal{G};\mathcal{I};\mathbf{X}^{(t-T^{\prime}+1)},\mathbf{X}^{(t-T^{\prime})},\cdots,\mathbf{X}^{(t)})\\\
\xrightarrow{f}(\hat{\mathbf{X}}^{(t+1)},\cdots,\hat{\mathbf{X}}^{(t+T-1)},\cdots,\hat{\mathbf{X}}^{(t+T)})$
(1)

In this work, we focus on using historical one-hour speed data for 15min,
30min, and 60min ahead speed forecasting $(P=Q=1)$. More specifically, when
the speed readings are aggregated into 5 minutes windows, it turns out to be:

$(\mathcal{G};\mathcal{I};\mathbf{X}^{(t-11)},\mathbf{X}^{(t-10)},\cdots,\mathbf{X}^{(t)})\\\
\xrightarrow{f}(\hat{\mathbf{X}}^{(t+3)},\hat{\mathbf{X}}^{(t+6)},\hat{\mathbf{X}}^{(t+12)})$
(2)

###  III-B Causal discovery

We need to clarify that, though causal discovery comes from statistical
ideology, it indeed surpasses the latter. One ordinary example is that,
despite the strong positive correlation between ice cream sales on the beach
and the number of drownings, taking ice cream does not lead to drowning.
Factors behind this correlation are the season and the temperature. Therefore,
causal structure discovery is necessary to reveal the faithful relation behind
observational data.

The causal structure discovery analyzes the causal relation between different
variables. In this work, we adopt ICD [18] as our causal discovery program.
The input of ICD is a set of variable values, and the output is the graph
adjacency matrix, in which each element represents the importance of the
causal relation between variables. The algorithm is only implemented on the
PEMS-BAY dataset since there are too many missing values in METR-LA (8.11%)
than in PEMS-BAY (0.003%). The dataset will be further described in Section
V-A.

Defining variables is a crucial step for causal discovery. Compared to the
speed itself, the speed variation overall represents better the changes in
traffic volume as well as the occurrence of external events, as it more
clearly expresses the changes in traffic. Figure 1 represents road sensors as
directed graph nodes. The speed detected by each sensor is a time series, and
speed changes are obtained by taking the difference between adjacent time
slices. What’s more, Appendix C shows that speed variation is a natural
Gaussian distribution, which is necessarily required for ICD.

Figure 1: Example of deriving speed variation of each node with 5 timestamps
as input.

Given node $v_{i}$, we note $\mathbb{I}_{1}(i),\mathbb{O}_{1}(i)$ the set of
its first-order in-neighbors/out-neighbors and
$\mathbb{I}_{2}(i),\mathbb{O}_{2}(i)$ its second-order in-neighbors/out-
neighbors. They can be formally formulated as below.

$\displaystyle\mathbb{I}_{1}(i)=\\{j|(j,i)\in\mathcal{E},i\neq j\\}$ (3)
$\displaystyle\mathbb{O}_{1}(i)=\\{j|(i,j)\in\mathcal{E},i\neq j\\}$ (4)
$\displaystyle\mathbb{I}_{2}(i)=\\{j|\exists
k,k\not\in\\{i,j\\},(j,k)\in\mathcal{E},(k,i)\in\mathcal{E},i\neq j\\}$ (5)
$\displaystyle\mathbb{O}_{2}(i)=\\{j|\exists
k,k\not\in\\{i,j\\},(i,k)\in\mathcal{E},(k,j)\in\mathcal{E},i\neq j\\}$ (6)

At the graph level, the input is denoted as $\mathbf{X}$, the integrated
features along the adjacency matrix of first-order in-neighbors/out-neighbors
are denoted as $\mathbf{I}_{1},\mathbf{O}_{1}$ , which can be mathematically
expressed as below.

$\displaystyle\mathbf{W}_{I1}=\mathbf{W}^{T}-\mathbf{I},\mathbf{W}_{O1}=\mathbf{W}-\mathbf{I}$
(7) $\displaystyle\begin{aligned}
\mathbf{I}_{1}=(\text{diag}(\mathbf{W}_{I1}\mathbf{1})^{-1}\mathbf{W}_{I1})\mathbf{X}\end{aligned}$
(8)
$\displaystyle\mathbf{O}_{1}=(\text{diag}(\mathbf{W}_{O1}\mathbf{1})^{-1}\mathbf{W}_{O1})\mathbf{X}$
(9)

For second order in-neighbors/out-neighbors, they are denoted as
$\mathbf{I}_{2},\mathbf{O}_{2}$:

$\displaystyle\begin{aligned}
\mathbf{W}_{I2}=\mathbf{W}_{I1}\mathbf{W}_{I1},\text{diag}(\mathbf{W}_{I2})=0\end{aligned}$
(10)
$\displaystyle\mathbf{W}_{O2}=\mathbf{W}_{O1}\mathbf{W}_{O1},\text{diag}(\mathbf{W}_{O2})=0$
(11) $\displaystyle\begin{aligned}
\mathbf{I}_{2}=(\text{diag}(\mathbf{W}_{I2}\mathbf{1})^{-1}\mathbf{W}_{I2})\mathbf{X}\end{aligned}$
(12)
$\displaystyle\mathbf{O}_{2}=(\text{diag}(\mathbf{W}_{O2}\mathbf{1})^{-1}\mathbf{W}_{O2})\mathbf{X}$
(13)

Given time $t$ and node $i$, for each timestamp $t_{0}$ between $t$ and
$t+30\min$, we define variables to be the speed variation of
$X^{(t_{0})}_{i},I_{1i}^{(t_{0})},O_{1i}^{(t_{0})},I_{2i}^{(t_{0})},O_{2i}^{(t_{0})}$.
Since the PEMS-BAY dataset aggregates traffic speed readings into 5 minutes
windows, causal variables will be stored in a set of 6-time slices, resulting
in 30 causal inferred variables below. The output of ICD is thus a $30\times
30$ matrix $\mathbf{C}$, where $C_{i,j}$ represents the causal relation
between variable $i$ and variable $j$.

$\displaystyle\mathbb{V}=$ $\displaystyle\\{V_{0},V_{1},\dots,V_{30}\\}$
$\displaystyle=$
$\displaystyle\\{X^{(t)}_{i},I_{1i}^{(t)},O_{1i}^{(t)},I_{2i}^{(t)},O_{2i}^{(t)},$
$\displaystyle\cdots,X^{(t+5)}_{i},I_{1i}^{(t+5)},O_{1i}^{(t+5)},I_{2i}^{(t+5)},O_{2i}^{(t+5)}\\}$
(14)

We randomly sample data in the dataset and then process it to obtain the above
30 causal variables. To reduce accidental error, we sample 2000 samples each
time as input of ICD, and repeat 100 times to sum the output matrix
$\mathbf{C}$.

A strong link between $v_{i},\mathbb{I}_{1}(i)$ and $\mathbb{O}_{1}(i)$ is
revealed from the causal discovery program, which is double-checked with
statistical correlation and case analysis in Appendix A.

##  IV Methodology

In this section, we start by formulating the studied problem, then present an
overview of our model structure, and end with a detailed introduction of its
crucial component layers.

###  IV-A Model structure overview

Figure 2 is a brief illustration of GT-CausIn structure. The main components
of GT-CausIn model consist of the causal insight layer, graph diffusion layer,
TCN layer, merge layer, and inherent feature layer. A detailed explanation of
each layer can be found in Section IV-B.

Figure 2: System architecture for the Graph Spatial-Temporal Network Based on
Causal Insight. Time series are fed into the causal insight layer and a series
of GT blocks. Skip connections are used at the end of each GT block to guard
useful information. Periodic features are finally merged to make predictions.

Graph diffusion layers facilitate information propagation through spatial
nodes, but not inter-time series communication. TCN layer, which focuses on
convolution along the temporal dimension, becomes complementary to graph
diffusion layers. For the sake of expression simplicity, we name GT block as
the combination of one graph diffusion layer and one TCN layer. Serializing
multiple GT blocks can thus disclose inter-spatiotemporal correlations.

The input speed, together with its differentiated variant, is fed into a
causal insight layer, where the neighbor-level information is carefully
investigated. Two fully connected layers are stacked to better fit vector
space after the first GT block.

The number of GT block $L$ is not limited to a certain value, bigger $L$ means
a bigger receptive field, but model complexity is increased at the same time
and important long-sequence information may be dispersed, the effects of $L$
will be further studied in Section V-D. Skip connections are added to stack
hidden states learned at different scales, which eases feature propagation for
short-term and long-term prediction.

###  IV-B Component layers

Causal insight layer

Figure 3: A brief illustration of the causal insight layer. Node features are
first concatenated with neighbor embedding features, then pass to an attention
layer to reveal interaction influence. Three matrices are followed to learn
individual influence on each node.

From Section III-B, we get the information that special attention on the
first-order in-neighbors/out-neighbors need to be paid. The causal insight
layer is therefore designed to embed this prior knowledge. To be cautious with
the possible mutual influence between $\mathbb{I}_{1}(i)$,
$\mathbb{O}_{1}(i)$, and different nodes $v_{i}$, we adopt a self-attention
mechanism to learn inherent relations between them. Neighbor-level and node-
level embedding spaces are followed to convert to individual node spaces.

The speed variation is concatenated with speed to form the input of the causal
insight layer. With the speed variation of each node and road connection
matrix, we can deduce the embedding of neighbor levels. Therefore, the whole
process is formulated as below and is depicted in Figure 3.

$\displaystyle\mathbf{I}_{1}=\mathbf{D}_{I1}^{-1}$
$\displaystyle(\mathbf{W}^{T}-\mathbf{I})\mathbf{X},\mathbf{O}_{1}=\mathbf{D}_{O1}^{-1}(\mathbf{W}-\mathbf{I})\mathbf{X}$
(15) $\displaystyle\mathcal{S}$
$\displaystyle=\mathbf{X}\bigoplus\mathbf{I}_{1}\bigoplus\mathbf{O}_{1}$ (16)
$\displaystyle\mathcal{S}^{{}^{\prime}}$
$\displaystyle=\mathbf{X}^{{}^{\prime}}\bigoplus\mathbf{I}_{1}^{{}^{\prime}}\bigoplus\mathbf{O}_{1}^{{}^{\prime}}$
$\displaystyle=\text{attention}(\mathcal{S}\mathbf{W}_{q},\mathcal{S}\mathbf{W}_{k},\mathcal{S}\mathbf{W}_{v})\mathbf{W}^{O}$
(17) output
$\displaystyle=\mathbf{W}_{x}\mathbf{X}^{{}^{\prime}}+\mathbf{W}_{i}\mathbf{I}_{1}^{{}^{\prime}}+\mathbf{W}_{o}\mathbf{O}_{1}^{{}^{\prime}}$
(18)

where
$\mathbf{D}_{I1}=\text{diag}((\mathbf{W}^{T}-\mathbf{I})\mathbf{1}),\mathbf{D}_{O1}=\text{diag}((\mathbf{W}-\mathbf{I})\mathbf{1})$,
$\mathbf{1}$ is the all one vector.

Graph diffusion layer Spatial dependency is modeled in this part. Consider a
directed graph $\mathcal{G}$, the diffusion process is modeled by a Markov
process with the transition matrix $\mathbf{D}_{O}^{-1}\mathbf{W}$ and a
starting probability $\alpha\in[0,1]$
($\mathbf{D}_{O}=\text{diag}(\mathbf{W}\mathbf{1})$). Teng et al. [28] show
that the ultimate distribution is converged to a stationary contribution:

$\mathcal{P}=\sum_{k=0}^{\infty}\alpha(1-\alpha)^{k}(\mathbf{D}_{O}^{-1}\mathbf{W})^{k}$
(19)

where $k$ is the diffusion step. Based on this theory, Li et al. [10] present
a formal expression of the graph signal diffusion process, which truncates the
diffusion process to be $K$ steps and assigns trainable weights for each step.

For each input feature $p\in\\{1,\cdots,P\\}$ and output feature
$q\in\\{1,\cdots,Q\\}$,

$\displaystyle\mathbf{H}_{:,q}=\mathbf{X}_{:,p\star\mathcal{G}}f_{\mathrm{\theta}_{q,p,:,:}}$
$\displaystyle=\sum_{k=0}^{K-1}(\theta_{q,p,k,1}(\mathbf{D}_{O}^{-1}\mathbf{W})^{k}+\theta_{q,p,k,2}(\mathbf{D}_{I}^{-1}\mathbf{W}^{T})^{k})\mathbf{X}_{:,p}$
(20)

where $\mathbf{X}$ is the layer input, $\mathbf{H}$ the layer output,
$\bm{\theta}\in\mathbb{R}^{Q\times P\times K\times 2}$ the filter parameter,
$\mathbf{D}_{I}=\text{diag}(\mathbf{1}\mathbf{W})$ and $\mathbf{W}^{T}$ the
transpose of graph adjacency matrix $\mathbf{W}$.

Different from the popular spectral GCN method proposed by Kipf et al. [29],
which is only applicable for undirected graphs, graph diffusion layers
simulate well the diffusion process of directed graphs, and become a natural
choice for traffic scenarios.

TCN layer We adopt the dilated causal convolution described in Yu et al. [30]
as our TCN layer. In contrast to RNN-based networks, dilated causal
convolution can capture information over long-range sequences while getting
rid of the gradient explosion problem. Concerning dilated causal convolution,
stacking dilated layers ensures the exponential receptive field, and zero
padding guarantees the temporal causal order since only historical information
is involved to predict the current time step. Considering input feature
dimension as $P$, output dimension as $Q$, given a time-series input
$\mathbf{X}\in\mathbb{R}^{N\times P\times T}$ and a filter with trainable
parameter $\mathrm{\theta}\in\mathbb{R}^{Q\times P\times K}$, the dilated
causal convolution at time $t$ can be formulated as:

$\mathbf{H}_{:,q,t}=\mathbf{X}_{:,p,:}\star
f_{\theta_{q,p,k}}(t)=\sum_{k=0}^{K-1}\theta_{q,p,k}x_{:,p,t-d*k}$ (21)

where $K$ is the kernel size, $d$ is the dilation factor, which decides the
skipping distance.

Merge layer Although GT block extracts dependencies in the temporal and
spatial dimension, information may still be lost for such a long series
forecasting task. As a result, we keep the embedding after each GT block and
concatenate them to enhance information flow over great lengths.

Inherent feature layer We notice that speed changes periodically for stations
as shown in Figure 4. Hence, we extract the day of the week, the month of the
year, and historic speed as inherent features. Day of the week and month of
the year features are represented as one-hot vectors while historic speed is
the average speed of past days at the same moment. Since people’s travel
routines are very different on weekdays and weekends, we only use the past
five weekday records or the past two weekend records to calculate historic
speed.

(a) Station 400714

(b) Station 400147

(c) Station 401403

Figure 4: Periodic speed changes of stations of PEMS-BAY within one week.

The inherent features are integrated as follows: each feature goes through a
dense layer to get embedding $\mathbf{E}_{D}$ (day of the week),
$\mathbf{E}_{M}$ (month of the year), and $\mathbf{E}_{H}$ (historic speed).
Together with embedding $\mathbf{E}_{m}$ obtained from merge layer, all
features are fused into a long feature matrix by concatenation, e.g.,
$\mathbf{F}=\mathbf{E}_{m}\bigoplus\mathbf{E}_{D}\bigoplus\mathbf{E}_{M}\bigoplus\mathbf{E}_{H}$,
then we stack two fully connected layers upon $\mathbf{F}$ to ensure that the
output is of the desired shape.

##  V Experiments

###  V-A Dataset description

We conduct dense experiments on two public datasets: PEMS-BAY [10] and METR-LA
[31]. PEMS-BAY is collected by California Transportation Agencies (CalTrans)
Performance Measurement System (PeMS). It contains 325 highway sensors of the
Bay Area from Jan 1st 2017 to Jun 30th 2017. For METR-LA, 207 loop detectors
on the highway of Los Angeles County are selected from Mar 1st 2012 to Jun
30th 2012. For both datasets, traffic speed readings are aggregated into 5
minutes windows. Based on the road connection map, we adopt the same way as Li
et al. [10] to get the graph adjacency matrix. With pairwise road network
distances between sensors, the adjacency matrix representing connection from
any sensor $i$ to any sensor $j$ is expressed as:

$W_{i,j}=\exp{(\frac{-d(v_{i},v_{j})^{2}}{\sigma^{2}})}\text{if
}d(v_{i},v_{j})\leq\kappa,\text{otherwise }0$ (22)

where $\sigma$ is the standard deviation, $\kappa$ the threshold, and
$d(v_{i},v_{j})$ the distance between station $v_{i}$ and station $v_{j}$. For
both datasets, 80% of the data serves as the training set, 10% as the test
set, and the remaining 10% as the validation set.

###  V-B Experiment setting

All models are implemented with PyTorch [32] and trained with an Adam
optimizer with an annealing learning rate. Experiments on dataset METR-LA are
carried out on two Nvidia Tesla V100 servers and PEMS-BAY on two Nvidia Tesla
A100 servers. The detailed parameter setting for all experiments executed is
available in Appendix F.

All methods are evaluated based on three commonly used metrics in traffic
forecasting, which are Mean Absolute Error (MAE), Mean Absolute Percentage
Error (MAPE), and Root Mean Squared Error (RMSE). The formula of all metrics
can be found in Appendix E. MAE is used as the loss function. Missing values
are linearly filled for training, and are excluded for results evaluation.

###  V-C Comparison with baselines

TABLE I: Performance comparison with baseline models for traffic speed
forecasting.

Data | Method | 15min | 30min | 60min
---|---|---|---|---
MAE | RMSE | MAPE | MAE | RMSE | MAPE | MAE | RMSE | MAPE
METR-LA | DCRNN | 2.77 | 5.38 | 7.30% | 3.15 | 6.45 | 8.8% | 3.60 | 7.59 | 10.5%
Graph Wavenet | 2.69 | 5.15 | 6.90% | 3.07 | 6.22 | 8.37% | 3.53 | 7.37 | 10.01%
GMAN | 4.04 | 8.53 | 10.2% | 4.59 | 9.85 | 11.69% | 5.33 | 11.21 | 13.60%
ST-GRAT | 2.60 | 5.07 | 6.61% | 3.01 | 6.21 | 8.15% | 3.49 | 7.42 | 10.01%
SLCNN | 2.53 | 5.18 | 6.70% | 2.88 | 6.15 | 8.00% | 3.30 | 7.20 | 9.70%
DGCRN | 2.62 | 5.01 | 6.63% | 2.99 | 6.05 | 8.02% | 3.44 | 7.19 | 9.73%
DMSTGCN | 2.85 | 5.54 | 7.54% | 3.26 | 6.56 | 9.19% | 3.72 | 7.55 | 10.96%
PGCN | 2.70 | 5.16 | 6.98% | 3.08 | 6.22 | 8.38% | 3.54 | 7.36 | 9.94%
GT-NoCausIn | 2.82 | 5.52 | 7.18% | 3.27 | 6.71 | 8.90% | 3.95 | 8.24 | 11.67%
GT-CausIn | 2.61 | 5.07 | 6.60% | 2.73 | 5.34 | 7.09% | 3.06 | 6.11 | 8.30%
Improvement | -3.2% | -1.2% | 0.2% | 5.2% | 11.7% | 11.4% | 7.3% | 15.0% | 14.4%
PEMS-BAY | DCRNN | 1.38 | 2.95 | 2.90% | 1.74 | 3.97 | 3.90% | 2.07 | 4.74 | 4.90%
Graph Wavenet | 1.30 | 2.74 | 2.73% | 1.63 | 3.70 | 3.67% | 1.95 | 4.52 | 4.63%
GMAN | 1.34 | 2.82 | 2.81% | 1.62 | 3.72 | 3.63% | 1.86 | 4.32 | 4.31%
ST-GRAT | 1.29 | 2.71 | 2.67% | 1.61 | 3.69 | 3.63% | 1.95 | 4.54 | 4.64%
SLCNN | 1.44 | 2.90 | 3.00% | 1.72 | 3.81 | 3.90% | 2.03 | 4.53 | 4.80%
DGCRN | 1.28 | 2.69 | 2.66% | 1.59 | 3.63 | 3.55% | 1.89 | 4.42 | 4.43%
DMSTGCN | 1.33 | 2.83 | 2.80% | 1.67 | 3.79 | 3.81% | 1.99 | 4.54 | 4.78%
PGCN | 1.30 | 2.73 | 2.72% | 1.62 | 3.67 | 3.63% | 1.92 | 4.45 | 4.55%
GT-NoCausIn | 1.43 | 2.85 | 3.00% | 1.79 | 3.84 | 4.01% | 2.25 | 4.90 | 5.44%
GT-CausIn | 1.30 | 2.54 | 2.68% | 1.47 | 2.92 | 3.10% | 1.70 | 3.37 | 3.68%
Improvement | -1.6% | 5.6% | -0.8% | 7.5% | 19.6% | 12.7% | 8.6% | 22.0% | 14.6%

We compare GT-CausIn with various popular traffic forecasting models in recent
years, including (1) DCRNN: Diffusion Convolutional Recurrent Neural Network
[10]; (2) Graph Wavenet [33], which develops an adaptive dependency matrix
based on node embedding; (3)GMAN: Graph Multi-Attention Network for Traffic
Prediction [11], of which the main structure is an encoder and a decoder both
with several attention blocks; (4) ST-GRAT: Spatio-temporal Graph Attention
Networks for Accurately Forecasting Dynamically Changing Road Speed [34],
which includes spatial attention, temporal attention, and spatial sentinel
vectors; (5) SLCNN: Structure Learning Convolution Neural Network [35], which
learns the graph structure information with convolutional methods; (6) DGCRN:
Dynamic Graph Convolutional Recurrent Network (DGCRN) [36], which filters node
embedding to generate dynamic graph at each time step; (7) DMSTGCN: Dynamic
and Multi-faceted Spatio-temporal Deep Learning for Traffic Speed Forecasting
[37], which provides a multi-faceted fusion module to incorporate the hidden
states learned at different stages; (8) PGCN: Progressive Graph Convolutional
Networks for Spatial-temporal Traffic Forecasting [23], which constructs
progressive adjacency matrices by learning in training and test phases.

Our approach is different from all the approaches above, we first use a causal
discovery program to discover a general rule between node neighbors of
different orders, then build a model named GT-CausIn to enhance prediction
performance. Although the causality relation is discovered with dataset PEMS-
BAY, it generalizes as well on the other dataset METR-LA. Besides, the
proposed causal insight layer is integrated with directed graph diffusion
layers to model the spatial dependency. TCN layers and skip connections are
further used to discover the temporal dependency.

A comparison of GT-CausIn and other baseline models is shown in Table I. The
causal insight layer is taken out from GT-CausIn and the rest model is named
as GT-NoCausIn. All the baseline results are taken from its original paper if
it is available and from [23] if not111Only GMAN and DMSTGCN authors implement
experiments on different datasets..

From Table I, we can observe a large margin between GT-NoCausIn and GT-CausIn
, which shows that causal insights can significantly improve prediction
accuracy. We also notice that models with dynamic training graph generally
outperforms other models, however, GT-CausIn achieves the best overall result,
especially for mid-term (30min) and long-term (60min) prediction. We need to
remark that DCRNN [10] shares similar graph diffusion layers and GMAN [11]
widely adopt the attention mechanism, yet our model exceeds them on all
metrics, which justifies our model structure.

###  V-D Ablation study

Causal effects In GT-CausIn, we put special attention to relations between
node $v_{i}$, $\mathbb{I}_{1}(i)$ and $\mathbb{O}_{1}(i)$, nonetheless, this
conclusion drawn from the causal discovery program may be questioned.
Therefore, we replace Equation (16, 17, 18) as below to check its
effectiveness, naming this model GT-BadCausIn.

$\displaystyle\mathcal{S}$
$\displaystyle=\mathbf{X}\bigoplus\mathbf{X}\bigoplus\mathbf{X}$ (23)
$\displaystyle\mathcal{S}^{{}^{\prime}}$
$\displaystyle=\mathbf{X}^{{}^{\prime}}_{0}\bigoplus\mathbf{X}_{1}^{{}^{\prime}}\bigoplus\mathbf{X}_{2}^{{}^{\prime}}$
$\displaystyle=\text{attention}(\mathcal{S}\mathbf{W}_{q},\mathcal{S}\mathbf{W}_{k},\mathcal{S}\mathbf{W}_{v})\mathbf{W}^{O}$
(24) output
$\displaystyle=\mathbf{W}_{x_{0}}\mathbf{X}_{0}^{{}^{\prime}}+\mathbf{W}_{x_{1}}\mathbf{X}_{1}^{{}^{\prime}}+\mathbf{W}_{x_{2}}\mathbf{X}_{2}^{{}^{\prime}}$
(25)

Therefore, GT-BadCausIn integrates no neighbor embedding in causal insight
layer.

TABLE II: Ablation study on causal effects.

Data | Method | 15min | 30min | 60min
---|---|---|---|---
MAE | RMSE | MAPE | MAE | RMSE | MAPE | MAE | RMSE | MAPE
METR-LA | GT-BadCausIn | 2.65 | 5.17 | 6.72% | 2.86 | 5.71 | 7.58% | 3.18 | 6.50 | 8.85%
GT-CausIn | 2.61 | 5.07 | 6.60% | 2.73 | 5.34 | 7.09% | 3.06 | 6.11 | 8.30%
Improvement | 1.5% | 1.9% | 1.8% | 4.5% | 6.5% | 6.5% | 3.8% | 6.0% | 6.2%
PEMS-BAY | GT-BadCausIn | 1.35 | 2.63 | 2.80% | 1.56 | 3.13 | 3.34% | 1.81 | 3.64 | 4.01%
GT-CausIn | 1.30 | 2.54 | 2.68% | 1.47 | 2.92 | 3.10% | 1.70 | 3.37 | 3.68%
Improvement | 3.7% | 3.4% | 4.3% | 5.8% | 6.7% | 7.2% | 5.0% | 7.4% | 8.2%

Experiments results are summarised in Table II. The impact brought by neighbor
embedding in the causal insight layer is positive for all metrics on both
datasets, especially for mid-term (30min) and long-term (60min) predictions.
The reason behind may be the propagation time of causal effects. Besides, we
can observe that the improvement is smaller on METR-LA than on PEMS-BAY, which
may result from: (1) a larger percentage of missing values in METR-LA (8.11%)
than in PEMS-BAY (0.003%); (2) causal discovery program is only executed on
the PEMS-BAY dataset.

Number of stacking GT blocks In Section IV-A, we introduced that the number of
GT blocks $L$ is not limited to a certain value. Larger $L$ roughly
corresponds to a larger spatiotemporal receptive field and enables the model
to capture broader spatiotemporal dependency. We compare model performance
with different $L\in\\{3,4,5,6\\}$ on MAE of different time horizons in Figure
5. We can observe a steady rise when $L$ goes up. Since the gap between $L=3$
and $L=4$ is more important than between $L=4$ and $L=5$, taking the increased
model complexity into account, we adopt $L=4$ as our best parameter.

(a) MAE/15min

(b) MAE/30min

(c) MAE/60min

Figure 5: Parameter analysis of $L$ (number of stacking GT layers) on METR-LA
dataset

###  V-E Model interpretation

To better understand our model, we compare it with GT-BadCausIn on 60min ahead
prediction task. We shift the input window to get continuous predictions in
Figure 6, observing a notable performance margin between GT-CausIn and GT-
BadCausIn. As described in Section V-D, the only structural difference between
these two models is the dissimilar perspectives integrated. For GT-BadCausIn,
only node-level information $\mathbf{X}$ is investigated. For GT-CausIn,
neighbor-level information $\mathbf{I}$ and $\mathbf{O}$ are also included.
This structural change results in different token attention scores of the
causal insight layer and is responsible for this disparity, this is further
investigated in Appendix D.

(a) Station 71447

(b) Station 773062

(c) Station 716337

Figure 6: Figures (a), (b), and (c) offer a comparison between GT-CausIn and
GT-BadCausin for prediction task on METR-LA, differing on the target station
and target time. Each figure contains the ground truth of the considered
station and its neighbors, as well as its 60min ahead prediction. It is
notable that the predictions of GT-BadCausIn are heavily dependent on some
neighbors, and are therefore swayed off the ground truth. The same cannot be
said for GT-causIn, instead of being influenced by some stations, it has a
broader outlook of its neighbors, resulting in a much more faithful and smooth
prediction.

##  VI Conclusion

In this paper, we implemented knowledge learned from a causal discovery
program with deep learning models and proposed Graph Spatial-Temporal Network
Based on Causal Insight to capture spatiotemporal dependencies. Specifically,
we serialize graph diffusion layers and TCN layers to capture dependencies
between spatial flow and temporal flow, we also use skip connections to
guarantee information propagation for long sequences. We creatively design a
causal insight layer that focuses on stations, their first-order in-neighbors,
and their first-order out-neighbors. When evaluating our model on two real-
world traffic datasets, we achieved significantly better overall prediction
than baselines. We also conducted ablation studies to show the effectiveness
of causal knowledge and the influence of the stacking layer number. In the
future, we plan to (1) apply the proposed model to other spatial-temporal
forecasting tasks; (2) investigate scalable methods to apply to large-scale
datasets.

## References

  * [1] A. M. Nagy and V. Simon, “Survey on traffic prediction in smart cities,” _Pervasive and Mobile Computing_ , vol. 50, pp. 148–163, 2018. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1574119217306521
  * [2] E. Cascetta, _Transportation systems engineering: theory and methods_. Springer Science & Business Media, 2013, vol. 49.
  * [3] M. M. Hamed, H. R. Al-Masaeid, and Z. M. B. Said, “Short-term prediction of traffic volume in urban arterials,” _Journal of Transportation Engineering_ , vol. 121, no. 3, pp. 249–254, 1995.
  * [4] B. M. Williams and L. A. Hoel, “Modeling and forecasting vehicular traffic flow as a seasonal arima process: Theoretical basis and empirical results,” _Journal of transportation engineering_ , vol. 129, no. 6, pp. 664–672, 2003\.
  * [5] A. J. Smola and B. Schölkopf, “A tutorial on support vector regression,” _Statistics and computing_ , vol. 14, no. 3, pp. 199–222, 2004.
  * [6] T. Evgeniou, M. Pontil, and T. Poggio, “Regularization networks and support vector machines,” _Advances in computational mathematics_ , vol. 13, no. 1, pp. 1–50, 2000.
  * [7] X. Ma, Z. Dai, Z. He, J. Ma, Y. Wang, and Y. Wang, “Learning traffic as images: a deep convolutional neural network for large-scale transportation network speed prediction,” _Sensors_ , vol. 17, no. 4, p. 818, 2017.
  * [8] S. Guo, Y. Lin, S. Li, Z. Chen, and H. Wan, “Deep spatial–temporal 3d convolutional neural networks for traffic data forecasting,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 20, no. 10, pp. 3913–3926, 2019.
  * [9] X. Geng, Y. Li, L. Wang, L. Zhang, Q. Yang, J. Ye, and Y. Liu, “Spatiotemporal multi-graph convolution network for ride-hailing demand forecasting,” in _Proceedings of the AAAI conference on artificial intelligence_ , vol. 33, no. 01, 2019, pp. 3656–3663.
  * [10] Y. Li, R. Yu, C. Shahabi, and Y. Liu, “Diffusion convolutional recurrent neural network: Data-driven traffic forecasting,” in _6th International Conference on Learning Representations, ICLR 2018, Vancouver, BC, Canada, April 30 - May 3, 2018, Conference Track Proceedings_. OpenReview.net, 2018. [Online]. Available: https://openreview.net/forum?id=SJiHXGWAZ
  * [11] C. Zheng, X. Fan, C. Wang, and J. Qi, “Gman: A graph multi-attention network for traffic prediction,” in _Proceedings of the AAAI conference on artificial intelligence_ , vol. 34, no. 01, 2020, pp. 1234–1241.
  * [12] J. Pearl, _Causality_. Cambridge university press, 2009.
  * [13] S. Pranay and N. Nagaraj, “Causal discovery using compression-complexity measures,” _Journal of Biomedical Informatics_ , vol. 117, p. 103724, 2021\.
  * [14] C. Cundy, A. Grover, and S. Ermon, “Bcd nets: Scalable variational approaches for bayesian causal discovery,” _Advances in Neural Information Processing Systems_ , vol. 34, pp. 7095–7110, 2021.
  * [15] P. Spirtes, C. Glymour, R. Scheines, S. Kauffman, V. Aimale, and F. Wimberly, “Constructing bayesian network models of gene expression networks from microarray data,” 2000.
  * [16] M. Zečević, D. Dhami, A. Karanam, S. Natarajan, and K. Kersting, “Interventional sum-product networks: Causal inference with tractable probabilistic models,” _Advances in Neural Information Processing Systems_ , vol. 34, pp. 15 019–15 031, 2021.
  * [17] S. Akbari, E. Mokhtarian, A. Ghassami, and N. Kiyavash, “Recursive causal structure learning in the presence of latent variables and selection bias,” _Advances in Neural Information Processing Systems_ , vol. 34, pp. 10 119–10 130, 2021.
  * [18] R. Y. Rohekar, S. Nisimov, Y. Gurwicz, and G. Novik, “Iterative causal discovery in the possible presence of latent confounders and selection bias,” _Advances in Neural Information Processing Systems_ , vol. 34, pp. 2454–2465, 2021.
  * [19] D. Cao, Y. Wang, J. Duan, C. Zhang, X. Zhu, C. Huang, Y. Tong, B. Xu, J. Bai, J. Tong _et al._ , “Spectral temporal graph neural network for multivariate time-series forecasting,” _Advances in neural information processing systems_ , vol. 33, pp. 17 766–17 778, 2020.
  * [20] A. Derrow-Pinion, J. She, D. Wong, O. Lange, T. Hester, L. Perez, M. Nunkesser, S. Lee, X. Guo, B. Wiltshire _et al._ , “Eta prediction with graph neural networks in google maps,” in _Proceedings of the 30th ACM International Conference on Information & Knowledge Management_, 2021, pp. 3767–3776.
  * [21] B. Du, X. Hu, L. Sun, J. Liu, Y. Qiao, and W. Lv, “Traffic demand prediction based on dynamic transition convolutional neural network,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 22, no. 2, pp. 1237–1247, 2020.
  * [22] L. Han, B. Du, L. Sun, Y. Fu, Y. Lv, and H. Xiong, “Dynamic and multi-faceted spatio-temporal deep learning for traffic speed forecasting,” in _Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining_, 2021, pp. 547–555.
  * [23] Y. Shin and Y. Yoon, “Pgcn: Progressive graph convolutional networks for spatial-temporal traffic forecasting,” _arXiv preprint arXiv:2202.08982_ , 2022.
  * [24] L. Lin, Z. He, and S. Peeta, “Predicting station-level hourly demand in a large-scale bike-sharing network: A graph convolutional neural network approach,” _Transportation Research Part C: Emerging Technologies_ , vol. 97, pp. 258–276, 2018.
  * [25] H. Yao, F. Wu, J. Ke, X. Tang, Y. Jia, S. Lu, P. Gong, J. Ye, and Z. Li, “Deep multi-view spatial-temporal network for taxi demand prediction,” in _Proceedings of the AAAI conference on artificial intelligence_ , vol. 32, no. 1, 2018.
  * [26] J. Jin, D. Rong, T. Zhang, Q. Ji, H. Guo, Y. Lv, X. Ma, and F.-Y. Wang, “A gan-based short-term link traffic prediction approach for urban road networks under a parallel learning framework,” _IEEE Transactions on Intelligent Transportation Systems_ , 2022.
  * [27] M. Lv, Z. Hong, L. Chen, T. Chen, T. Zhu, and S. Ji, “Temporal multi-graph convolutional network for traffic flow prediction,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 22, no. 6, pp. 3337–3348, 2020.
  * [28] S.-H. Teng _et al._ , “Scalable algorithms for data and network analysis,” _Foundations and Trends® in Theoretical Computer Science_ , vol. 12, no. 1–2, pp. 1–274, 2016.
  * [29] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” in _5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings_. OpenReview.net, 2017\. [Online]. Available: https://openreview.net/forum?id=SJU4ayYgl
  * [30] F. Yu and V. Koltun, “Multi-scale context aggregation by dilated convolutions,” in _4th International Conference on Learning Representations, ICLR 2016, San Juan, Puerto Rico, May 2-4, 2016, Conference Track Proceedings_ , Y. Bengio and Y. LeCun, Eds., 2016. [Online]. Available: http://arxiv.org/abs/1511.07122
  * [31] H. V. Jagadish, J. Gehrke, A. Labrinidis, Y. Papakonstantinou, J. M. Patel, R. Ramakrishnan, and C. Shahabi, “Big data and its technical challenges,” _Commun. ACM_ , vol. 57, no. 7, p. 86–94, jul 2014. [Online]. Available: https://doi.org/10.1145/2611567
  * [32] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Köpf, E. Z. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, “Pytorch: An imperative style, high-performance deep learning library,” in _Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada_ , H. M. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché-Buc, E. B. Fox, and R. Garnett, Eds., 2019, pp. 8024–8035. [Online]. Available: https://proceedings.neurips.cc/paper/2019/hash/bdbca288fee7f92f2bfa9f7012727740-Abstract.html
  * [33] Z. Wu, S. Pan, G. Long, J. Jiang, and C. Zhang, “Graph wavenet for deep spatial-temporal graph modeling,” in _Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, IJCAI-19_. International Joint Conferences on Artificial Intelligence Organization, 7 2019, pp. 1907–1913. [Online]. Available: https://doi.org/10.24963/ijcai.2019/264
  * [34] C. Park, C. Lee, H. Bahng, Y. Tae, S. Jin, K. Kim, S. Ko, and J. Choo, “St-grat: A novel spatio-temporal graph attention networks for accurately forecasting dynamically changing road speed,” ser. CIKM ’20. New York, NY, USA: Association for Computing Machinery, 2020, p. 1215–1224. [Online]. Available: https://doi.org/10.1145/3340531.3411940
  * [35] Q. Zhang, J. Chang, G. Meng, S. Xiang, and C. Pan, “Spatio-temporal graph structure learning for traffic forecasting,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 34, no. 01, 2020, pp. 1177–1185.
  * [36] F. Li, J. Feng, H. Yan, G. Jin, F. Yang, F. Sun, D. Jin, and Y. Li, “Dynamic graph convolutional recurrent network for traffic prediction: Benchmark and solution,” _ACM Transactions on Knowledge Discovery from Data (TKDD)_ , 2021\.
  * [37] L. Han, B. Du, L. Sun, Y. Fu, Y. Lv, and H. Xiong, “Dynamic and multi-faceted spatio-temporal deep learning for traffic speed forecasting,” in _Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining_, ser. KDD ’21. New York, NY, USA: Association for Computing Machinery, 2021, p. 547–555. [Online]. Available: https://doi.org/10.1145/3447548.3467275
  * [38] D. Freedman and R. Pisani, “R. purves: Statistics,” 2007.
  * [39] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin, “Attention is all you need,” in _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , ser. NIPS’17. Red Hook, NY, USA: Curran Associates Inc., 2017, p. 6000–6010.

##  Appendix A Causal result validation

In this section, we carefully validate causal results obtained in Appendix
III-B in two ways: (1) coherence with the correlation between variables; (2)
case analysis.

Correlation between variables We apply the sampling algorithm mentioned in
Appendix III-B to the random data and the accident data, and adopt Pearson
correlation coefficient [38] to calculate the correlation between the
different variables. The result can be plotted in Figure 7, in which the size
of the circle indicates the strength of its relations.

Figure 7: ICD result vs data correlation, on the left: accident data, on the
right: random data. Labels represent the variables defined in Equation (14).
Red circles indicate results deduced from ICD. Gray circles indicate the
Pearson correlation coefficients between variables mined from incident data
(left) and random data (right).

Four phenomena can be observed in Figure 7: (1) The relation obtained by
causal discovery holds for both random data and accident data, indicating that
the speed changes at one station do have an effect on neighboring stations and
that this law is universal, without the need to include an additional judgment
about whether an accident has occurred; (2) The causality is mainly reflected
in spatial and temporal space, there is no clear spatiotemporal intertwined
causality; (3) There is another interesting effect that the station $v_{i}$,
its out-neighbor $\mathbb{O}_{1}(i)$ and in-neighbor $\mathbb{I}_{1}(i)$ are
causally related to each other; (4) For relations in the temporal dimension,
variables are more connected to themselves from $t$ to $t+1$ and $t+2$.

Case analysis In the real world there may be many reasons behind speed
changes, such as sudden traffic jams, weather phenomena, and car accidents. In
light of these unpredictable external factors, and to ensure that the results
and conclusions on causality seen before are founded in reality, it is
essential to plot the speed of a node and its neighbors over time and analyze
its behaviors.

To do so, a dataset of notable events is created. These events are
characterized by a sudden and lasting spike in speed variation, positive or
negative. Due to the large scale of the PEMS-BAY dataset, it is necessary to
sample random mini-batches and calculate the weighted speed change in the
station and its neighbors. Sorting through these average velocity
fluctuations, the most influential ones are selected, i.e. highest and lowest
peaks.

(a) Event at station 400714

(b) Event at station 407185

(c) Event at station 400209

Figure 8: Notable events displaying the correlation between the station and
its neighbors. Curves $I_{1}(i)$ and $O_{1}(i)$ denote the weighted sum
formulated in Equation (9) and shadow gives the speed variation interval of
the first-order in-neighbors/out-neighbors.

According to these notable events seen in Figure 8, we can confirm the mutual
impacts of speed changes between a station, its in-neighbors, and its out-
neighbors. In most cases, these three entities share a similar pattern for
speed variation, as shown in Figures 8a and 8b. We also observe a delayed
impact in some cases, as shown in Figure 8c, where the speed variation is
spread from its out-neighbors.

##  Appendix B Notation

TABLE III: Core concept definitions

Concept | Description
---|---
$\mathcal{G}(\mathcal{V},\mathcal{E})$ | Graph with node set $\mathcal{V}$ and edge set $\mathcal{E}$, $|\mathcal{V}|=N$
$v_{i}$ | The $i$-th node
$\mathbf{W}$ | Adjacency matrix of graph, $\mathbf{W}\in\mathbb{R}^{N\times N}$
$\mathbf{D}_{I},\mathbf{D}_{O}$ | In-degree/out-degree matrix
$\mathbf{1}$ | All one vector
$\mathbf{X},\hat{\mathbf{X}}$ | Graph signal and the predicted graph signal
$\mathbf{I}$ | Identity matrix
$\mathbb{I}_{1}(i),\mathbb{O}_{1}(i)$ | Set of first-order in/out-neighbors of $i$-th node
$\mathbb{I}_{2}(i),\mathbb{O}_{2}(i)$ | Set of second-order in/out-neighbors of $i$-th node
$\mathbf{I}_{1},\mathbf{O}_{1}$ | Integrated embedding of first-order in/out-neighbors
$\mathbf{I}_{2},\mathbf{O}_{2}$ | Integrated embedding of second-order in/out-neighbors
$\text{diag}(\mathbf{A})$ | Column vector containing diagonal elements of $\mathbf{A}$ ($\mathbf{A}$ is a matrix)
$\text{diag}(A)$ | Diagonal matrix whose entries are elements of $A$ ($A$ is a vector)

##  Appendix C Data visualization

Below are the data distributions for PEMS-BAY dataset.

(a) Speed variation distribution

(b) Speed distribution

Figure 9: PEMS-BAY data distribution.

##  Appendix D Attention score analysis

In this section we focus on the 60min ahead prediction task on station 717477
of METR-LA at 11:15:00. As seen in Figure 6a, there is a performance gain for
GT-CausIn in comparison to GT-BadCausIn model. This is due to a structural
difference seen in V-D which includes neighbor-level information in the causal
insight layer.

To look into what is happening inside the causal insight layer, we take out
the score matrices. We revise that the attention scores are calculated by
processing the score matrix on a softmax layer [39], which guarantees that the
sum of each row equals one. Given input $\mathbf{X}$, the attention layer is
formulated as below.

$\displaystyle\mathbf{Q}=\mathbf{X}\mathbf{W}_{q},\mathbf{K}=\mathbf{X}\mathbf{W}_{k},\mathbf{V}=\mathbf{X}\mathbf{W}_{v}$
(26)
$\displaystyle\text{attention}=\text{softmax}\frac{(\mathbf{Q}\mathbf{K}^{T})}{\sqrt{d_{k}}}\mathbf{V}$
(27)

where $\mathbf{W}_{q},\mathbf{W}_{k}$ and $\mathbf{W}_{v}$ are trainable
weights, $\frac{1}{\sqrt{d_{k}}}$ is the scalar factor, $\mathbf{Q}$ is query
matrix, $\mathbf{K}$ is the key matrix and $\mathbf{V}$ is the value matrix.
The attention score is
$\mathcal{S}=\frac{(\mathbf{Q}\mathbf{K}^{T})}{\sqrt{d_{k}}}$, and its element
$S_{i,j}$ weights the attention $\mathbf{X}_{i}$ pays to $\mathbf{X}_{j}$.

Figure 10: Heat map of attention score in causal insight layer for GT-
BadCausIn Figure 11: Heat map of attention score in causal insight layer for
GT-CausIn

In our case, the attention input is
$\mathcal{S}=\mathbf{X}\bigoplus\mathbf{I}_{1}\bigoplus\mathbf{O}_{1}$ for GT-
CausIn and $\mathcal{S}=\mathbf{X}\bigoplus\mathbf{X}\bigoplus\mathbf{X}$ for
GT-BadCausIn, as indicated in Equation (17) for GT-CausIn and Equation (24)
for GT-BadCausIn. We take rows corresponding to the station $v_{i}$ to focus
on how it is impacted. Afterward, we subdivide the columns into three parts,
each part contains tokens of all graph nodes from different perspectives, as
noted by the labels in the bottom.

With these tokens in hand, the attention scores are visualized with heat maps.
Figure 10 shows attention scores inside the GT-BadCausIn model, that includes
solely $\mathbf{X}$ tokens. In contrast, in Figure 11, the attention scores
inside GT-CausIn are not constrained to $\mathbf{X}$, but also including
$\mathbf{I}_{1}$ and $\mathbf{O}_{1}$ tokens. Finally the performance
difference seen in Figure 6a can be explained by different attention scores in
the Figures 10 and 11. In the first case, the attention is mainly distributed
on some nodes. In the second case, the attention has a global view on all the
nodes thanks to neighbor-level tokens, e.g., $\mathbf{I}_{1},\mathbf{O}_{1}$.

##  Appendix E Metrics

Note $\mathrm{x}=\\{x_{1},\cdots,x_{n}\\}$ the ground truth,
$\hat{\mathrm{x}}=\\{\hat{x}_{1},\cdots,\hat{x}_{n}\\}$ predicted value, and
$\Omega$ indices of observed samples, the metrics are defined as follows.

$\displaystyle\text{MAE}(\mathrm{x},\hat{\mathrm{x}})$
$\displaystyle=\frac{1}{|\bf\Omega|}\sum_{i\in{\bf\Omega}}|x_{i}-\hat{x}_{i}|$
(28) $\displaystyle\text{RMSE}(\mathrm{x},\hat{\mathrm{x}})$
$\displaystyle=\sqrt{\frac{1}{|\bf\Omega|}\sum_{i\in{\bf\Omega}}|(x_{i}-\hat{x}_{i})^{2}}$
(29) $\displaystyle\text{MAPE}(\mathrm{x},\hat{\mathrm{x}})$
$\displaystyle=\frac{1}{|\bf\Omega|}\sum_{i\in{\bf\Omega}}|\frac{x_{i}-\hat{x}_{i}}{x_{i}}|$
(30)

##  Appendix F Detailed experimental settings

For all models, we stack $L=4$ GT blocks (except ablation study on $L$) and
the output dimension of each block is 8. In each graph diffusion layer, the
maximum number of steps of random walks is 3. The kernel size is 3 for all TCN
layers. We use a step learning rate with $\gamma=0.5$ for all experiments and
adopt a grid search to find the best initial learning rate for each
experiment, as shown in Table IV.

TABLE IV: Learning rate setting Data | Model | Initial lr\start step\step size
---|---|---
PEMS-BAY | GT-NoCausIn | 0.002\50\20
GT-CausIn | 0.006\50\20
GT-BadCausIn | 0.007\50\20
METR-LA | GT-NoCausIn | 0.008\180\50
GT-CausIn | 0.004\180\50
GT-BadCausIn | 0.004\180\50
GT-CausIn$(L=3)$ | 0.005\180\50
GT-CausIn$(L=5)$ | 0.006\180\50
GT-CausIn$(L=6)$ | 0.004\180\50


# The “magic” angle in the self-assembly of colloids suspended in a nematic
host phase

Sergej Schlotthauer  Stranski-Laboratorium für Physikalische und Theoretische
Chemie, Fakultät für Mathematik und Naturwissenschaften, Technische
Universität Berlin, Straße des 17. Juni 115, 10623 Berlin, GERMANY  Tillmann
Stieger  Stranski-Laboratorium für Physikalische und Theoretische Chemie,
Fakultät für Mathematik und Naturwissenschaften, Technische Universität
Berlin, Straße des 17. Juni 115, 10623 Berlin, GERMANY  Michael Melle
Stranski-Laboratorium für Physikalische und Theoretische Chemie, Fakultät für
Mathematik und Naturwissenschaften, Technische Universität Berlin, Straße des
17. Juni 115, 10623 Berlin, GERMANY  Marco G. Mazza  Max-Planck-Institut für
Dynamik und Selbstorganisation, Am Faßberg 17, 37077 Göttingen, GERMANY
Martin Schoen  Stranski-Laboratorium für Physikalische und Theoretische
Chemie, Fakultät für Mathematik und Naturwissenschaften, Technische
Universität Berlin, Straße des 17. Juni 115, 10623 Berlin, GERMANY  Department
of Chemical and Biomolecular Engineering, 911 Partners Way, North Carolina
State University, Raleigh, NC 27695, U.S.A.

###### Abstract

Using extensive Monte Carlo (MC) simulations of colloids immersed in a nematic
liquid crystal we compute an effective interaction potential via the local
nematic director field and its associated order parameter. The effective
potential consists of a local Landau-de Gennes (LdG) and a Frank elastic
contribution. Molecular expressions for the LdG expansion coefficients are
obtained via classical density functional theory (DFT). The DFT result for the
LdG parameter $A$ is improved by locating the phase transition through finite-
size scaling theory. We consider effective interactions between a pair of
homogeneous colloids with Boojum defect topology. In particular, colloids
attract each other if the angle between their center-of-mass distance vector
and the far-field nematic director is about $30^{\circ}$ which settles a long-
standing discrepancy between theory and experiment. Using the effective
potential in two-dimensional MC simulations we show that self-assembled
structures formed by the colloids are in excellent agreement with experimental
data.

###### pacs:

61.30.-v,61.30.Jf,82.70Dd,05.10.Ln

If a liquid crystal is in the nematic phase the overall orientation of its
molecules (i.e., mesogens) can be described quantitatively by the non-local
unit vector (i.e., the far-field nematic director) $\bm{\widehat{n}}_{0}$ [1].
Immersing a colloidal particle in this nematic host phase gives rise to a
director field $\bm{\widehat{n}}\left(\bm{r}\right)$ such that sufficiently
close to the colloid’s surface, $\bm{\widehat{n}}\left(\bm{r}\right)$ and
$\bm{\widehat{n}}_{0}$ may differ. The deviation between
$\bm{\widehat{n}}\left(\bm{r}\right)$ and $\bm{\widehat{n}}_{0}$ is caused by
the specific anchoring of mesogens at the surface of the colloid. Depending on
details of the host phase $\bm{\widehat{n}}\left(\bm{r}\right)$ can be of such
dazzling complexity that experts are just beginning to unravel its structural
details [2].

The mismatch between $\bm{\widehat{n}}\left(\bm{r}\right)$ and
$\bm{\widehat{n}}_{0}$ also gives rise to effective interactions between
several colloids that are mediated by the nematic host [3]. These interactions
may therefore be used to self-assemble the colloids into supramolecular
entities in a controlled (i.e., directed) manner. This way ordered assemblies
of colloids of an enormously complex structure with rich symmetries may be
built that would not exist without the ordered structure of the host phase [4,
5].

The complex self-assembled structures formed by the colloids are also of
practical importance. For instance, taking as a specific example dielectric
colloids it could be demonstrated that the propagation of light through a
self-assembled ordered colloidal arrangement is affected in a way similar to
the propagation of electrons in a semiconductor crystal [6]. Hence, ordered
periodic assemblies of colloids are already discussed within the framework of
novel photonic devices with fascinating properties [7].

Clearly, to use the effective interaction potential for the self-assembly of
colloids in a nematic host phase the molecular origin of the potential itself
must be understood. Our motivation to contribute to such an improved
understanding goes back to an observation made some time ago by Poulin and
Weitz [8]. They found experimentally that in a nematic phase the colloidal
center-to-center distance vector $\bm{r}_{12}$ forms a “magic” angle of
$\theta\approx 30^{\circ}$ with $\bm{\widehat{n}}_{0}$ if the mesogens at the
surfaces of the colloidal pair are anchored in a locally planar fashion.
Hence, near an isolated colloid a so-called Boojum defect would arise under
these conditions [9].

This experimental observation has resisted a quantitative theoretical
explanation to date. In previous theoretical attempts a much larger angle of
about $50^{\circ}$ is usually found [8, 10]. This number is based upon
calculations where one employs the electrostatic analog of the Boojum defect
topology [8]. In fact, as stated explicitly by Poulin and Weitz “This
theoretical value is different from the experimentally observed value for
$\theta$ $\ldots$ since the theory is a long-range description that does not
account for short-range effects” [8]. Another motivation for our work is the
more recent experimental observation that between a pair of colloids in a
nematic host repulsive and attractive forces act depending on $\theta$ [10].
For example, at $\theta\approx 30^{\circ}$ the colloids attract each other
whereas at $\theta=0^{\circ}$ and $90^{\circ}$ repulsion between the colloids
is observed.

To unravel the persisting discrepancy between theory and experiment we employ
a combination of Monte Carlo (MC) simulations in the isothermal-isobaric
ensemble, two-dimensional (2D) MC simulations in the canonical ensemble,
classical density functional theory (DFT), concepts of finite-size scaling
(FSS), and Landau-de Gennes (LdG) theory to investigate the effective
interaction between a pair of spherical, chemically homogeneous colloids
mediated by a nematic host phase.

To model the host phase we adopt the so-called Hess-Su model. In this model
mesogen-mesogen interactions are described by an isotropic core where
$\varepsilon$ and $\sigma$ set energy and length scale, respectively.
Superimposed to the isotropic core are anisotropic attractions of respective
strengths $\varepsilon_{1}\varepsilon$ and $\varepsilon_{2}\varepsilon$ where
the dimensionless anisotropy parameters
$2\varepsilon_{1}=-\varepsilon_{2}=0.08$ throughout this work. Under these
conditions the Hess-Su model exhibits isotropic-nematic (IN) phase transitions
[11].

The colloid-mesogen interaction is modeled via short-range repulsive
interactions and an attractive Yukawa tail where we take its inverse Debye
screening length $\lambda\sigma=0.50$ [12]. Mesogens at the colloids’ surfaces
are anchored in a locally planar fashion. This setup is then placed between
structureless, planar solid substrates separated by a distance
$s_{\mathrm{z}}=24\sigma$. Mesogens at the substrates are anchored such that
their longer axes point along the $x$-axis $\bm{\widehat{e}}_{\mathrm{x}}$.
Under these conditions a Boojum defect topology emerges at a single, isolated
colloid. Colloids are immersed in the host phase such that their center-to-
center distance vector is given by $\bm{r}_{12}=\left(x_{12},y_{12},0\right)$.

We employ dimensionless units, that is length is given in units of $\sigma$,
energy in units of $\varepsilon$, and temperature in units of
$\varepsilon/k_{\mathrm{B}}$ ($k_{\mathrm{B}}$ Boltzmann’s constant). Other
derived units are then expressed as combinations of these basic ones as usual
[12]. In particular, we set temperature $T=0.90$ and pressure $P=1.80$ such
that the host phase is nematic at a mean number density $\rho\approx 0.90$.
The hard-core radius of each colloid is $r_{0}=3.00$. Other conditions of the
MC simulations are exactly the same as in Ref. 12 where additonal details of
the model can also be found.

Figure 1: (Color online) (a)–(c) Defect topologies for a colloidal pair with
locally planar surface anchoring of mesogens separated by $\bm{r}_{12}$;
$\cos\theta=\bm{r}_{12}\cdot\bm{\widehat{n}}_{0}/r_{12}$,
$\bm{\widehat{n}}_{0}\cdot\bm{\widehat{e}}_{x}=1$, and
$r_{12}=\left|\bm{r}_{12}\right|$. (d)–(f) As (a)–(c) but projected onto the
$x$–$y$ plane. Attached color bars give $S\left(x,y\right)$ and dashes
indicate $\bm{\widehat{n}}\left(\bm{r}\right)$.

Results of our MC simulations shown in Fig. 1(a) indicate that parts of the
Boojum defects interact forming a torus. As $\theta$ increases, the torus is
“ripped apart” [Fig. 1(b)]. Eventually, a handle-like defect topology emerges
at $\theta\simeq 90^{\circ}$ [Fig. 1(c)].

Defect regions around the colloids are visualized by shading them if the local
nematic order parameter $S\left(\bm{r}\right)\leq 0.20$. We obtain
$S\left(\bm{r}\right)$ numerically as the largest eigenvalue of the local
alignment tensor [12]. The eigenvector $\bm{\widehat{n}}\left(\bm{r}\right)$
associated with $S\left(\bm{r}\right)$ is the director field.

The latter is illustrated in Figs. 1(d)–1(f). The plots indicate relatively
localized regions of low $S\left(\bm{r}\right)$ in the vicinity of the
colloids and that $\bm{\widehat{n}}\left(\bm{r}\right)$ is bent in ways that
depend on the specific defect topology (i.e., on $\theta$).

Naturally, the reduction of $S\left(\bm{r}\right)$ and the bending of
$\bm{\widehat{n}}\left(\bm{r}\right)$ causes the free-energy density
$f\left(\bm{r}\right)$ of the system to increase locally relative to that of
the pure host phase without the colloids. Consequently, we adopt

$\displaystyle\Delta f\left(\bm{r}\right)\\!$ $\displaystyle=$
$\displaystyle\\!A\left(T,\rho\right)S^{2}\left(\bm{r}\right)\\!+\\!B\left(T,\rho\right)S^{3}\left(\bm{r}\right)\\!+\\!C\left(T,\rho\right)S^{4}\left(\bm{r}\right)$
(1)
$\displaystyle+\frac{K}{2}\left\\{\left[\bm{\nabla}\cdot\bm{\widehat{n}}\left(\bm{r}\right)\right]^{2}+\left[\bm{\nabla}\times\bm{\widehat{n}}\left(\bm{r}\right)\right]^{2}\right\\}-f_{0}$

where the first three terms on the right side correspond to a local LdG free-
energy density $f_{\mathrm{LdG}}\left(\bm{r}\right)$,
$f_{0}=AS^{2}+BS^{3}+CS^{4}$ is the LdG free-energy density obtained under the
same thermodynamic conditions but in the absence of the colloids, and $\Delta
f_{\mathrm{LdG}}\left(\bm{r}\right)\equiv
f_{\mathrm{LdG}}\left(\bm{r}\right)-f_{0}$. Coefficients $A$, $B$, and $C$ are
coefficients in the LdG expansion and $S$ is the global nematic order
parameter.

The two terms on the second line of Eq. (1) correspond to the local Frank
free-energy density $f_{\mathrm{el}}\left(\bm{r}\right)$ that accounts for
elastic distortions of the director field where $K$ is an elastic constant. We
consider here the so-called one-constant approximation in which it is assumed
that splay, twist, and bend deformations of
$\bm{\widehat{n}}\left(\bm{r}\right)$ contribute equally to
$f_{\mathrm{el}}\left(\bm{r}\right)$. It has recently been shown [13] that the
one-constant approximation is an excellent approximation for the present model
system because of the small aspect ratio of the mesogens. Under the present
thermodynamic conditions, $K=1.66$. We then obtain $f_{\mathrm{el}}$ by
numerically differentiating $\bm{\widehat{n}}\left(\bm{r}\right)$ [12].

We assume that both the local LdG contribution in Eq. (1) and $f_{0}$ are
governed by the same set $A$, $B$, and $C$. Moreover, $\Delta
f\left(\bm{r}\right)$ in Eq. (1) does not account for either fluctuations in
$S\left(\bm{r}\right)$ or $\bm{\widehat{n}}\left(\bm{r}\right)$ and therefore
constitutes a mean-field expression. Notice also that using in Eq. (1)
$S\left(\bm{r}\right)$ and $\bm{\widehat{n}}\left(\bm{r}\right)$ from MC is
advantageous because then both quantities correspond to an equilibrium
situation. Conventionally, $S\left(\bm{r}\right)$ and
$\bm{\widehat{n}}\left(\bm{r}\right)$ are treated as variational functions in
the ansatz in Eq. (1) which bears the risk that the numerical minimization of
the functional $\Delta
f\left[S\left(\bm{r}\right),\bm{\widehat{n}}\left(\bm{r}\right)\right]$ may
miss the true equilibrium solution.

To use Eq. (1), $A$, $B$, and $C$ are required. Whereas these quantities are
notoriously difficult to compute for reasons described by Eppenga and Frenkel
a long time ago [14], Gupta and Ilg have devised a new approach that works
reliably for mesogens with a relatively large aspect ratio [15]. In practice,
however, we observed that the method of Gupta and Ilg does not work well for
our model fluid where mesogens have a rather small aspect ratio of only
$1.26$.

Because Eq. (1) constitutes a mean-field expression we resort to mean-field
DFT alternatively where [16]

$\beta\Delta
f_{\mathrm{or}}=\rho\int\limits_{-1}^{1}\mathrm{d}x\,\overline{\alpha}\left(x\right)\ln\left[\overline{\alpha}\left(x\right)\right]+\rho^{2}\sum\limits_{\begin{subarray}{c}l=2\\\
l\text{ even}\end{subarray}}^{\infty}S_{l}^{2}u_{l}$ (2)

is the difference in free-energy density of the nematic relative to the
isotropic phase. In Eq. (2), $x=\cos\vartheta$ where $\vartheta$ is the
azimuthal angle, $\beta=1/k_{\mathrm{B}}T$, $\overline{\alpha}\left(x\right)$
is the orientation distribution function, and members of the set
$\left\\{u_{l}\right\\}$ account for the contribution of anisotropic mesogen-
mesogen interactions to the free-energy density. Because of the uniaxial
symmetry of the nematic phase we expand

$\overline{\alpha}\left(x\right)=\frac{1}{2}+\sum\limits_{\begin{subarray}{c}l=2\\\
l\text{
even}\end{subarray}}^{\infty}\frac{2l+1}{2}S_{l}P_{l}\left(x\right)\equiv\frac{1}{2}+\xi\left(x\right)$
(3)

in terms of Legendre polynomials $\left\\{P_{l}\left(x\right)\right\\}$. We
assume $\bm{\widehat{n}}_{0}\cdot\bm{\widehat{e}}_{\mathrm{z}}=1$ and $0\leq
S_{l}\leq 1$ are order parameters.

We then insert the expression on the far right side of Eq. (3) into Eq. (2)
and expand the integrand in terms of $\xi$ around $\xi=0$ (i.e., at the IN
phase transition). Retaining in this expansion only the leading term of $\xi$
for $l=2$ and neglecting terms proportional to $S_{2}^{n}$ ($n\geq 5$) allows
us to rewrite Eq. (2) as

$\Delta
f_{\mathrm{or}}=a\left(\rho\right)\left(T-T^{\ast}\right)S_{2}^{2}-\frac{8\rho
k_{\mathrm{B}}T}{105}S_{2}^{3}+\frac{4\rho k_{\mathrm{B}}T}{35}S_{2}^{4}$ (4)

where $a\left(\rho\right)=2\rho k_{\mathrm{B}}/5$ and $T^{\ast}=-5\rho
u_{2}/2k_{\mathrm{B}}$ is the temperature at which the nematic phase becomes
thermodynamically stable. Assuming that $S=S_{2}$ we equate terms of equal
power in $S$ in $f_{0}$ and Eq. (4) which yields molecular expressions for the
LdG constants $A$, $B$, and $C$. In particular, $A$ changes sign at
$T=T^{\ast}$, $B<0$, and $C>0$ as they must at a first-order phase transition
[1].

One also notices that the value of $T^{\ast}$ depends on $u_{2}$ where the
precise form of $u_{2}$ is a consequence of the level of sophistication at
which pair correlations are treated within mean-field DFT [11]. For example,
at simple mean-field (SMF) level,
$u_{2}=-32\pi\varepsilon_{1}\varepsilon\sigma^{3}/15$ is a constant. At the
more elaborate modified mean-field (MMF) level, $u_{2}$ becomes a function of
$T$ [see Eqs. (3.7) and (3.8) of Ref. 11]. It turns out that at SMF level,
$T^{\ast}$ is underestimated whereas at MMF level it is overestimated.

To overcome this problem we determine $T^{\ast}$ via FSS. Following Ref. 17 we
first calculate the coexistence temperature $T_{\mathrm{IN}}\simeq 1.02$ at
the IN phase transition. It is given as the intersection of the second-order
Binder cumulants of $S$ for different system sizes [17]. From the expression
$T^{\ast}=T_{\mathrm{IN}}-2B^{2}/9aC$ [18] and using $B$, $C$, and $a$ from
DFT, $T^{\ast}$ can easily be determined. Notice also that
$S_{\mathrm{IN}}=-2B/3C=\frac{4}{9}$ irrespective of $T_{\mathrm{IN}}$ [1]
whereas MMF DFT predicts this value of $S_{\mathrm{IN}}$ only to be a
threshold reached for sufficiently high $T_{\mathrm{IN}}$ (Fig. 2 of Ref. 11)
thus pointing to a certain deficiency of LdG theory.

Figure 2: (Color online)
$\mathcal{F}_{\mathrm{el}}=\int\mathrm{d}\bm{r}\,f_{\mathrm{el}}\left(\bm{r}\right)$
($\bullet$) (left ordinate) and
$\Delta\mathcal{F}_{\mathrm{LdG}}=\int\mathrm{d}\bm{r}\,\Delta
f_{\mathrm{LdG}}\left(\bm{r}\right)$ ($\blacksquare$) (right ordinate) as
functions of $\theta$ for $r_{12}=2r_{0}$ (see Fig. 1).

Plots of $\mathcal{F}_{\mathrm{el}}$ and $\Delta\mathcal{F}_{\mathrm{LdG}}$ in
Fig. 2 illustrate the impact of a colloidal pair on the free energy of the
host phase. Both quantities vary nonmonotonically with the angle $\theta$ and
exhibit minima at $\theta\simeq 30^{\circ}$ in agreement with the experimental
findings of Poulin and Weitz [8]. Because deformations of
$\bm{\widehat{n}}\left(\bm{r}\right)$ cost free energy,
$\mathcal{F}_{\mathrm{el}}>0$. Similarly, the presence of the colloids reduces
$S\left(\bm{r}\right)$ such that in some regions $S>S\left(\bm{r}\right)$ (see
Fig. 1). Because the host phase without the colloids is deep in the nematic
phase, $f_{0}<0$ such that $\Delta\mathcal{F}_{\mathrm{LdG}}>0$ as well.

That both $\mathcal{F}_{\mathrm{el}}$ and $\Delta\mathcal{F}_{\mathrm{LdG}}$
become minimal at about the same $\theta$ indicates that destortions of
$\bm{\widehat{n}}\left(\bm{r}\right)$ and a local reduction of nematic order
are coupled. However, deformations of $\bm{\widehat{n}}\left(\bm{r}\right)$
turn out to be more important than reduction of nematic order because
$\mathcal{F}_{\mathrm{el}}$ exceeds $\Delta\mathcal{F}_{\mathrm{LdG}}$ by
between one and two orders of magnitude over the entire range of $\theta$’s.
This conclusion is drawn on the basis of plots in Fig. 2 and by noticing that
for both curves the ground state is the same, namely
$\bm{\widehat{n}}\left(\bm{r}\right)=\bm{\widehat{n}}_{0}$
($\mathcal{F}_{\mathrm{el}}=0$) and $S\left(\bm{r}\right)=S$
($\Delta\mathcal{F}_{\mathrm{LdG}}=0$).

Figure 3: (Color online)
$\Delta\mathcal{F}_{\mathrm{eff}}/\Delta\mathcal{F}_{\mathrm{B}}$ as a
function of relative positions of the colloids in the $x$–$y$ plane (see
attached color bar). The white semicircle at the center represents a reference
colloid.

Results presented in Fig. 2 have been obtained for two colloids in contact
with each other. However, the general physical picture reflected by Fig. 2 is
preserved if besides $\theta$, $r_{12}$ is varied, too. To that end we realize
from Eq. (1) that $\lim_{r_{12}\to\infty}\Delta f\left(\bm{r}\right)=2\Delta
f_{\mathrm{B}}\left(\bm{r}\right)$ where $\Delta
f_{\mathrm{B}}\left(\bm{r}\right)$ is the local free energy density of two
isolated Boojum defects relative to the same ground state used above. Taking
$\Delta\mathcal{F}_{\mathrm{B}}=\int\mathrm{d}\bm{r}\,\Delta
f_{\mathrm{B}}\left(\bm{r}\right)$ allows us to introduce
$\Delta\mathcal{F}_{\mathrm{eff}}\equiv\mathcal{F}_{\mathrm{el}}+\Delta\mathcal{F}_{\mathrm{LdG}}-2\Delta\mathcal{F}_{\mathrm{B}}$
as the effective potential acting between a pair of colloids and mediated by
the nematic host.

A map of $\Delta\mathcal{F}_{\mathrm{eff}}$ in Fig. 3 shows that for
$\theta=0^{\circ}$, $\Delta\mathcal{F}_{\mathrm{eff}}$ is strongly repulsive
in a relatively localized region. This is a consequence of the merger of parts
of the Boojum defect illustrated by Figs. 1(a) and 1(d). In agreement with
plots in Fig. 2 we see that $\Delta\mathcal{F}_{\mathrm{eff}}$ is attractive
if $r_{12}$ is sufficiently small where the absolute minimum of
$\Delta\mathcal{F}_{\mathrm{eff}}$ is found at $\theta\approx 30^{\circ}$.

One also notices from Fig. 3 a small repulsive barrier in
$\Delta\mathcal{F}_{\mathrm{eff}}$ as $\theta$ approaches $90^{\circ}$ and
$7\lesssim r_{12}\lesssim 10$. Hence, a pair of colloids at $\theta\approx
0^{\circ}$ and at sufficiently large $r_{12}$ and
$75^{\circ}\lesssim\theta\lesssim 90^{\circ}$ would repel each other whereas
those forming an angle of $\theta\approx 30^{\circ}$ would attract each other.
These findings are in excellent agreement with experimental observations [Fig.
2(b) of Ref. 10].

Figure 4: (Color online) 2D MC configurations
($\bm{\widehat{n}}_{0}\cdot\bm{\widehat{e}}_{\mathrm{x}}=1$). (a)
$\phi=N_{\mathrm{coll}}\pi r_{0}^{2}/s_{\mathrm{x}}s_{\mathrm{y}}=0.065$, (b)
$\phi=0.234$ ($s_{\mathrm{x}}=s_{\mathrm{y}}=50$).

Taking $\Delta\mathcal{F}_{\mathrm{eff}}$ as an effective, pairwise additive
potential we perform standard Metropolis 2D MC simulations of
$N_{\mathrm{coll}}$ colloids modeling the nematic host phase implicitly.
Technically, $\Delta\mathcal{F}_{\mathrm{eff}}$ is stored at nodes of a
regularly spaced grid in the $x$–$y$ plane; the actual value of
$\Delta\mathcal{F}_{\mathrm{eff}}$ at $\bm{r}_{12}$ is obtained by bilinear
interpolation between the four nearest nodes. The simulations are carried out
in the canonical ensemble. Results in Fig. 4(a) show that at low packing
fraction $\phi$ the colloids tend to form linear chains of an angle of about
$30^{\circ}$ with $\bm{\widehat{n}}_{0}$. At higher $\phi$ the snapshot in
Fig. 4(b) reveals more extended two-dimensional structures. Plots in both
parts of Fig. 4 are in excellent qualitative agreement with experimental
findings (see Fig. 1 of Ref. 10).

To summarize we used a combination of MC simulations, FSS, and mean-field DFT
to compute the effective interaction potential between a pair of colloids
immersed in a nematic liquid crystal. The colloids are chemically homogeneous
and anchor mesogens in a locally planar fashion at their surface. On accound
of the mismatch between this local alignment and $\bm{\widehat{n}}_{0}$ a
Boojum defect topology emerges at an isolated colloid. If two such colloids
approach each other the Boojum defects interact such that the precise topology
changes with the angle $\theta$ formed between the distance vector connecting
the centers of the colloidal pair and $\bm{\widehat{n}}_{0}$.

As a result of the topological change repulsive and attractive effective
interactions arise. These are dominated by the distortion of
$\bm{\widehat{n}}\left(\bm{r}\right)$ whereas the accompanying reduction of
local nematic order is negligible. Most notably, the distribution of regions
in which the effective interaction potential
$\Delta\mathcal{F}_{\mathrm{eff}}$ is attractive or repulsive matches
experimental results reported by Smalyukh et al. despite their much larger
colloids [10].

It is particularly gratifying that the most favorable angle we find is
$\theta\approx 30^{\circ}$ in agreement with the work by Poulin and Weitz [8]
and Smalyukh et al. [10]. Our work therefore offers the first quantitative
theoretical explanation of earlier experimental observations. Moreover, we
show that it is the relatively short-range effects that are responsible for
the observed attraction and repulsion between nematic colloids thereby
confirming the earlier conjecture by Poulin and Weitz [8].

###### Acknowledgements.

We acknowledge financial support from Deutsche Forschungsgemeinschaft through
the International Graduate Research Training Group 1524. S. S. and M. M. are
grateful for discussions with Prof. C. K. Hall (NCSU).

## References

  * [1] P. G. de Gennes and J. Prost, The physics of liquid crystals, (Oxford Science Publications, Oxford, 1995).
  * [2] V. S. R. Jampani et al., Phys. Rev. E 84, 031703 (2011).
  * [3] K. Izaki and Y. Kimura, Phys. Rev. E 87, 062507 (2013).
  * [4] U. Ognysta et al., Langmuir 25, 12092 (2009).
  * [5] H. Qi and T. Hegmann, J. Mater. Chem. 16, 4197 (2006).
  * [6] M. Humar et al., Nat. Photonics 3, 595 (2009).
  * [7] I. Muševic et al., J. Phys.: Condens. Matter 23, 284112 (2011).
  * [8] P. Poulin and D. A. Weitz, Phys. Rev. E 57, 626 (1998).
  * [9] P. Poulin et al., Science 275, 1770 (1997).
  * [10] I. I. Smalyukh et al., Phys. Rev. Lett. 95, 157801 (2005).
  * [11] S. Giura and M. Schoen, Phys. Rev. E 90, 022507 (2014).
  * [12] M. Melle et al., J. Chem. Phys. 136, 194703 (2012).
  * [13] T. Stieger et al., J. Chem. Phys. 140, 054905 (2014).
  * [14] R. Eppenga and D. Frenkel, Mol. Phys. 52, 1303 (1984).
  * [15] B. Gupta and P. Ilg, Polymers, 5, 328 (2013).
  * [16] S. Giura et al., Phys. Rev. E 87, 012313 (2013).
  * [17] M. Greschek and M. Schoen, Phys. Rev. E 83, 011704 (2011).
  * [18] L. Senbetu and C.-W. Woo, Mol. Cryst. Liq. Cryst. 84, 101 (1982).