# A Computational Note on the Graphical Ridge in High-dimension

A. Bekkera,b, A. Kheyria and M. Arashic CONTACT A. Bekker. Email:
<EMAIL_ADDRESS>aDepartment of Statistics, Faculty of Natural and
Agricultural Sciences, University of Pretoria, Pretoria, South Africa
bCentre of Excellence in Mathematical and Statistical Sciences, Johannesburg,
South Africa cDepartment of Statistics, Faculty of Mathematical Sciences,
Ferdowsi University of Mashhad, Mashhad, Iran

###### Abstract

This article explores the estimation of precision matrices in high-dimensional
Gaussian graphical models. We address the challenge of improving the accuracy
of maximum likelihood-based precision estimation through penalization.
Specifically, we consider an elastic net penalty, which incorporates both
$L_{1}$ and Frobenius norm penalties while accounting for the target matrix
during estimation. To enhance precision matrix estimation, we propose a novel
two-step estimator that combines the strengths of ridge and graphical lasso
estimators. Through this approach, we aim to improve overall estimation
performance. Our empirical analysis demonstrates the superior efficiency of
our proposed method compared to alternative approaches. We validate the
effectiveness of our proposal through numerical experiments and application on
three real datasets. These examples illustrate the practical applicability and
usefulness of our proposed estimator.

###### keywords:

Graphical model; Graphical lasso estimation; High-dimensional Gaussian
graphical models; Penalization; Precision matrix estimation; Ridge estimation

##  1 Introduction

In recent years, there has been a substantial amount of research dedicated to
investigating the estimation of the covariance matrix and its inverse,
commonly known as the precision matrix. The motivation behind these studies
arises from the critique that commonly employed likelihood-based approaches
yield inaccurate estimates for the covariance and precision matrices,
particularly in scenarios where the number of variables exceeds the available
data points. One commonly employed approach is incorporating an $L_{1}$
penalty into the log-likelihood function. Nearly concurrently, [20], [3], [8]
all investigated the estimation of the precision matrix, represented by a
graphical lasso. This method has attracted a lot of interest and turned into
an active topic of study due to the sparse solution ([7]; [5]; [18]; [12];
[6]; [2]).
While sparsity in precision matrix estimation is often desirable, there are
situations where precise representations of the high-dimensional precision
matrix are fundamentally advantageous. Moreover, it is important to note that
the true precision matrix does not necessarily need to exhibit sparsity,
meaning it may not have a significant number of zero elements. In such
scenarios, a regularization approach is employed to effectively shrink the
estimated elements of the precision matrix. Several studies have explored the
use of ridge regularization in the log-likelihood function by incorporating a
Frobenius norm penalty. Relevant works in this area include [11]; [13]; [17];
[16]; [15], among others.
A challenge in the regularized precision estimation is to find a closed-form
solution dealing with the $L_{1}$ penalty. Indeed element-wise estimator is
not the case, unlike vector-variate estimation. This paper considers a
combination of alternative ridge and $L_{1}$ penalty in a 2-step approach. The
main goal of this paper is to present and compare the efficiency gained by
considering our proposed methodology for estimating the precision matrix in
the Gaussian graphical model compared to some existing methods. Secondly, we
propose a closed-form estimator for the precision matrix when we have an
elastic net penalty. Thus, we organize the rest of this article as follows. In
Section 2, we review two types of the archetypal ridge and the alternative
ridge estimators. Then, we propose our 2-step approach with elastic net type
penalty to find a non-sparse solution for the precision matrix in 3. In the
numerical study, Section 4 confirms that the proposed estimation is comparable
with some alternative estimations, and in most cases, it behaves better.
Section 5 provides examples of real-world data applications strengthening our
methodology and argument. To provide a quick summary, we discuss our key
findings briefly in Section 6.

##  2 Existing Methods

Gaussian graphical models are widely popular in graphical modeling due to
their ease of handling. One key advantage of Gaussian graphical models is that
both their conditional and unconditional dependence structures are fully
determined by the precision matrix of the underlying distribution. Estimating
the covariance matrix in such models can be achieved using standard
methodologies, such as maximum likelihood (ML) estimation or its regularized
variants.
Mathematically, a Gaussian graphical model is represented by a pair
$(\mathbf{z},\mathcal{G})$, where $\mathbf{z}$ is a random vector following a
multivariate normal distribution $\mathcal{N}_{p}(\bm{\mu},\bm{\Sigma})$.
Here, $\bm{\mu}\in\mathbb{R}^{p}$ represents an arbitrary mean vector, and
$\bm{\Sigma}$ is a symmetric and positive definite covariance matrix denoted
as $\bm{\Sigma}\succ 0$. The graph $\mathcal{G}$ corresponds to the structural
representation of the conditional independence relationships among the
variables.
In Gaussian graphical models, the graph $\mathcal{G}$ and the precision matrix
$\bm{\Theta}:=\bm{\Sigma}^{-1}$ are closely interconnected. The zero entries
in the precision matrix, denoted by $\theta_{ij}=0$, correspond to pairs of
variables $\mathbf{z}_{i}$ and $\mathbf{z}_{j}$ that are conditionally
independent given all the other variables, i.e.,
$\mathbf{z}_{i}\mathbf{z}_{j}|\mathbf{z}_{\\{1,\ldots,p\\}\setminus i,j}$. In
other words, there is no edge between these two variables (nodes) in the graph
$\mathcal{G}$. Additionally, two variables $\mathbf{z}_{i}$ and
$\mathbf{z}_{j}$ are independent if and only if there is no path of edges
connecting them in the graph $\mathcal{G}$.

The maximum likelihood (ML) estimation is an analytical technique used to
estimate the precision matrix $\bm{\Theta}$. The ML estimator for
$\bm{\Theta}$ is typically expressed as

$\hat{\bm{\Theta}}=\arg\underset{\bm{\Theta}\succ
0}{\max}\\{\log\det(\bm{\Theta})-\mathop{\rm
tr}\nolimits(\bm{S}\bm{\Theta})\\},$ (1)

Where $\bm{S}=\mathbf{z}^{T}\mathbf{z}/n$ represents the sample covariance
matrix, the ML estimator of $\bm{\Theta}$ is given by $\bm{S}^{-1}$. However,
when applying this technique to estimate $\bm{\Theta}$, two issues may arise.
Firstly, in the high-dimensional case where the number of variables $p$
exceeds the number of observations $n$, the empirical covariance matrix
$\bm{S}$ becomes singular, preventing the inversion required to estimate
$\bm{\Theta}$. Even if $p$ and $n$ are approximately equal, and $\bm{S}$ is
non-singular, the maximum likelihood estimate for $\bm{\Theta}$ tends to have
a relatively high variance. Secondly, in a graphical model, it is often
desirable to identify pairs of variables that are disconnected and
conditionally independent, which corresponds to zero entries in the precision
matrix $\bm{\Theta}$. However, the maximum likelihood estimation given by
equation (1) generally does not yield any zero elements in the estimated
$\bm{\Theta}$.

###  2.1 Ridge Estimators

To address the limitations mentioned above, the likelihood function is
augmented with a penalty function in high-dimensional settings, leading to the
regularized maximum likelihood estimator. In the case of non-sparse solutions,
a ridge estimator based on archetypes was proposed in [11]. The corresponding
penalized log-likelihood can be expressed as follows:

$\log\det(\bm{\Theta})-(1-\lambda)\mathop{\rm
tr}\nolimits(\bm{S}\bm{\Theta})-\lambda\mathop{\rm
tr}\nolimits(\bm{\Theta}\bm{\Gamma}),$ (2)

where $\lambda\in\text{(0,1]}$ is the tuning parameter, and $\bm{\Gamma}$ is a
positive definite target matrix. Refer to [1] for more information about the
ridge notion. A common target selection is a diagonal matrix, where its
elements are the same as the corresponding elements of the sample covariance
matrix $\bm{S}$ in (2).
Solving the optimization problem (2), we have the first archetypal ridge
estimator given by

$\hat{\bm{\Theta}}^{1}(\lambda)=[(1-\lambda)\bm{S}+\lambda\bm{\Gamma}]^{-1}.$
(3)

In the second archetypal ridge approach, Schafer and Strimmer [13] introduced
a specific form of penalized log-likelihood. They considered the following
expression:

$\log\det(\bm{\Theta})-\mathop{\rm
tr}\nolimits(\bm{S}\bm{\Theta})-\lambda\mathop{\rm tr}\nolimits(\bm{\Theta}),$
(4)

where $\lambda\in(0,\infty)$ and the corresponding estimator has form

$\hat{\bm{\Theta}}^{2}(\lambda)=[\bm{S}+\lambda I_{p}]^{-1}.$

Contrary to our expectations, the penalties in equations (2) and (4) do not
correspond to the Frobenius norm penalty, which is defined as the square root
of the sum of squares of the elements of the precision matrix.
The alternative ridge approach introduced by [16] and the ridge-type operator
for precision matrix estimation (ROPE) proposed by [10] both utilize the
Frobenius norm penalty. These methods consider two different types of ridge
precision estimators, where the penalized log-likelihood function is given by

$\log\det(\bm{\Theta})-\mathop{\rm
tr}\nolimits(\bm{S}\bm{\Theta})-\frac{1}{2}\lambda||\bm{\Theta}-\bm{T}||_{F}^{2}\\},$
(5)

with $\bm{T}$ denoting a symmetric positive definite target matrix and the
penalty parameter $\lambda\in(0,\infty)$. The type-I ridge precision estimator
(or ROPE), the consequence of the optimization problem (5), is given by

$\hat{\bm{\Theta}}^{\@slowromancap
i@a}(\lambda)=\left[\left(\lambda\bm{I}_{p}+\frac{1}{4}(\bm{S}-\lambda\bm{T})^{2}\right)^{1/2}+\frac{1}{2}(\bm{S}-\lambda\bm{T})\right]^{-1}.$
(6)

The RHS and LHS limits of $\hat{\bm{\Theta}}^{\@slowromancap i@a}$ are the
inverses of the ML estimator $\bm{S}$ (if $p<n$) and the target matrix,
respectively. Target matrix $\bm{T}$ is considered positive definite to avoid
having wrong estimates in the limit.

######  Remark 1.

The type-I ridge estimator can also be represented without inversion as

$\hat{\bm{\Theta}}^{Ia}(\lambda)=\frac{1}{\lambda}\left[\left(\lambda\bm{I}_{p}+\frac{1}{4}(\bm{S}-\lambda\bm{T})^{2}\right)^{1/2}-\frac{1}{2}(\bm{S}-\lambda\bm{T})\right].$
(7)

Also, the type-II ridge estimator is obtained by considering $\bm{T}=\bm{0}$
in (6) (or (7)) as follows

$\hat{\bm{\Theta}}^{\@slowromancap
ii@a}(\lambda)=\left[\left(\lambda\bm{I}_{p}+\frac{1}{4}\bm{S}^{2}\right)^{1/2}+\frac{1}{2}\bm{S}\right]^{-1}.$
(8)

The RHS and LHS limits of $\hat{\bm{\Theta}}^{\@slowromancap ii@a}$, like the
second archetypal estimator, are the inverse of the ML estimator $\bm{S}$
(where $p<n$) and the null matrix, respectively. It should be mentioned that
[19] studied the estimator (8) in a different context.
In a recent study by [15], a generalized ridge precision estimator was
proposed, allowing for entry-wise penalization. This approach enables the
estimator to shrink towards a user-specified, nonrandom target matrix while
maintaining positive definiteness and consistency. Furthermore, the authors
derived a generalization of the graphical lasso estimator and its elastic net
counterpart. It is worth noting that Kovács et al., [9] considered elastic net
type penalization for precision matrix estimation, specifically in the
presence of a diagonal target matrix. They demonstrated that adapting the
iterative procedure proposed by [15] for the elastic net problem is feasible
but may not be computationally efficient.
All of the above set the platform for the motivation of this paper. Here, in
this paper, we consider the elastic net type penalty and propose a novel non-
sparse estimate technique for precision matrix using a 2-step strategy.

##  3 Constructing the Proposed Estimator

In this section, we begin by providing the motivation for our methodology.
Considering the continuity property of ridge estimators and disregarding the
influence of sparsity, there is a possibility that an estimator of
$\bm{\Theta}$ may deviate significantly from the true value $\Theta$; for
instance, the estimators mentioned in Section 2. Hence, it is of major
importance to find a closer estimator. Therefore, we can penalize the deviance
between any primary or initial target estimator and $\bm{\Theta}$ and tune it
to improve the final estimation using the Frobenius norm penalty. Then we add
the $L_{1}$ penalty, returning to the well-known elastic net regularization
framework, but with refinement.

Given a sample $\bm{X}_{n\times p}$ of $n$ realizations from a $p$-variate
Gaussian distribution with zero mean and positive definite covariance matrix
$\bm{\Sigma}$. Consider the following optimization problem aimed at estimating
the unknown precision matrix based on $n$ observations:

$\hat{\bm{\Theta}}(\alpha,\lambda,\bm{T})=\arg\underset{\bm{\Theta}\succ
0}{\max}\\{\log\det(\bm{\Theta})-\mathop{\rm
tr}\nolimits(\bm{S}\bm{\Theta})-\lambda(\alpha||\bm{\Theta}||_{1}+\frac{(1-\alpha)}{2}||\bm{\Theta}-\bm{T}||_{F}^{2})\\},$
(9)

where $\bm{T}$ is a known positive semi-definite diagonal target matrix. Also
$\lambda\geq 0$ and $\alpha\in[0,1]$ are tuning parameters. The target
$\bm{T}$ can be any of the estimators mentioned in Section 2 or a sparse
version like the graphical lasso.

We can rewrite the optimization problem (9) as

$\underset{\bm{\Theta}\succ 0}{max}\underset{\parallel
U\parallel_{\infty}\leq\lambda\alpha}{min}\log\det(\bm{\Theta})-\mathop{\rm
tr}\nolimits(\bm{\Theta},\bm{S}+\bm{U})-\frac{\lambda(1-\alpha)}{2}||\bm{\Theta}-\bm{T}||_{2}^{2},$
(10)

where we express the $L_{1}$ norm as

$||\bm{\Theta}||_{1}=\underset{\parallel\bm{U}\parallel_{\infty}\leq
1}{max}tr(\bm{\Theta}\bm{U}),$ (11)

and $||\bm{U}||_{\infty}$ represents the maximum absolute value element of the
symmetric matrix $\bm{U}$. To achieve a closed-form solution for (10), we have
the following results. All the proofs are provided in Appendix A.

######  Proposition 3.1.

The optimization problem (10) can be reformulated as an equivalent expression,
which is given by

$\displaystyle\underset{\parallel\bm{U}\parallel_{\infty}\leq\lambda\alpha}{min}\log\det\left(\frac{1}{\lambda(1-\alpha)}\left[\left(\mathcal{\bm{A}}^{2}+\lambda(1-\alpha)\bm{I}_{p}\right)^{1/2}-\mathcal{\bm{A}}\right]\right)$
$\displaystyle-\mathop{\rm
tr}\nolimits\left(\frac{1}{\lambda(1-\alpha)}\left[\left(\mathcal{\bm{A}}^{2}+\lambda(1-\alpha)\bm{I}_{p}\right)^{1/2}-\mathcal{\bm{A}}\right]\mathcal{\bm{A}}\right).$
(12)

where $\mathcal{\bm{A}}=\frac{1}{2}(\bm{S}+\bm{U}-\lambda(1-\alpha)\bm{T})$.

######  Corollary 3.2.

Since for any square nonsingular matrix $\bm{A}$ we have
$\log\det(\bm{A})=\mathop{\rm tr}\nolimits(\log(\bm{A}))$ so the optimization
problem (3.1) can also represented as follow

$\displaystyle\underset{\parallel\bm{U}\parallel_{\infty}\leq\lambda\alpha}{min}\mathop{\rm
tr}\nolimits(\log(\mathcal{\bm{B}})-\mathcal{\bm{B}}\mathcal{\bm{A}}),$ (13)

where
$\mathcal{\bm{B}}=\frac{1}{\lambda(1-\alpha)}[\left(\mathcal{\bm{A}}^{2}+\lambda(1-\alpha)\bm{I}_{p}\right)^{1/2}-\mathcal{\bm{A}}]$.

Also, we can express proposition 3.1 concerning eigenvalues of matrix
$\bm{S}+\bm{U}$. We summarize the result in the following proposition.

######  Proposition 3.3.

Let $\bm{S}+\bm{U}=\bm{P}\bm{D}\bm{P}^{-1}$ where $\bm{D}$ is a diagonal
matrix containing the eigenvalues of $\bm{S}+\bm{U}$ on the diagonal, and
$\bm{P}$ is a matrix that consists of the corresponding eigenvectors as its
columns. Also consider $\bm{T}=\gamma\bm{I}_{p}$ with $\gamma\in(0,\infty)$,
Substituting this decomposition, we can rewrite the optimization problem (3.1)
as follows

$\displaystyle\underset{\parallel\bm{U}\parallel_{\infty}\leq\lambda\alpha}{max}\sum_{i=1}^{p}\log\left(\sqrt{(b_{i})^{2}+\lambda(1-\alpha)}+b_{i}\right)+\frac{1}{\lambda(1-\alpha)}b_{i}\left(\sqrt{(b_{i})^{2}+\lambda(1-\alpha)}-b_{i}\right),$

where $b_{i}=1/2(\lambda_{i}-\lambda(1-\alpha)\gamma$ and
$\lambda_{1},\cdots,\lambda_{p}$ are eigenvalues of the matrix
$\bm{S}+\bm{U}$.

###  3.1 Conceptualizing the Estimator through Conjecture

Our problem statement starts from where we want to find an analytical solution
for Proposition 3.1. Indeed, the solution can be obtained numerically. Here,
we present an intuitive way to obtain a closed-form solution.
Noticing the proof of Proposition 3.1, in step one, we solve the objective
function subject to $||\bm{U}||_{\infty}\leq\lambda\alpha$ by setting the
tuning parameter equal to $\alpha\lambda$ resulting in the glasso. Then, in
step two, we substitute the glasso into the objective function and solve for
$\bm{\Theta}$. Subsequently with the following closed-form solution is
obtained

$\hat{\bm{\Theta}}^{2s}(\alpha,\lambda)=\left[\left(\lambda\bm{I}_{p}+\frac{1}{4}(\bm{W}-\lambda(1-\alpha)\bm{T})^{2}\right)^{1/2}+\frac{1}{2}(\bm{W}-\lambda(1-\alpha)\bm{T})\right]^{-1},$
(14)

where $\bm{W}$ is the covariance matrix estimate obtained from the glasso with
the tuning parameter $\alpha\lambda$. In Appendix, we have conducted a toy
simulation study to show how much the numerical solution is close to the
2-step ridge estimator $\hat{\bm{\Theta}}^{2s}$ in (14).

###  3.2 Detailed Development of the Estimator

The estimator presented in (14) can be achieved by considering the following
penalized log-likelihood:

$\log\det(\bm{\Theta})-(1-\lambda_{1})\mathop{\rm
tr}\nolimits(\bm{S}\bm{\Theta})-\lambda_{1}\mathop{\rm
tr}\nolimits(\bm{\Theta}\bm{\Gamma})-\frac{1}{2}\lambda_{2}||\bm{\Theta}-\bm{T}||_{F}^{2},$
(15)

where $\lambda_{1}\in(0,1]$ and $\lambda_{2}$ is a non-negative tuning
parameter. The penalized log-likelihood (15) is a combination of first
archetype ridge penalized log-likelihood, (2) and alternative ridge in (5).
Setting $\lambda_{1}=0$ yields the alternative ridge estimator, whereas
$\lambda_{2}=0$ esults in the Ledoit-Wolf estimator, as detailed in Equation
(3). To address the corresponding convex optimization problem, we consider the
normal equations derived from the derivative of Equation (15) with respect to
$\bm{\Theta}$,

$\bm{\Theta}^{-1}-(1-\lambda_{1})\bm{S}-\lambda_{1}\bm{\Gamma}-\lambda_{2}(\bm{\Theta}-\bm{T}).$
(16)

Following the method outlined in Kuismin (2017), we arrive at the following
solution for the estimator:

$\hat{\bm{\Theta}}=\left[\left(\lambda_{2}\bm{I}_{p}+\frac{1}{4}((1-\lambda_{1})\bm{S}+\lambda_{1}\bm{\Gamma}-\lambda_{2}\bm{T})^{2}\right)^{1/2}+\frac{1}{2}((1-\lambda_{1})\bm{S}+\lambda_{1}\bm{\Gamma}-\lambda_{2}\bm{T})\right]^{-1}.$
(17)

This result indicates an estimator similar to the Alternative Ridge, but with
a key difference: $\bm{S}$ in the equation is replaced with a convex
combination of $\bm{S}$ and $\bm{\Gamma}$. Notably, setting $\lambda_{1}=1$
and using the target matrix $\Gamma$ via the graphical lasso approach, we
achieve the estimator outlined in Equation (14). The estimator mentioned in
Equation (17), to the best of my knowledge, has not yet been explored in
existing literature. This opens a potential avenue for future research in this
domain. However, the primary focus of this paper remains on our proposed
estimator. In the subsequent section, we will illustrate the effective
performance of this estimator, underscoring its practical utility and
significance in the field.

##  4 Simulation Study

In this section, we evaluate and compare the performance of our proposed
2-step estimator $\hat{\bm{\Theta}}^{2s}$ given in (14) with other existing
methods, namely, graphical lasso or Glasso ([8]), type-I and II ridge or
Alternative ridge estimators ([10] and [16]), and the graphical elastic net or
Gelnet precision matrix estimator ([9]). The simulation structure used in our
evaluation is based on the framework proposed by [10].

We we generate simulated data from a multivariate Normal distribution
$\mathcal{N}_{p}(\mathbf{0},\bm{\Sigma})$, where
$\bm{\Sigma}=[\sigma_{k,k^{{}^{\prime}}}]$ and
$\bm{\Theta}=[\theta_{k,k^{{}^{\prime}}}]$ are positive definite matrices of
size $p\times p$. In our comparison, we consider six different models
representing network structures. The methods are tested on these structures,
considering both sparse and non-sparse scenarios. We specifically consider
network structures with $p=20,50,100$ nodes, and a sample size of $n=50$.

  * •

Network 1. One of the models used in the comparison is a compound symmetry
model, where $\sigma_{k,k}=1$ for all $k$, indicating equal variances across
variables, and $\sigma_{k,k^{\prime}}=0.6^{2}$ for $k\neq k^{\prime}$,
indicating a constant covariance between any pair of variables. This
covariance matrix exhibits a structured and non-sparse pattern.

  * •

Network 2. The precision matrix $\bm{\Theta}$ is derived from the prototype
matrix $\bm{\Theta}_{0}$ by standardizing it to have a unit diagonal. The
prototype matrix $\bm{\Theta}_{0}$ is defined as
$\bm{\Theta}_{0}=\bm{A}+a\bm{I}_{p}$, where each off-diagonal entry in
$\bm{A}$ is generated independently. Specifically, each off-diagonal entry is
set to $0.5$ with a probability of $0.1$, indicating a non-zero relationship,
or set to $0$ with a probability of $0.9$, indicating no relationship. The
parameter $a$ is chosen such that the condition number of the matrix is equal
to $p$, ensuring a well-conditioned matrix. This precision matrix exhibits
both unstructured and sparse patterns

  * •

Network 3. The precision matrix $\bm{\Theta}$ is estimated as
$\bm{\Theta}=\frac{1}{n}\bm{Y}^{T}\bm{Y}$, where $\bm{Y}=[y_{i,j}]$ is an
$n\times p$ matrix with $n=10,000$ observations and $p$ variables. Each
element $y_{i,j}$ is drawn independently from a standard normal distribution
$\mathcal{N}(0,1)$. The resulting precision matrix $\bm{\Theta}$ is an
estimate of the true precision matrix and represents the pairwise dependencies
among the variables. In this case, the precision matrix is considered
unstructured and non-sparse, as there are no specific patterns or constraints
imposed on its elements.

  * •

Network 4. The precision matrix $\bm{\Theta}$ follows a star model, where
$\theta_{k,k}=1$ for all variables. Additionally, there is a specific pattern
of connections between the first variable and all other variables, with
$\theta_{1,k}=\theta_{k,1}=0.1$, indicating a connection between the first
variable and the rest. All other off-diagonal elements $\theta_{k,k^{\prime}}$
are set to zero, implying that there are no connections between variables
other than the first variable. This results in a structured and sparse
precision matrix, as it exhibits a specific pattern of nonzero elements along
with many zero elements.

  * •

Network 5. The covariance matrix $\bm{\Sigma}$ follows a moving average (MA)
model, where each diagonal element $\sigma_{k,k}$ is equal to 1, indicating
that each variable has a unit variance. The off-diagonal elements follow a
specific pattern: $\sigma_{k,k-1}=\sigma_{k-1,k}=0.2$, representing a weaker
correlation between adjacent variables, and
$\sigma_{k,k-2}=\sigma_{k-2,k}=0.2^{2}$, representing an even weaker
correlation between variables that are two steps apart. All other off-diagonal
elements are set to zero, indicating no correlation beyond the specified
patterns. This results in a structured and sparse covariance matrix, as it
exhibits a specific pattern of nonzero elements along with many zero elements.

  * •

Network 6. The covariance matrix $\bm{\Sigma}$ follows a diagonally dominant
model. It is constructed based on a matrix $\bm{B}$ derived from another
matrix $\bm{A}$ by $\bm{B}=\frac{1}{2}(\bm{A}+\bm{A}^{T})$, where $\bm{A}$ is
a $p\times p$ matrix. Each off-diagonal element $a_{k,k^{{}^{\prime}}}$ (for
$k\neq k^{{}^{\prime}}$) in $\bm{A}$ is randomly drawn from a uniform
distribution $U(0,1)$, and the diagonal elements $a_{k,k}$ are set to zero. To
obtain the matrix $\bm{D}$, which is used to scale $\bm{B}$, each element of
$\bm{B}$ is divided by a scaling factor $\gamma$. The scaling factor $\gamma$
is determined as the largest row sum of the absolute values of the elements in
$\bm{B}$. Thus, $\bm{D}$ is computed as $\bm{D}=\frac{1}{\gamma}\bm{B}$,
resulting in a matrix $\bm{D}$ with values between 0 and 1. The final step
involves constructing the covariance matrix $\bm{\Sigma}$. The off-diagonal
elements $\sigma_{k,k^{{}^{\prime}}}$ are set to the corresponding elements
$d_{k,k^{{}^{\prime}}}$ from $\bm{D}$. The diagonal elements $\sigma_{k,k}$
are computed as $1+e_{i}$, where $e_{i}$ is randomly drawn from a uniform
distribution $U(0,0.1)$. This ensures that the diagonal elements have a
slightly larger value than 1, introducing some variability. Overall, this
construction results in an unstructured and non-sparse covariance matrix,
where the off-diagonal elements follow a specific pattern based on the values
in $\bm{D}$ and the diagonal elements have a slight perturbation.

To assess the performance of all methods, we conduct 100 independent
simulations for each network structure and calculate the average of various
loss functions. We determine the optimal tuning parameter $\lambda$ for each
method by performing five-fold cross-validation on each simulation run. For
the Glasso, Gelnet, and Alternative ridge methods, we utilize the five-fold
cross-validation provided by the R-package ”GLassoElnetFast” with $\lambda$
values ranging from $0$ to $10$ and $\alpha$ values ranging from $0$ to $1$.
As target matrices, we consider both the identity matrix $\bm{I}_{p}$ and a
scalar matrix $\nu\bm{I}_{p}$, where $\nu$ is calculated as
$\nu=\frac{p^{2}}{\text{tr}(\bm{S})}$.

###  4.1 Performance Measures

For the evaluation of the performance of a given estimator
$\hat{\bm{\Theta}}$, we consider three loss functions which have been widely
used (see, e.g., [2], [10], [4] and [9]). Here are the four loss functions:

  1. 1.

Kullback-Leibler loss (KL): This loss function measures based on the
likelihood function. It is defined as
$KL=\text{tr}(\bm{\Sigma}\hat{\bm{\Theta}})-\text{logdet}(\bm{\Sigma}\hat{\bm{\Theta}})-p$.

  2. 2.

L2 loss (L2): The L2 loss calculates the Frobenius norm of the difference
between the estimated precision matrix and the true precision matrix. It is
given by $L2=||\bm{\Theta}-\hat{\bm{\Theta}}||_{F}$.

  3. 3.

Quadratic loss (QL): It is defined as
$QL=\text{tr}((\bm{\Sigma}\hat{\bm{\Theta}}-\bm{I}_{p})^{2})$, where
$\bm{I}_{p}$ is the identity matrix of size $p$.

  4. 4.

Spectral norm loss (SP): The spectral norm loss measures the largest singular
value of the matrix $(\bm{\Theta}-\hat{\bm{\Theta}})^{2}$.

These loss functions provide quantitative measures to evaluate the performance
of the estimator by capturing different aspects of the dissimilarity between
the estimated precision matrix and the true precision matrix.

To compare the methods, the averages of the aforementioned loss functions are
computed based on $100$ simulations. Figures 1–4 present the final results of
the evaluation. In these figures, the columns represented by small black dots
represent the mean values of the loss functions, while the bars positioned on
top of each column indicate the standard errors (shown as mean ± SD). These
visual representations provide an overview of the performance of each method,
with the mean values indicating the average performance and the error bars
providing an indication of the variability or uncertainty associated with the
estimates. The plot layout is taken from [10].

Figure 1: The mean of the Kullback–Leibler loss for different models and
methods based on $100$ replications. Figure 2: The mean of $L2$ loss for
different models and methods based on $100$ replications. Figure 3: The mean
of quadratic loss for different models and methods based on $100$
replications. Figure 4: The mean quadratic loss for different models and
methods based on 100 replications.

The results for different networks are displayed in Figures 1-4. We summarize
some observations based on the results as follows:

  * •

In general, incorporating a target matrix in the estimation methods can be
advantageous because each approach tends to perform better when paired with
the appropriate target.

  * •

The Alternative ridge with target performs better than Glasso with no target
(as [10] mentioned before). Also, the 2-step ridge outperforms the Alternative
ridge in most considered measures and networks. (in Kullback–Leibler loss, the
2-step ridge shows better results than the alternative ridge in all networks.)

  * •

The effectiveness of different target matrices remains a question, but based
on our simulations, it is generally observed that the Identity target matrix
yields better results in most cases.

##  5 Application

In this section, we investigate the performance of our proposed estimator in
analysing three different datasets.

### The Ionosphere dataset

In the present section, our focus is on the ionosphere dataset (available at
http://archive.ics.uci.edu/ml/datasets/Ionosphere). Here, we analyze the
performance of linear discriminant analysis (LDA) by computing
misclassification rates. This analysis involves the utilization of various
estimates for covariance and precision matrices, alongside the implementation
of the Mahalanobis distance to partition samples into two distinct groups.
The dataset was gathered by a system located in Goose Bay, Labrador, focusing
on capturing information about free electrons present in the ionosphere. Radar
returns classified as ”good” exhibit noticeable structure within the
ionosphere, while ”bad” returns indicate signals that penetrate through the
ionosphere without any evident features. The dataset comprises a total of 351
observations and 35 variables, among which one variable is nominal with the
values ”good” or ”bad,” while the remaining 34 variables are continuous.
To compare the considered estimators, we adopted the methodology used by [10]
in their study on the impact of covariance matrix and precision matrix
estimations on the performance of Linear Discriminant Analysis (LDA) for
classification. They demonstrated the effect of different estimators on the
classification accuracy of LDA using this dataset. In our current study, we
aim to replicate their methodology and apply it to the estimators we are
considering. By doing so, we can evaluate how well our proposed estimator
performs in terms of its impact on LDA classification accuracy.

According to their work, the dataset consists of two groups: $G_{1}$
representing good radar returns and $G_{2}$ representing bad radar returns.
The population mean vectors for $G_{1}$ and $G_{2}$ are denoted as
$\bm{\mu}_{1}$ and $\bm{\mu}_{2}$, respectively. The precision matrix of the
population, denoted as $\bm{\Theta}$, is assumed to be the inverse of the
covariance matrix $\bm{\Sigma}$. The dataset is divided into a training set
and a test set. The training set consists of 40 observations, which is close
to the dimension of the data ($p=32$). The remaining 311 observations are
included in the test set. The goal is to estimate $\bm{\mu}_{1}$,
$\bm{\mu}_{2}$, and $\bm{\Theta}$. The sample means
$\overline{\mathbf{x}}_{1}$ and $\overline{\mathbf{x}}_{2}$ are used to
estimate $\bm{\mu}_{1}$ and $\bm{\mu}_{2}$, respectively. Different estimators
are employed to estimate $\bm{\Theta}$. Each observation $\mathbf{x}$ in the
test set is classified into either $G_{1}$ or $G_{2}$ using the LDA rule. The
classification is based on the squared Mahalanobis distance, which is computed
as $\mathbf{a}^{T}(\mathbf{x}-\bm{\mu})$, where
$\mathbf{a}=\bm{\Theta}(\bm{\mu}_{1}-\bm{\mu}_{2})$ and
$\bm{\mu}=0.5(\bm{\mu}_{1}+\bm{\mu}_{2})$. The misclassification rate is
estimated by applying the classification rule obtained from the training set
to the observations in the test set.

In our study, we consider Glasso, Gelnet, and Alternative ridge as alternative
estimators for the precision matrix $\Theta$. These estimators are used to
estimate $\Theta$, which is then utilized to allocate an observation
$\mathbf{x}$ to either group $G_{1}$ or $G_{2}$. To determine the group
allocation, we calculate the estimated discriminant
$\hat{\mathbf{a}}^{T}(\mathbf{x}-\hat{\bm{\mu}})$, where
$\hat{\mathbf{a}}=\hat{\Theta}(\overline{\mathbf{x}}_{1}-\overline{\mathbf{x}}_{2})$
and $\hat{\bm{\mu}}=0.5(\overline{\mathbf{x}}_{1}+\overline{\mathbf{x}}_{2})$.
If the calculated discriminant value is greater than zero, the observation
$\mathbf{x}$ is assigned to group $G_{1}$.

By employing these alternative estimators for $\Theta$, we aim to assess how
they influence the classification performance of the LDA model. This approach
allows us to compare the performance of the estimators and determine their
effectiveness in accurately allocating observations to the appropriate groups.
As target matrices, we used a zero matrix, identity matrix $I$ and a scalar
matrix $\nu I$, where $\nu=p^{2}/tr(S)$. The associated misclassification
rates are the mean values obtained from 100 iterations of the previously
indicated setting.

In Figure 5, we present the misclassification rates obtained for all the
estimators considered in our study, across different values of the tuning
parameter. The tuning parameter, denoted as $\rho$, is varied within a series
of length 50, ranging from $0.1$ to $0.9$. For the estimators Glasso and
Alternative ridge, the tuning parameter $\rho$ is directly applied. As for
Gelnet and 2-step ridge, we set $\alpha$ to $0.2$, and the value of the tuning
parameter $\lambda$ is set equal to $\rho$.

For further analysis, the data is randomly split into three sets: a training
set with $40$ observations, a validation set with $40$ more observations, and
a test set with the remaining observations. Cross-validation for Glasso,
Gelnet, Alternative ridge, and the 2-step ridge is used to choose the tuning
parameter from the validation set. The precision matrix estimates are then
derived using this tuning parameter value from the training set. Finally, we
compute the misclassification errors using the test set as our basis. This
process is repeated $100$ times. Figure 5 displays the boxplots of the $100$
misclassification errors made using this approach.

---
(a)
(b)
Figure 5:  Ionosphere dataset: (a) The misclassification rates of LDA and
standard errors (shaded areas)(b) boxplots of test misclassification errors

For the small value of $\rho$, 2-step ridge, Gelnet, and Alternative ridge
behave similarly and give a lower misclassification rate than Glasso, but when
$\rho$ increases, the Alternative ridge gives a lower rate than 2-step ridge,
and 2-step ridge gives a lower rate than Gelnet. The 2-step ridge has the
smallest standard error, slightly different from the Alternative ridge. When
the tuning parameter is chosen with cross-validation from the validation set
(5), 2-step ridge and Alternative ridge give the lowest median calculated for
the misclassification error: $0.164$ (zero target matrix), $0.164$ (Identity
target matrix), $0.162$ (scalar target matrix). In this example, the
Alternative ridge estimator gives a lower misclassification rate; considering
a small value for tuning parameter $\alpha$, the 2-step ridge can get close to
this rate. Now we consider another real example and show that our proposed
estimator has a lower misclassification rate.

### Wisconsin breast cancer dataset

This dataset is available at
https://archive.ics.uci.edu/ml/datasets/breast+cancer+wisconsin+(diagnostic)
and contains 569 observations and 31 variables, including one nominal variable
with the values ”malignant” or ”benign” and 30 continuous variables. The
problem is distinguishing cancerous, malignant, from non-cancerous, benign. We
follow the same strategy previously to compare the misclassification rate for
Alternative ridge and 2-step ridge estimators.

---
(a)
(b)
Figure 6: Wisconsin breast cancer dataset: The misclassification rates of LDA
and standard errors (shaded areas)when (a) $\alpha=0.2$. (b) $\alpha=0.5$.

Figure 6 shows the misclassification rates for Alternative ridge for each
tuning parameter value $\rho$ from a series of length $50$ ranging from $0.1$
to $0.9$ and 2-step ridge for $\alpha=0.2$ and $\lambda=\rho$. Figure 6
displays misclassification rates when $\alpha=0.5$. Figure 6 shows the
misclassification rates of the 2-step ridge are lower than the Alternative
ridge estimators in this example for both consider tunning parameter $\alpha$.

### Top 100 cryptocurrencies historical dataset

Here, we consider the top 100 cryptocurrencies’ historical dataset, which is
available on Kaggle. The dataset contains the open, high, low, close and
volume values of the cryptocurrencies priced based on USD. We examine a total
of 30 cryptocurrencies over a span of three years, specifically from 2019 to
2021. The names and corresponding abbreviations of these 30 cryptocurrencies
can be found in Table 1.

Table 1: List of abbreviations of the considered cryptocurrencies. Cryptocurrencies |  |
---|---|---
ADA | Cardano | EOS | Electro-Optical System | NEM | New Economy Movement
BAT | Basic Attention Token | ETC | Ethereum Classic | QTUM | Qtum
BCH | Bitcoin Cash | ETH | Ethereum | TRX | TRON
BNB | Binance Coin | FIL | Filecoin | USDT | Tether
BTC | Bitcoin | GNO | Gnosis | WAVES | WAVES
BTG | Bitcoin Gold | IOTA | IOTA | XLM | Stellar
DASH | Digital Cash | KCS | KuCoin Shares | XMR | Monero
DCR | Decred | LRC | Loopring | XRP | Ripple
DOGE | Dogecoin | LTC | Litecoin | XTZ | Tezos
ENJ | Enigma | MANA | Decentraland | ZEC | ZCash

We compute the daily log returns of the closed prices for the remaining 30
cryptocurrencies.

We utilize our proposed estimator to estimate the networks using both three-
year and annual data separately. In the estimation process, the Identity
matrix is used as the target. To select the tuning parameters, we employ five-
fold cross-validation. The estimated networks obtained from these estimations
are presented in Figure 7. The edges that correspond to negative partial
correlations are indicated in orange, and the edges that belong to positive
partial correlations are indicated in blue. The width of the edge in the
network visualization is proportional to the magnitude of the partial
correlation estimates.

|
---|---
(a) Estimated network based on 2019 | (b) Estimated network based on 2020
|
(c) Estimated network based on 2021 | (d) Estimated network based on 2019-2021
Figure 7: Estimated cryptocurrency networks based on 2019 to 2021 using 2-step ridge. |
---|---
(a) Ascending-descending trend of partial correlations | (b) Stable relations based on 2019-2021
Figure 8: Ascending-descending trend and Stable relations.

Figure 8 shows the increasing (in green) and decreasing (in red) trend of
partial correlations. Figure 8 shows the stable relationships between all
dependencies based on a three-year estimate. So, for instance, the partial
correlation between Bitcoin and Ethereum is increasing across 2019, 2020, and
2021, while the partial correlation between Bitcoin and Bitcoin Cash is
decreasing during the same time. Also, about 8, for example, the connection
between Bitcoin and Monero is stable across 2019, 2020 and 2021, while the
connection between Bitcoin and Digital Cash is not stable because it is only
detected in 2019.
As a network measure, we provide the mean strength value. The strength of a
node is equal to the sum of the weights (partial correlations) of all
connected edges. This measure shows the intensity or strength of the
connections between currencies. Figure 9 presents the mean value of the
strength estimated using a rolling window of one year, with shifts of one
month. This rolling window analysis allows us to capture the changes in
strength over time and observe any patterns or trends in the data.

Figure 9: The mean value of strength estimated on a rolling window of one year
with shifts of one month.

We observe a decrease in the mean values of strength between 2019 and 2020,
followed by an increase in the first month of 2021 and a more noticeable
increase in the latest months of 2021. The high mean strength values in 2021
can explain the similarity of the estimated network 2021 to the three-year
estimated network.

##  6 Conclusions

Including a reasonable target matrix seems to improve the estimation. It is
also mentioned by [9]. Furthermore, we can not choose one estimator as the
best. According to the simulation study and real data examples, the better
performance of the 2-step ridge precision estimator is visible than the
alternative ridge estimator with corresponding target matrices. However, the
alternative ridge is superior to our suggested estimate in terms of
computational time. Finally, We consider two real data examples to compare the
proposed estimator to the alternative ridge. Followed by an empirical
examination of the network across 30 cryptocurrencies. We consider the daily
closed prices and the 2-step ridge to estimate their interactions and, in
turn, the perceived network.

## Acknowledgement

This work was based upon research supported in part by the National Research
Foundation of South Africa Ref.: SRUG2204203965, and RA171022270376, UID
119109; the Department of Research and Innovation at the University of
Pretoria (SA), STATOMET, as well as the Centre of Excellence in Mathematical
and Statistical Sciences, based at the University of the Witwatersrand (SA).
The opinions expressed and conclusions arrived at are those of the authors and
are not necessarily to be attributed to the NRF. The second author’s research
(M. Arashi) is supported by a grant from Ferdowsi University of Mashhad
(N.2/58091).

## References

  * [1] M. Arashi, A.M.E. Saleh, and B.G. Kibria, _Theory of ridge regression estimation with applications_ , John Wiley & Sons, 2019.
  * [2] V. Avagyan, A.M. Alonso, and F.J. Nogales, _Improving the graphical lasso estimation for the precision matrix through roots of the sample covariance matrix_ , Journal of Computational and Graphical Statistics 26 (2017), pp. 865–872.
  * [3] O. Banerjee, L. El Ghaoui, and A. d’Aspremont, _Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary data_ , The Journal of Machine Learning Research 9 (2008), pp. 485–516.
  * [4] D. Bernardini, S. Paterlini, and E. Taufer, _New estimation approaches for graphical models with elastic net penalty_ , arXiv preprint arXiv:2102.01053 (2021).
  * [5] J. Bien and R.J. Tibshirani, _Sparse estimation of a covariance matrix_ , Biometrika 98 (2011), pp. 807–820.
  * [6] P. Danaher, P. Wang, and D.M. Witten, _The joint graphical lasso for inverse covariance estimation across multiple classes_ , Journal of the Royal Statistical Society: Series B (Statistical Methodology) 76 (2014), pp. 373–397.
  * [7] J. Fan, Y. Feng, and Y. Wu, _Network exploration via the adaptive lasso and scad penalties_ , The annals of applied statistics 3 (2009), p. 521.
  * [8] J. Friedman, T. Hastie, and R. Tibshirani, _Sparse inverse covariance estimation with the graphical lasso_ , Biostatistics 9 (2008), pp. 432–441.
  * [9] S. Kovács, T. Ruckstuhl, H. Obrist, and P. Bühlmann, _Graphical elastic net and target matrices: Fast algorithms and software for sparse precision matrix estimation_ , arXiv preprint arXiv:2101.02148 (2021).
  * [10] M. Kuismin, J. Kemppainen, and M. Sillanpää, _Precision matrix estimation with rope_ , Journal of Computational and Graphical Statistics 26 (2017), pp. 682–694.
  * [11] O. Ledoit and M. Wolf, _A well-conditioned estimator for large-dimensional covariance matrices_ , Journal of multivariate analysis 88 (2004), pp. 365–411.
  * [12] R. Mazumder and T. Hastie, _The graphical lasso: New insights and alternatives_ , Electronic journal of statistics 6 (2012), p. 2125.
  * [13] J. Schäfer and K. Strimmer, _A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics_ , Statistical applications in genetics and molecular biology 4 (2005).
  * [14] A. Shirilord and M. Dehghan, _Closed-form solution of non-symmetric algebraic riccati matrix equation_ , Applied Mathematics Letters 131 (2022), p. 108040.
  * [15] W.N. van Wieringen, _The generalized ridge estimator of the inverse covariance matrix_ , Journal of Computational and Graphical Statistics 28 (2019), pp. 932–942.
  * [16] W.N. Van Wieringen and C.F. Peeters, _Ridge estimation of inverse covariance matrices from high-dimensional data_ , Computational Statistics & Data Analysis 103 (2016), pp. 284–303.
  * [17] D.I. Warton, _Penalized normal likelihood and ridge regularization of correlation and covariance matrices_ , Journal of the American Statistical Association 103 (2008), pp. 340–349.
  * [18] D.M. Witten, J.H. Friedman, and N. Simon, _New insights and faster computations for the graphical lasso_ , Journal of Computational and Graphical Statistics 20 (2011), pp. 892–900.
  * [19] D.M. Witten and R. Tibshirani, _Covariance-regularized regression and classification for high dimensional problems_ , Journal of the Royal Statistical Society: Series B (Statistical Methodology) 71 (2009), pp. 615–636.
  * [20] M. Yuan and Y. Lin, _Model selection and estimation in the gaussian graphical model_ , Biometrika 94 (2007), pp. 19–35.

## Appendix A: Proofs

Lemma A1(Theorem 3.1 of [14]): Let $\bm{A},\bm{B},\bm{C},\bm{D}$ be real or
complex value $n\times n$ matrices. Suppose that $\bm{C}$ is nonsingular and
the matrices $\bm{A},\bm{B}$ and $\bm{C}$ satisfy $\bm{B}\bm{C}=\bm{C}\bm{A}$.
Denote $\bm{\Delta}:=\bm{A}^{2}+\bm{D}\bm{C}$, be a nonsingular matrix. Then a
closed-form expression for the solution of a non-symmetric algebraic Riccati
equation, $\bm{A}\bm{X}+\bm{X}\bm{B}+\bm{X}\bm{C}\bm{X}=\bm{D}$ is given by
the following formula

$\bm{X}=(\sqrt{\bm{\Delta}}-\bm{A})\bm{C}^{-1},$ (18)

where $\sqrt{\bm{M}}$ is square root of matrix $\bm{M}$.

Proof of Proposition 3.1: Consider the dual optimization of problem (10) by
exchanging the max and the min. The resulting inner problem in $\bm{\Theta}$
can be solved analytically by setting the gradient of the objective to zero
and solving for $\bm{\Theta}$, we have

$\displaystyle\bm{\Theta}^{-1}-\bm{S}-\bm{U}+\lambda(1-\alpha)(\bm{\Theta}-\bm{T})=\bm{0},$
$\displaystyle\bm{\Theta}^{-1}-\bm{S}-\bm{U}+\lambda(1-\alpha)\bm{T}-\lambda(1-\alpha)\bm{\Theta}=\bm{0}.$
(19)

Multiplying the LHS of (Appendix A: Proofs) by $\bm{\Theta}$, we get

$(\bm{S}+\bm{U}-\lambda(1-\alpha)\bm{T})\bm{\Theta}+\lambda(1-\alpha)\bm{\Theta}^{2}=\bm{I}_{p}.$
(20)

We use Lemma A1 to solve the quadratic matrix equation (20), therefore we
consider the matrices as follows;
$\bm{A}=\bm{B}=\frac{1}{2}(\bm{S}+\bm{U}-\lambda(1-\alpha)\bm{T})$,
$\bm{C}=\lambda(1-\alpha)\bm{I}$ and $\bm{D}=\bm{I}$, we have

$\bm{\Delta}=(\frac{1}{2}(\bm{S}+\bm{U}-\lambda(1-\alpha)\bm{T}))^{2}+\lambda(1-\alpha)\bm{I}.$
(21)

Since $\bm{C}$ and $\bm{\Delta}$ are non-singular and
$\bm{B}\bm{C}=\bm{C}\bm{A}$ then the solution of equation (20) is given by
following formula

$\bm{\Theta}^{*}=\frac{1}{\lambda(1-\alpha)}\left[\left(\mathcal{\bm{A}}^{2}+\lambda(1-\alpha)\bm{I}\right)^{1/2}-\mathcal{\bm{A}})\right],$
(22)

where $\mathcal{\bm{A}}=\frac{1}{2}(\bm{S}+\bm{U}-\lambda(1-\alpha)\bm{T})$.
Subsequently, we can rewrite (10) as

$\underset{\parallel\bm{U}\parallel_{\infty}\leq\lambda\alpha}{min}\log\det(\bm{\Theta}^{*})-tr(\bm{\Theta}^{*},\bm{S}+\bm{U})+\frac{(1-\alpha)}{2}||\bm{\Theta}^{*}-T||_{2}^{2}.$
(23)

The second and third parts of (23) can be expressed of the following form

$\displaystyle tr(\bm{\Theta}^{*}(\bm{S}+\bm{U}))$ $\displaystyle+$
$\displaystyle\frac{(1-\alpha)}{2}||\bm{\Theta}^{*}-\bm{T}||_{2}^{2}$ (24)
$\displaystyle=$ $\displaystyle
tr(\bm{\Theta}^{*}(2\mathcal{\bm{A}}+\frac{\lambda(1-\alpha)}{2}\bm{\Theta}^{*}))+\frac{\lambda(1-\alpha)}{2}tr(\bm{T}^{T}\bm{T}),$

with substituting (22) in (24) and doing some calculations, we have

$\displaystyle\mathop{\rm
tr}\nolimits(\bm{\Theta}^{*}(2\mathcal{\bm{A}}+\frac{\lambda(1-\alpha)}{2}\bm{\Theta}^{*}))=\frac{1}{\lambda(1-\alpha)}tr((\mathcal{\bm{A}}^{2}+\lambda(1-\alpha)\bm{I})^{1/2}-\mathcal{\bm{A}})\mathcal{\bm{A}})+\frac{p}{2}.$
(25)

Therefore the optimization problem (23) is equivalent to

$\displaystyle\underset{\parallel\bm{U}\parallel_{\infty}\leq\lambda\alpha}{\min}\log\det\left(\frac{1}{\lambda(1-\alpha)}\left[\left(\mathcal{\bm{A}}^{2}+\lambda(1-\alpha)\bm{I}\right)^{1/2}-\mathcal{\bm{A}}\right]\right)$
$\displaystyle-\mathop{\rm
tr}\nolimits\left(\frac{1}{\lambda(1-\alpha)}\left[\left(\mathcal{\bm{A}}^{2}+\lambda(1-\alpha)\bm{I}\right)^{1/2}-\mathcal{\bm{A}}\right]\mathcal{\bm{A}}\right).$

which completes the proof.

Proof of Proposition 3.3: Since $\bm{S}+\bm{U}=\bm{P}\bm{D}\bm{P}^{-1}$ and
$\bm{T}=\gamma\bm{I}_{p}$ we can rewrite $\mathcal{\bm{A}}$ in equation (3.1)
as
$\mathcal{\bm{A}}=\bm{P}[\bm{D}-\lambda(1-\alpha)\gamma\bm{I}_{p}]\bm{P}^{-1}$.
After some calculations, we have

$\displaystyle\underset{\parallel\bm{U}\parallel_{\infty}\leq\lambda\alpha}{\max}\mathop{\rm
tr}\nolimits\left(\bm{P}\left[\log(\sqrt{\bm{B}^{2}+\lambda(1-\alpha)\bm{I}}+\bm{B})+\frac{1}{\lambda(1-\alpha)}\bm{B}(\sqrt{\bm{B}^{2}+\lambda(1-\alpha)\bm{I}}-\bm{B})\right]\bm{P}^{-1}\right),$

where $\bm{B}=1/2(\bm{D}-\lambda(1-\alpha)\bm{T})$ and according to this fact
that the sum of the eigenvalues of any matrix is the same as the trace of a
matrix, the proof is complete.

## Appendix B: A toy simulation study

For the given conjecture in (14), we simulate the data from a multivariate
Gaussian distribution $\mathcal{N}_{p}(\bm{\mu},\bm{\Sigma})$, where
$\bm{\mu}$ is zero vector and $\bm{\Sigma}=[\sigma_{k,k^{{}^{\prime}}}]$ with
$\sigma_{k,k}=1$ and $\sigma_{k,k^{{}^{\prime}}}=0.6^{2}$ for $k\neq
k^{{}^{\prime}}$. Then we compare the solution of two optimization problems
(3.1) and $\max\log\det(\bm{S}+\bm{U})$ subject to
$||\bm{U}||_{\infty}\leq\lambda\alpha$ for some sample covariance matrices
simulated from multivariate Gaussian distribution. It should be noted for the
case $p=2$. They have the same solution. Two objective functions for
$\alpha=0.4$ and $\lambda=0.6$ plot in Figure 10, and it is clear they have
the same optimal point.

Figure 10: Objective functions for two considered optimization problems where
the green one is related to problem (3.1). $\alpha=0.4$ and $\lambda=0.5$ are
considered, and two vertical lines show the acceptable area.

For $p=3$, if we consider $\alpha=0.4$ and $\lambda=0.6$ and matrix
$U=[u_{k,k^{{}^{\prime}}}]$, Table 2 shows the difference between the
estimated non-diagonal elements $u_{12},u_{13}$ and $u_{23}$ of the solution
of two optimization problems for $15$ iterations. The differences show the
results are very close to zero, as we expect in the proposed conjecture.

Table 2: The difference between the non-diagonal elements of solution of two optimization problems $Iteration$ | $p=3$
---|---
1 | 4.191486e-11, 1.373684e-08, 5.860205e-11
2 | 2.810245e-09, 1.942890e-16, 2.380667e-09
3 | 1.134227e-06, 6.891155e-07, 2.775558e-17
4 | 1.354186e-10, 1.921179e-08, 8.326673e-17
5 | 1.897935e-07, 1.223164e-07, 0.000000e+00
6 | 1.767697e-07, 3.756689e-08, 3.608225e-16
7 | 2.021518e-07, 6.238265e-08, 4.440892e-16
8 | 6.032975e-08, 2.808495e-08, 1.665335e-16
9 | 1.342834e-06, 7.493853e-07, 5.551115e-17
10 | 8.326673e-17, 7.227610e-10, 9.018657e-10
11 | 5.024786e-08, 5.419546e-06, 3.608225e-16
12 | 4.345127e-07, 1.283257e-07, 5.551115e-17
13 | 3.330669e-16, 8.649000e-10, 2.368810e-10
14 | 5.007422e-11, 2.035536e-09, 4.718448e-16
15 | 4.617334e-12, 8.790247e-03, 1.942890e-16


# Long-time relaxation dynamics in nematic and smectic liquid crystals of
soft-repulsive colloidal rods

Daniela Cywiak<EMAIL_ADDRESS>División de Ciencias e Ingenierías,
Universidad de Guanajuato, Campus León, México  Alessandro Patti
<EMAIL_ADDRESS>Department of Chemical Engineering and
Analytical Science, University of Manchester, United Kingdom  Alejandro Gil-
Villegas<EMAIL_ADDRESS>División de Ciencias e Ingenierías, Universidad
de Guanajuato, Campus León, México

(August 28, 2024)

###### Abstract

Understanding the relaxation dynamics of colloidal suspensions is crucial to
identify the elements that influence the mobility of their constituents,
assess their macroscopic response across the relevant time and length scales,
and thus disclose the fundamentals underpinning their exploitation in
formulation engineering. In this work, we specifically assess the impact of
long-ranged ordering on the relaxation dynamics of suspensions of soft-
repulsive rod-like particles, which are able to self-organise into nematic and
smectic liquid-crystalline phases. By performing Dynamic Monte Carlo
simulations, we analyse the effect of translational and orientational order on
the diffusion of the rods along the relevant directions imposed by the
morphology of the background phases. To provide a clear picture of the
resulting dynamics, we assess their dependence on temperature, which can
dramatically determine the response time of the system relaxation and the
self-diffusion coefficients of the rods. The computation of the van Hove
correlation functions allows us to identify the existence of rods that diffuse
significantly faster than the average and whose concentration can be
accurately adjusted by a suitable choice of temperature.

††preprint: APS/123-QED

##  I Introduction

The study of the dynamical properties of colloidal suspensions has generated
an exceptional insight into multiple fields of science and technology,
shedding light on processes such as self-assembly and nucleation Sacanna  _et
al._ (2013a), and systems, such as active matter Huber  _et al._ (2018) and
liquid crystals (LCs) Fraden  _et al._ (1989); Dogic and Fraden (2006); Pouget
_et al._ (2011), whose key mechanisms are dramatically influenced by the
particle dynamics over time. The relatively recent development of experimental
techniques for the synthesis of a fascinating spectrum of particle shapes
Sacanna  _et al._ (2013b) has boosted the interest in extending to anisotropic
particles methods predominantly used to study dynamic arrest and glass
transition in suspensions of spherical particles Murphy  _et al._ (2020).
These methods include the use of hard-core models to mimic the essential
features of a spectrum of colloidal systems, including LCs. As Onsager
demonstrated in his seminal theory, these relatively simple models are indeed
able to predict the isotropic-to-nematic transition of infinitely long rod-
like particles, as later on verified by computer simulation Frenkel (1988);
McGrother  _et al._ (1996). Nevertheless, Van der Waals and long-range
electrostatic forces modify the phase diagram and introduce important effects
in the thermodynamic and structural properties of these systems, particularly
in the location of phase boundaries, such as isotropic-nematic and nematic-
smectic transitions McGrother  _et al._ (1998); van Duijneveldt  _et al._
(2000).

A wide variety of hard- and soft-core models have been used to gain a better
insight into the behavior of thermotropic (molecular) and colloidal LCs Allen
(2019). In the case of spherical particles, Medina-Noyola and coworkers
demonstrated that structural properties and diffusion coefficients of spheres
interacting with a repulsive Sutherland potential can be mapped onto the
corresponding properties of hard-sphere (HS) systems, using effective
diameters dependent on density and temperature Guevara-Rodríguez and Medina-
Noyola (2003); Ramírez-González  _et al._ (2011). More recently, Jackson and
coworkers obtained an effective soft-repulsive potential based on the Mie
model in order to reproduce structural properties of the HS system Jover  _et
al._ (2012), a procedure that has been extended to square-well potentials
Zerón  _et al._ (2018) and applied successfully in molecular simulations of
colloidal systems Báez  _et al._ (2018). The approach of using potential
models with soft repulsive and attractive interactions of variable range, such
as the Mie model for chain molecules Davies  _et al._ (1998); Lafitte  _et
al._ (2013), have had important implications in the prediction of a wide
variety of phase diagrams for molecular fluids Dufal  _et al._ (2015).

When investigating the dynamics of hard-core particles, the usual simulation
techniques of choice, namely Molecular dynamics (MD) and Brownian dynamics
(BD), cannot be directly applied, as the form of the interaction potential
does not allow one to integrate the equations of motion. Additionally, while
MD can accurately reproduce the deterministic dynamics of atoms and molecules,
it is not suitable to mimic the Brownian motion of colloids unless both
dispersed and continuous phases are explicitly incorporated, resulting in a
very computationally demanding system to study. On the other hand, BD relies
on stochastic equations of motion that implicitly incorporate the presence of
a solvent, but faces, along with MD, the limitation of an integration time-
step that should be sufficiently small to guarantee precision at the cost of
missing the long-time relaxation dynamics. This is not a major problem in
sufficiently dilute colloidal suspensions, but it is indeed a challenge in
dense colloids, including LCs, whose structural relaxation completely unfolds
over relatively long time scales. These limitations can be bypassed by the
Dynamic Monte Carlo (DMC) method, which is able to capture the Brownian
dynamics of colloids without employing stochastic or deterministic equations
of motion. Based on the standard Metropolis algorithm Metropolis  _et al._
(1953), DMC can quantitatively and qualitatively reproduce BD simulation
results in the limit of small displacements Sanz and Marenduzzo (2010); Romano
_et al._ (2011); Patti and Cuetos (2012); Cuetos and Patti (2015); Corbett
_et al._ (2018); Chiappini  _et al._ (2020); García Daza  _et al._ (2020);
García Daza  _et al._ (2022).

In this work, we apply the DMC method to study the dynamics of nematic (N) and
smectic (Sm) LCs of colloidal rod-like particles. Rods are here represented as
soft repulsive spherocylindrical (SRS) particles interacting via the Kihara
potential, which has been used in the past to investigate the behavior of
prolate Cuetos  _et al._ (2005); Cuetos and Martínez-Haya (2015) and oblate
Cuetos  _et al._ (2002) spherocylinders. We are specifically interested in
assessing the effect of temperature on the long-time relaxation dynamics of
SRS particles in N and Sm phases. While the former are characterised by a mere
orientational order with the rods approximately aligned along a common
direction, but randomly positioned, the latter also exhibit translational
ordering with the rods arranged in contiguous layers parallel to each other.

This paper is organised as follows. In Section II, we describe the SRS model,
provide the main details of the DMC simulation method applied in this work and
introduce the observables calculated to characterize the dynamics in N and Sm
LCs and its dependence on long-range ordering and temperature. In Section III,
we discuss the ability of SRS to diffuse in these LC phases and the structural
relaxation of the systems over time. Finally, we draw our conclusions in
Section IV.

##  II Model and simulations

As mentioned above, colloidal rods are here modelled as SRS particles
interacting via the Kihara potential. In particular, we studied systems of
$N_{r}=1000$ prolate spherocylinders of length-to-diameter ratio
$L^{\star}\equiv L/\sigma=5$, where $L$ and $\sigma$ are, respectively, the
length and diameter of a cylindrical body capped by two hemispheres of
identical diameter. The total length of the rod is thus $L+\sigma$. The
interaction between particles is described by a repulsive, truncated and
shifted Kihara potential that reads

$U_{ij}=\left\\{\begin{array}[]{ll}4\epsilon\left[\left(\frac{\sigma}{d_{m}}\right)^{12}-\left(\frac{\sigma}{d_{m}}\right)^{6}+\frac{1}{4}\right]&\quad
d_{m}\leq\sqrt[6]{2}\sigma\\\ \\\ 0&\quad
d_{m}>\sqrt[6]{2}\sigma\end{array},\right.$ (1)

where $U_{ij}=U_{ij}(\textbf{r}_{ij},\hat{\bf u}_{i},\hat{\bf u}_{j})$. The
subscripts $i$ and $j$ refer to a pair of interacting spherocylinders,
$\textbf{r}_{ij}$ is the centre-to-centre distance between them, $\epsilon$
their interaction strength, $\hat{\bf u}_{i}$ and $\hat{\bf u}_{j}$ indicate
the particles’ orientation, and $d_{m}=d_{m}(\textbf{r}_{ij},\hat{\bf
u}_{i},\hat{\bf u}_{j})$ is the minimum distance between $i$ and $j$. The
interested reader is referred to Ref. Vega and Lago (1994) for additional
details on the computation of the minimum distance between prolate
spherocylinders. All simulations have been performed in the canonical ensemble
at constant number of particles, temperature ($T$) and volume ($V$). The
density, $\rho$, of N and Sm phases was set according to the phase diagrams of
Kihara spherocylinders reported in Ref. Cuetos and Patti (2015). In
particular, $\rho_{\rm N}^{\star}=0.12$ and $\rho_{\rm Sm}^{\star}=0.15$ for N
and Sm phases, respectively, with $\rho^{\star}=N_{r}\sigma^{3}/V$. To
equilibrate N and Sm phases at these values of density, we melted initial
configurations of perfect crystals at the desired temperature by performing
standard MC simulations. In particular, rotational and translational movements
of randomly selected particles were accepted or rejected according to the
Metropolis algorithm Metropolis  _et al._ (1953). The systems were considered
to be at equilibrium when the potential energy of the system had reached a
steady value within moderate statistical fluctuations.

After equilibration, we performed DMC simulations to study the relaxation
dynamics of N and Sm phases. In DMC simulations, to realistically mimic the
Brownian motion of colloidal particles, unphysical moves, such as swaps or
cluster moves, are not performed. An insightful description of the DMC method
applied in this work is available elsewhere Patti and Cuetos (2012); Cuetos
and Patti (2015); Corbett  _et al._ (2018); García Daza  _et al._ (2020);
Chiappini  _et al._ (2020); García Daza  _et al._ (2022). Here we only provide
a brief overview of the main features of the method and refer the interested
reader to these works for details. DMC simulations were performed in N and Sm
phases at scaled temperatures $T^{\star}=5,8,10,12,15$ and 20, with
$T^{\star}=k_{B}T/\epsilon$ and $k_{B}$ the Boltzmann constant. We selected
$\sigma$, $\epsilon$ and $\tau=\sigma^{2}/D_{0}$ as units of length, energy
and time, respectively, where $D_{0}=k_{B}T/\mu\sigma$ is a diffusion constant
and $\mu$ the viscosity of the implicit solvent. One DMC cycle consists of
$N_{r}$ attempts of simultaneously displacing and rotating a randomly selected
particle. These moves are accepted according to the probability
$\min[1,\exp(-\Delta U/k_{B}T)]$, where $\Delta U$ is the energy difference
between new and old configurations. The magnitude of elementary displacements
and rotations is defined according to the Einstein equations and the particle
diffusion coefficients at infinite dilution. More specifically, displacements
in the direction of $\hat{\bf u}_{i}$ and perpendicular to it where randomly
selected from uniform distributions that satisfy the conditions
$|X_{\parallel}|\leq\sqrt{2D_{\parallel}\delta t_{\rm MC}}$ and
$|X_{\perp}|\leq\sqrt{2D_{\perp}\delta t_{\rm MC}}$, respectively, where
$\delta t_{\rm MC}$ is the arbitrarily set MC time step, while $D_{\parallel}$
and $D_{\perp}$ correspond to the translational diffusion coefficients at
infinite dilution along $\hat{\bf u}_{i}$ and perpendicularly to it,
respectively. A similar approach was used to calculate particle rotations,
where the vector $\hat{\bf u}_{i}$ changes to $\hat{\bf u}_{i}+\delta\hat{\bf
u}_{i}$, with $\delta\hat{\bf u}_{i}=Y_{\varphi,1}\hat{\bf
w}_{j,1}+Y_{\varphi,2}\hat{\bf w}_{j,2}$. The two randomly chosen vectors
$\hat{\bf w}_{j,m}$ are perpendicular to each other and to $\hat{\bf u}_{i}$.
Finally, $Y_{\varphi,1}$ and $Y_{\varphi,2}$ are random numbers selected from
uniform distributions that satisfy the condition
$|Y_{\varphi,m}|\leq\sqrt{2D_{\varphi}\delta t_{\rm MC}}$ where $D_{\varphi}$
is the rotational self-diffusion coefficient of the rod.

Translational and rotational diffusion coefficients of rod-like particles at
infinite dilution have been estimated by applying the analytical expressions
based on the induced-forces method by Bonet Avalos and coworkers Bonet Avalos
_et al._ (1994). In particular:

$\displaystyle\frac{D_{\perp}}{D_{0}}$
$\displaystyle=\frac{\ln(2/\gamma)-1/2-I^{tt}}{2\pi/\gamma}$ (2)
$\displaystyle\frac{D_{\parallel}}{D_{0}}$
$\displaystyle=\frac{\ln(2/\gamma)-3/2-I^{tt}}{\pi/\gamma}$ (3)
$\displaystyle\frac{D_{\varphi}}{D_{0}}$
$\displaystyle=3\frac{\ln(2/\gamma)-11/6-I^{rr}}{\pi\sigma^{2}/\left(2\gamma\right)^{3}}$
(4)

where $1/\gamma=2\left(L^{\star}+1\right)$,
$I^{tt}\equiv\frac{1}{2}\int_{-1}^{1}dx\ln h(x)$,
$I^{rr}\equiv\frac{3}{2}\int_{-1}^{1}dxx^{2}\ln h(x)$ and
$h(x)=(1-x^{2n})^{1/2n}$. The function $h(x)$, applied to model particles with
revolution symmetry, approximates very well the shape of a spherocylinder when
$n=8$. Under these conditions, $I^{tt}\simeq-0.0061$ and $I^{rr}\simeq-0.017$.
The translational and rotational diffusion coefficients at infinite dilution
used in our DMC simulations are given in Table I for different values of the
reduced temperature. We note that the value of the three sets of diffusion
coefficients are different from those employed in Ref. Patti and Cuetos
(2012), especially so at $T^{\star}=20$. Such difference is due to the fact
that Ref. Patti and Cuetos (2012) made use of the analytical expressions
available for prolate spheroids, whose geometry is roughly similar, but not
identical to that of spherocylinders.

Table 1: Translational and rotational diffusion coefficients at infinite dilution as obtained from Eqs. (2-4) at different values of the reduced temperature. $T^{\star}$ | $D_{\parallel}/{D_{0}}$ | $D_{\perp}/{D_{0}}$ | $D_{\varphi}/{D_{0}}$
---|---|---|---
5 | 0.223 | 0.177 | 0.029
8 | 0.357 | 0.284 | 0.047
10 | 0.446 | 0.355 | 0.059
12 | 0.536 | 0.427 | 0.071
15 | 0.670 | 0.533 | 0.089
17 | 0.759 | 0.605 | 0.101
20 | 0.893 | 0.711 | 0.119

In order to obtain the actual time scale for Brownian dynamics, we have
rescaled the MC timescale via the acceptance rate using the following relation

$t_{\rm BD}=\frac{\mathcal{A}}{3}C_{\rm MC}\delta t_{\rm MC}$ (5)

where $t_{\rm BD}$ is the Brownian time after $C_{\rm MC}$ cycles,
$\mathcal{A}$ is the average acceptance rate over this number of cycles and
$\delta t_{MC}$ the corresponding MC time. This result allows one to rescale
the dynamical properties and consistently compare the structural relaxation of
N and Sm phases across the spectrum of temperatures studied here. More
specifically, we evaluated the mean-squared displacement (MSD), the self-part
of the Van Hove correlation function (s-VHF) and the self-part of the
intermediate scattering function (s-ISF). Each of these functions has been
calculated in the direction of the nematic director and in planes
perpendicular to it. In particular, the parallel and perpendicular MSD read

$\big{\langle}\Delta
r^{2}_{\parallel}(t)\big{\rangle}=\frac{1}{N_{r}}\bigg{\langle}\sum[(\mathbf{r}_{j}(t)-\mathbf{r}_{j}(0))\cdot\mathbf{n}]^{2}\bigg{\rangle}$
(6)

$\big{\langle}\Delta
r^{2}_{\perp}(t)\big{\rangle}=\frac{1}{N_{r}}\bigg{\langle}\sum[(\mathbf{r}_{j}(t)-\mathbf{r}_{j}(0)])\times\mathbf{n}]^{2}\bigg{\rangle}$
(7)

where $\big{\langle}...\big{\rangle}$ denotes ensemble average and
$\mathbf{n}$ is the nematic director. Similarly, the pair distribution
functions parallel and perpendicular to the director are given by

$g_{\parallel}(z)=\frac{1}{N_{r}^{2}}\bigg{\langle}\sum\sum\delta(|(\mathbf{r}_{i}-\mathbf{r}_{j})\cdot\mathbf{n}|-z)\bigg{\rangle}$
(8)
$g_{\perp}(r)=\frac{1}{N_{r}^{2}}\bigg{\langle}\sum\sum\delta(|(\mathbf{r}_{i}-\mathbf{r}_{j})\times\mathbf{n}|-r)\bigg{\rangle}$
(9)

where $z$ and $r\equiv\sqrt{x^{2}+y^{2}}$ are the particle center-to-center
distances in the direction of $\mathbf{n}$ and perpendicularly to it,
respectively, and $\delta$ is the Dirac delta function. Additionally, the
probability distribution of the particles at time $t_{0}$, given by the self-
part of the van Hove Function, reads

$G_{s}(z,t)=\frac{1}{N_{r}}\bigg{<}\sum_{j=1}^{N}\delta(z-[{z}_{j}(t+t_{0})-{z}_{j}(t_{0}))]\bigg{>}$
(10)

$G_{s}(r,t)=\frac{1}{N_{r}}\bigg{<}\sum_{j=1}^{N}\delta(r-[{r}_{j}(t+t_{0})-{r}_{j}(t_{0}))]\bigg{>}$
(11)

Finally, the s-ISF gives a measure of the structural relaxation of the system
over time and quantifies the decay of its density fluctuations. In this case,
we also estimate the parallel and perpendicular contributions as follows

$F_{s,z}(t)=\frac{1}{N_{r}}\bigg{<}\sum_{j=1}^{N_{r}}\exp[i\mathbf{q}\cdot[\mathbf{r}_{j}(t+t_{0})-\mathbf{r}_{j}(t_{0}))]\cdot\mathbf{n}]\bigg{>}$
(12)

$F_{s,xy}(t)=\frac{1}{N_{r}}\bigg{<}\sum_{j=1}^{N_{r}}\exp[i\mathbf{q}\cdot[\mathbf{r}_{j}(t+t_{0})-\mathbf{r}_{j}(t_{0}))\times\mathbf{n}]\bigg{>}$
(13)

where the wave vector $\mathbf{q}=\mathbf{q}_{\parallel}+\mathbf{q}_{\perp}$
is defined at the main peaks of the static structure factor, taking values of
$|\mathbf{q}_{\parallel}|\sigma\simeq 1$ and $|\mathbf{q}_{\perp}|\sigma\simeq
6$ for the parallel and perpendicular components with respect to the director,
respectively.

##  III Results

Before presenting and discussing the relevant observations on the dynamics of
SRS rods in N and Sm phases, we first analyse the structural properties that
help us distinguish between positionally and merely orientationally ordered
phases. To this end, we computed the pair-correlation functions $g_{\perp}(r)$
and $g_{\parallel}(z)$, shown in Figs. 1 and 2, respectively. Our results
indicate the existence of a typical fluid-like behavior in the perpendicular
direction of both N and Sm phases, with the first peak occurring at a distance
of approximately one diameter length and the amplitude of oscillations
decaying exponentially to 1 at relatively short distances. Upon increasing
temperature from $T^{\star}=5$ to 20, all peaks tend to flatten and
fluctuations decay at progressively shorter distances. As far as the parallel
pair correlation function, $g_{\parallel}(z)$, is concerned, the pronounced
periodic correlations of Fig. 2 confirm the layered structure typically
observed in Sm LCs. In particular, the position of each peak roughly
corresponds to the location of each smectic layer, where particle density is
maximal. By contrast, the inter-layer spacing is almost completely empty, the
probability of observing particles in between layers being very low.
Increasing temperature has a significant effect on the density distribution in
the Sm phase, but a practically negligible effect in the N phase, which shows
no evidence of positional order at all temperatures studied.

Figure 1: Pair correlation functions in the direction perpendicular to the
nematic director as computed in N (top) and Sm (bottom) phases at different
temperatures. Insets a) and b) in the bottom frame correspond, respectively,
to N ($\rho^{*}=0.12$) and Sm ($\rho^{*}=0.15$) phases at $T^{\star}=10$.
Figure 2: Pair correlation functions in the direction parallel to the nematic
director at $T^{\star}$=5, 10, 15 and 20. Lines and symbols refer to Sm and N
phases, respectively. Plots with larger amplitudes correspond to lower
temperatures. Equilibrium densities are $\rho^{\star}=0.12$ for the N phases
and $\rho^{\star}=0.15$ for the Sm phases. Figure 3: Mean square displacement
in the Sm phase at $T^{\star}=10$ in the direction parallel to the nematic
vector. Dotted, dashed and long-dashed lines refer to independent DMC
simulations at $\delta t_{\rm MC}/\tau=10^{-5}$, $10^{-3}$ and $10^{-2}$,
respectively. Solid lines have been obtained by rescaling the corresponding
dashed lines according to Eq. (5).

Having established the structural differences between the N and Sm phases
studied here, we now investigate how these distinctive morphological
attributes can in turn determine distinct dynamical signatures. To start with,
we first estimated the MSD in the direction parallel and perpendicular to the
nematic director and calculated the long-time diffusion coefficients. To this
end, we ran DMC simulations at different values of the MC time step, between
$\delta t_{\rm MC}/\tau=10^{-5}$ and $10^{-2}$, and then applied Eq. (5) to
rescale the results and recover the unique BD time scale Patti and Cuetos
(2012). The rescaled MSDs collapse into a single master curve as shown in Fig.
3, where we report the parallel MSD in the Sm phase at $T^{\star}=10$. The re-
scaling procedure to obtain the MSD across the whole spectrum of relevant time
scales is the same for all the remaining systems. The resulting master curve
has been obtained by superimposing four separate rescaled MSDs (dashed lines)
calculated over DMC simulations at $\delta t_{\rm MC}/\tau=10^{-5}$,
$10^{-4}$, $10^{-3}$ and $10^{-2}$. The so-calculated MSD exhibits an initial
diffusive regime, mostly determined by the particle geometry, followed by an
intermediate time regime where the presence of neighboring layers, forming a
sort of cage around the particles, slows down diffusion and, finally, a long-
time diffusive regime that fully develops at $t/\tau>1$. Very similar
tendencies have also been detected at $T^{\star}=5$, 8, 12, 15 and 20, with
some differences observed in the extension of the cage effect and the onset of
the long-time diffusive regime. In particular, the effect of temperature on
the MSD is clarified in Fig. 4, where we report parallel and perpendicular
MSDs in N (top frame) and Sm (bottom frame) LCs at $T^{\star}=5$ and 20, that
is the lowest and highest temperatures studied. For the sake of clarity, we do
not show the parallel and perpendicular MSDs at intermediate temperatures,
which exhibit a profile in between those reported in Fig. 4 for N and Sm
phases.

Figure 4: Mean square displacements in N (top frame) and Sm (bottom frame) LCs
at $T^{\star}$=5 (circles) and 20 (triangles). Open and solid symbols refer to
perpendicular and parallel direction, respectively.

We also notice that, in the N phase, at short time scales the relationship
between parallel and perpendicular MSDs is such that $\xi\equiv\Delta
r^{2}_{\perp}/\Delta r^{2}_{\parallel}>1$, while it inverts at intermediate
time scales, when the parallel MSD becomes larger and stays so up to the long-
time diffusive regime. This behavior has also been reported in Brownian
dynamics simulations of SRS rods at $T^{\star}=1.465$ Morillo  _et al._
(2019), a temperature at which the phase behavior of soft spherocylinders can
be mapped on that of hard spherocylinders Cuetos and Martínez-Haya (2015). The
dominant character of the long-time parallel diffusion has also been observed
experimentally in N phases of rod-like viruses Lettinga  _et al._ (2005), but
less clear is whether or not this tendency already exists at short time scales
as observed in simulations. To gain an insight into the effect of temperature
on the ratio between perpendicular and parallel diffusivities, we have
calculated the time when the crossover from $\xi>1$ to $\xi<1$ is observed.
Such an inversion time, referred to as $t_{i}$, does indeed change with
temperature and the results are shown in Fig. 5. We observe that $t_{i}$ is
relatively large at low temperatures and then gradually decreases following an
exponential law of the type $t_{i}/\tau=A+B\exp(-CT^{\star})$, where $A=0.40$,
$B=1.51$ and $C=0.37$ are fitting parameters. It is evident that at large
enough temperatures, $t_{i}$ tends to a constant value, approximately equal to
$0.40\tau$, that will not change significantly up to the I-to-N transition
temperature.

Figure 5: Inversion time between parallel and perpendicular MSD of the N
phase. Symbols are simulation results and the solid line is a fitting
exponential function.

The MSDs are instrumental to calculate the long-time self-diffusion
coefficients and their dependence on temperature. More specifically, parallel
and perpendicular self-diffusion coefficients in N and Sm phases were obtained
from the slope of the corresponding MSDs in the long-time diffusive regime:

$D_{\parallel,\perp}=\lim_{t\to\infty}\frac{1}{2l}\frac{d\bigg{\langle}\sum[\mathbf{r}_{l}(t)-\mathbf{r}_{l}(0)]^{2}\bigg{\rangle}}{dt}$
(14)

where $l=1$ or 2 denotes the dimensionality of particle dynamics associated to
the parallel or perpendicular MSD, respectively. By contrast, the total self-
diffusion coefficient has been calculated as
$D_{Tot}=(D_{\parallel}+2D_{\perp})/3$. The dependence of the self-diffusion
coefficients on temperature in N and Sm phases is presented in the two frames
of Fig. 7. In agreement with previous molecular dynamics simulation of rod-
like molecular liquid crystals Cifelli  _et al._ (2006), we find that the
three sets of long-time self-diffusion coefficients exhibit a dependence on
$T^{\star}$ that is well-described by an Arrhenius-like exponential law, that
reads $D/D_{0}\approx D^{\star}\exp(-E^{\star}/T^{\star})$, with the pre-
exponential factor $D^{\star}$ and activation energy $E^{\star}$ fitting
parameters. We also observe that the dependence of $D_{\perp}$ on temperature
is very similar in both N and Sm phases. At a given temperature, most likely
due to the packing difference between the two LC phases, the numerical value
of $D_{\perp}$ is slightly larger in the N phase than in the Sm phase, but
otherwise $D_{\perp}=D_{\perp}(T^{\star})$ exhibits the same exponential
trend, with very similar fitting parameters, in both frames of Fig. 7. On the
other hand, the diffusion along the director is significantly slower in the Sm
phase ($E^{\star}\approx 42$) than in the N phase ($E^{\star}\approx 11$) by
almost one order of magnitude, most likely due to the layered structure that
hampers the penetration of the rods and thus delays their diffusion along the
nematic director. These findings are in qualitative agreement with former
theoretical, simulation and experimental works that clarified the existence of
free-energy barriers hampering the diffusion of rod-like particles through Sm
layers van Roij  _et al._ (1995); Duijneveldt and Allen (1997); Lettinga and
Grelet (2007); Bier  _et al._ (2008); Patti  _et al._ (2009, 2010). Therefore,
while SRS particles in the N phase preferentially diffuse in the direction of
the nematic director, in the Sm phase they are essentially constrained in a
two-dimensional space, especially at $T^{\star}<10$, where $D_{\parallel}$ is
almost negligible.

Figure 6: Total ($D_{\rm Tot}$), parallel ($D_{\parallel}$) and perpendicular
($D_{\perp}$) self-diffusivities, represented with triangles, squares and
circles, respectively, as a function of the reduced temperature in the N (a)
and Sm (b) phases. Solid lines are exponential fits of the type
$D/D_{0}\approx D^{\star}\exp(-E^{\star}/T^{\star})$, with
$D^{\star}=\\{0.0621,\,0.0196,\,0.0328\\}$ and
$E^{\star}=\\{11.5901,\,12.1927,\,11.6009\\}$ for parallel, perpendicular and
total self-diffusivities, respectively, for the nematic states, and
$D^{\star}=\\{0.0308,\,0.0167,\,0.0129\\}$ and
$E^{\star}=\\{42.243,\,12.1584,\,12.948\\}$ for parallel, perpendicular and
total self-diffusivities, respectively, for the smectic states.

In the light of these considerations, we now turn our attention to the
probability of observing particles that displace significantly shorter or
longer distances than the average particles over the same period of time. The
existence of such particles, here referred to as fast or slow, is corroborated
by the computation of the s-VHFs along the nematic director and
perpendicularly to it, as given, respectively, in Eqs. (10) and (11). To
illustrate this, we show the s-VHFs at $t/\tau=10^{4}$, a time that is
sufficiently long to observe the relevant dynamical features of both N and Sm
phases across the whole spectrum of temperatures studied. In particular, the
parallel s-VHFs of Sm LCs, shown in the top frame of Fig. 7, display
periodically peaked profiles that follow the typical layered structure of this
phase. At increasing temperatures, from $T^{\star}=5$ to 20, these peaks
become less and less pronounced, suggesting a more uniform probability of
finding particles at any distance along the nematic director. Nevertheless, at
relatively low temperature, with the smectic layers well-defined and less
prone to density fluctuations, the profiles unambiguously suggest that
particles preferentially jump from layer to layer and almost no particles are
observed in between. For instance, at $T^{\star}=5$, while most particles are
still in their original layer (primary peak), there exist especially fast
particles that succeeded in diffusing, over the same period of time, to a
contiguous layer (secondary peak). This is also observed at larger
temperatures, but the difference between the height of primary and secondary
peaks gradually softens and eventually disappears at $T^{\star}=20$.

Figure 7: Self-part of the van Hove function at a time $t/\tau=10^{4}$ in N
(bottom) and Sm (top) phases along the nematic director at scaled temperatures
$T^{\star}$= 5 8, 10, 12, 15, 17, and 20.

The bottom frame of Fig. 7 reports similar s-VHFs for the N phase. In this
case, profiles with a maximum at $z=0$ and monotonic decay at relatively long
distances are observed. Most particles are therefore at or very close to their
original position, with few of them fast enough to be displaced substantially
larger distances over the same time window. With increasing temperature, more
and more particles are able to move longer distances and, correspondingly,
less and less are found at their original location. Finally, the perpendicular
s-VHFs shown in Fig. 8 for N (bottom frame) and Sm (top frame) phases reveal
the presence of an interesting variety of particles. At $t/\tau=10^{4}$, most
of them have left their initial position, as indicated by the peak of the
distribution. These particles coexist with others that either remained very
close to their original location or displaced significantly larger distances.
Upon increasing temperature, the probability of observing such slow and fast
particles becomes more and more uniform and would eventually become space-
independent at very large temperatures, at which the system would transform
into an isotropic phase.

Figure 8: Self-part of the van Hove function at a time $t/\tau=10^{4}$ in N
(bottom) and Sm (top) phases in planes perpendicular to the nematic director
for scaled temperatures $T^{\star}$= 5 8, 10, 12, 15, 17, and 20. Symbols as
in figure 7

Temperature also plays a key role in determining the time scale of the
structural relaxation of the system. This has been estimated by computing the
s-ISF in direction parallel and perpendicular to $\bf{n}$. The s-ISFs of N and
Sm phases are respectively shown in Figs. 9 and 10. In both cases, they have
been calculated at the wave vectors corresponding to the peak of the static
structure factor, which are $\mathbf{q}=(0,0,q_{z})$, with $q_{z}\sigma=1$ for
parallel s-ISFs and $\mathbf{q}=(q_{x},q_{y},0)$ with
$\sqrt{(q_{x}^{2}+q_{y}^{2})}\sigma=6$ for perpendicular s-ISFs. Both LC
phases exhibit a relevant difference between parallel and perpendicular
relaxation, with the former taking up to 2 to 3 extra time decades. In all the
cases studied, the decay of the s-ISFs closely follows a stretched-exponential
function of the form $\exp[-(t/t_{r})^{\beta}]$, typically observed in dense
fluids Brambilla  _et al._ (2009), with $t_{r}$ and $\beta$ fitting
parameters. In particular, the exponent $\beta$ is approximately between 0.6
and 0.7 for $F_{s,xy}$, and between 0.8 and 0.9 for $F_{s,z}$, suggesting a
more stretched decay in planes perpendicular to the nematic director than in
the direction parallel to it. These values agree well with those reported in
previous simulation of hard spherocylinders Matena  _et al._ (2010); Belli
_et al._ (2010). While the dependence of $\beta$ on temperature is relatively
mild, the relaxation time $t_{r}$, defined as the time at which $F_{s}=1/e$,
changes significantly with the temperature as can be inferred from Fig. 11,
where $\ln(t_{r}/\tau)$ is plotted as a function of $\ln(T^{\star})$. The so-
calculated relaxation time exhibits a power-law dependence on $T^{\star}$ that
holds in both Sm and N phases.

Figure 9: Temperature dependence of the self-part of the intermediate
scattering function in the N phase along the nematic director (top frame) and
perpendicular to it (bottom frame). Symbols indicate simulation results, while
solid lines are stretched-exponential fits. All s-ISFs have been calculated at
wave vectors $q_{z}\sigma\approx 1$ and
$\sqrt{q_{x}^{2}+q_{y}^{2}}\sigma\approx 6$. Dashed line corresponds to the
$1/e$ value where relaxation of iSF is commonly measured.  Figure 10:
Temperature dependence of the self-part of the intermediate scattering
function in the Sm phase along the nematic director (top frame) and
perpendicular to it (bottom frame). Symbols indicate simulation results, while
solid lines are stretched-exponential fits. All s-ISFs have been calculated at
wave vectors $q_{z}\sigma\approx 1$ and
$\sqrt{q_{x}^{2}+q_{y}^{2}}\sigma\approx 6$. Dashed line corresponds to the
$1/e$ value where relaxation of iSF is commonly measured. Figure 11:
Relaxation times for Sm (top frame) and N (bottom frame) phases. Open and
solid symbols refer respectively to structural relaxation perpendicular to the
nematic director and parallel to it. Solid lines are fitted functions of the
form $t_{r}/\tau=b\cdot{T^{*}}^{a}$ where $a$=-1.36 and $b$=2980

##  IV Conclusions

In summary, we have performed DMC simulations to investigate the dynamics of
soft-repulsive rod-like particles in nematic and smectic LCs. To this end, we
have calculated a spectrum of dynamical properties that helped us characterise
the long-time relaxation decay in the direction of the nematic director and
perpendicularly to it. In particular, the MSD was key to obtain the self-
diffusion coefficients and gain an insight into their dependence on
temperature, the s-VHF unveiled the existence of particles able to cover
significantly longer distances than most of the particles over the same time,
and the s-ISF clarified the timescales of the structural fluctuations’ decay.
Our results suggest a similar fluid-like diffusive behavior in the direction
perpendicular to the nematic director in both N and Sm phases. By contrast,
the crystal-like arrangement of the Sm phase along the nematic director
determines a significantly different dynamics in this direction as compared to
that observed in the N phase. In particular, rods struggle to diffuse from
layer to layer as the very low value of their long-time self-diffusion
coefficients reveals. The dependence of self-diffusivities on temperature
follows an Arrhenius-like exponential law that suggests an especially high
activation energy in the Sm phases at low temperatures. By contrast, at the
same temperatures, the activation energy is up to 4 times smaller in the N
phase. The dependence on temperature of the rods’ dynamics in the N phase is
also inferred from the inversion time, $t_{i}$, which is the time at which the
parallel diffusion becomes faster than the perpendicular diffusion. The
inversion time, which also displays an exponential dependence on temperature,
is larger at low temperatures, but then gradually shorter as the system
approaches the transition to the isotropic phase.

The analysis of the s-VHFs in the direction of the nematic director suggests
the presence of rods that are able to displace significantly longer distances
than the average particle over the same time scale. While such fast particles
are found in both the N and Sm phases, the probability distribution profiles
of their parallel displacements in these phases are not the same. More
specifically, periodically peaked s-VHFs are found in the Sm phases, with the
peaks becoming smoother and smoother at increasing temperature. By contrast,
in the N phase, no peaks are observed, but monotonically decreasing
probabilities that vanish at sufficiently long distances. To some extent, the
parallel s-VHFs of the N phase are very similar to those calculated
perpendicularly to the nematic director, whose profiles in the N and Sm phases
are almost completely indistinguishable. However, the perpendicular s-VHFs
exhibit a maximum at relatively short distances that suggests the existence of
especially slow particles that, even at sufficiently long times, have
displaced less than one rod diameter from their original position. The
simultaneous presence of slow and fast particles contributes to determine the
structural relaxation of the systems, which has been assessed by calculating
the s-ISF. As observed in dense liquids, the s-ISFs exhibit a stretched
exponential decay, which can take up to 3 time decades more along the nematic
director than in the directions perpendicular to it.

The dynamical properties of the soft repulsive model studied here can be
extended to introduce other effects that are relevant for the characterization
of soft matter phases, like polydispersity in viruses Pouget  _et al._ (2011)
and the introduction of electrostatic interactions in models for colloids and
proteins Falcón  _et al._ (2020). Since the Mie-potential has extended the
application of the Lennard-Jones model to describe properties of real
substances Dufal  _et al._ (2015), we can expect that in similar way shifted
Kihara systems of variable range, modifying the exponents 12-6, could be also
relevant in order to describe liquid crystalline phases.

##  V Acknowledgments

D.C. acknowledges support from CONACYT for funding her PhD scholarship. A.P.
acknowledges financial support from the Leverhulme Trust Research Project
Grant No. RPG-2018-415 and the Newton Mobility Grant NMG\R2\170137 awarded by
The Royal Society to fund his visit at the University of Guanajuato in León.

## References

  * Sacanna  _et al._ (2013a) S. Sacanna, D. J. Pine, and G.-R. Yi, Soft Matter 9, 8096 (2013a).
  * Huber  _et al._ (2018) L. Huber, R. Suzuki, T. Krüger, E. Frey, and A. R. Bausch, Science 361, 255 (2018).
  * Fraden  _et al._ (1989) S. Fraden, G. Maret, D. L. D. Caspar, and R. B. Meyer, Phys. Rev. Lett. 63, 2068 (1989).
  * Dogic and Fraden (2006) Z. Dogic and S. Fraden, Current Opinion in Colloid & Interface Science 11, 47 (2006).
  * Pouget  _et al._ (2011) E. Pouget, E. Grelet, and M. P. Lettinga, Phys. Rev. E 84, 041704 (2011).
  * Sacanna  _et al._ (2013b) S. Sacanna, M. Korpics, K. Rodriguez, L. Colón-Meléndez, S.-H. Kim, D. J. Pine, and G.-R. Yi, Nat. Commun. 4, 1688 (2013b).
  * Murphy  _et al._ (2020) R. P. Murphy, H. W. Hatch, A. Mahynski, Nathan, V. K. Shen, and N. J. Wagner, Soft Matter 16, 1279 (2020).
  * Frenkel (1988) D. Frenkel, J. Phys. Chem. 92, 3280 (1988).
  * McGrother  _et al._ (1996) S. C. McGrother, D. C. Williamson, and G. Jackson, j. Chem. Phys. 104, 6755 (1996).
  * McGrother  _et al._ (1998) S. C. McGrother, A. Gil-Villegas, and G. Jackson, Mol. Phys. 95, 657 (1998).
  * van Duijneveldt  _et al._ (2000) J. S. van Duijneveldt, A. Gil-Villegas, G. Jackson, and M. P. Allen, J. Chem. Phys. 112, 9092 (2000).
  * Allen (2019) M. P. Allen, Mol. Phys. 117, 2391 (2019).
  * Guevara-Rodríguez and Medina-Noyola (2003) F. J. Guevara-Rodríguez and M. Medina-Noyola, Phys. Rev. E 68, 011405 (2003).
  * Ramírez-González  _et al._ (2011) P. Ramírez-González, L. López-Flores, H. Acuña Campa, and M. Medina-Noyola, Phys. Rev. Lett. 107, 155701 (2011).
  * Jover  _et al._ (2012) J. Jover, A. J. Haslam, A. Galindo, G. Jackson, and E. A. Müller, J. Chem. Phys. 137, 144505 (2012).
  * Zerón  _et al._ (2018) I. M. Zerón, C. Vega, and A. L. Benavides, Mol. Phys. 116, 3355 (2018).
  * Báez  _et al._ (2018) C. A. Báez, A. Torres-Carbajal, R. Castañeda Priego, A. Villada-Balbuena, J. M. Méndez-Alcaraz, and S. Herrera-Velarde, J. Chem. Phys. 149, 164907 (2018).
  * Davies  _et al._ (1998) A. L. Davies, A. Gil-Villegas, and G. Jackson, Int. J. Thermophys. 19, 675 (1998).
  * Lafitte  _et al._ (2013) T. Lafitte, A. Apostolakou, C. Avendaño, A. Galindo, C. S. Adjiman, E. A. Muller, and G. Jackson, J. Chem. Phys. 139, 154504 (2013).
  * Dufal  _et al._ (2015) S. Dufal, T. Lafitte, A. Galindo, G. Jackson, and A. J. Haslam, AIChe J 61, 2891 (2015).
  * Metropolis  _et al._ (1953) A. W. Metropolis, N.and Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953).
  * Sanz and Marenduzzo (2010) E. Sanz and D. Marenduzzo, The Journal of Chemical Physics 132, 194102 (2010), https://doi.org/10.1063/1.3414827 .
  * Romano  _et al._ (2011) F. Romano, C. De Michele, D. Marenduzzo, and E. Sanz, The Journal of Chemical Physics 135, 124106 (2011), https://doi.org/10.1063/1.3629452 .
  * Patti and Cuetos (2012) A. Patti and A. Cuetos, Phys. Rev. E 86, 011403 (2012).
  * Cuetos and Patti (2015) A. Cuetos and A. Patti, Phys. Rev. E 92, 022302 (2015).
  * Corbett  _et al._ (2018) D. Corbett, A. Cuetos, M. Dennison, and A. Patti, Phys. Chem. Chem. Phys. 20, 15118 (2018).
  * Chiappini  _et al._ (2020) M. Chiappini, A. Patti, and M. Dijkstra, Phys. Rev. E 102, 040601 (2020).
  * García Daza  _et al._ (2020) F. A. García Daza, A. Cuetos, and A. Patti, Phys. Rev. E 102, 013302 (2020).
  * García Daza  _et al._ (2022) F. A. García Daza, A. M. Puertas, A. Cuetos, and A. Patti, Journal of Colloid and Interface Science 605, 182 (2022).
  * Cuetos  _et al._ (2005) A. Cuetos, B. Martínez-Haya, S. Lago, and L. F. Rull, J. Phys. Chem. B 109, 13729 (2005).
  * Cuetos and Martínez-Haya (2015) A. Cuetos and B. Martínez-Haya, Mol. Phys. 113, 1137 (2015), https://doi.org/10.1080/00268976.2014.996191 .
  * Cuetos  _et al._ (2002) A. Cuetos, B. Martínez-Haya, L. F. Rull, and S. Lago, J. Chem. Phys. 117, 2934 (2002).
  * Vega and Lago (1994) C. Vega and S. Lago, Computers & Chemistry 18, 55 (1994).
  * Bonet Avalos  _et al._ (1994) J. Bonet Avalos, J. M. Rubí, D. Bedeaux, and G. van der Zwan, Physica A 211, 193 (1994).
  * Morillo  _et al._ (2019) N. Morillo, A. Patti, and A. Cuetos, The Journal of Chemical Physics 150, 204905 (2019), https://doi.org/10.1063/1.5090975 .
  * Lettinga  _et al._ (2005) M. P. Lettinga, E. Barry, and Z. Dogic, EPL (Europhysics Letters) 71, 692 (2005).
  * Cifelli  _et al._ (2006) M. Cifelli, G. Cinacchi, and L. De Gaetani, J. Chem. Phys. 125, 164912 (2006), https://doi.org/10.1063/1.2359428 .
  * van Roij  _et al._ (1995) R. van Roij, P. Bolhuis, B. Mulder, and D. Frenkel, Phys. Rev. E 52, R1277 (1995).
  * Duijneveldt and Allen (1997) J. S. Duijneveldt and M. P. Allen, Molecular Physics 90, 243 (1997), https://doi.org/10.1080/002689797172723 .
  * Lettinga and Grelet (2007) M. P. Lettinga and E. Grelet, Phys. Rev. Lett. 99, 197802 (2007).
  * Bier  _et al._ (2008) M. Bier, R. van Roij, M. Dijkstra, and P. van der Schoot, Phys. Rev. Lett. 101, 215901 (2008).
  * Patti  _et al._ (2009) A. Patti, D. El Masri, R. van Roij, and M. Dijkstra, Phys. Rev. Lett. 103, 248304 (2009).
  * Patti  _et al._ (2010) A. Patti, D. El Masri, R. van Roij, and M. Dijkstra, The Journal of Chemical Physics 132, 224907 (2010), https://doi.org/10.1063/1.3432864 .
  * Brambilla  _et al._ (2009) G. Brambilla, D. El Masri, M. Pierno, L. Berthier, L. Cipelletti, G. Petekidis, and A. B. Schofield, Phys. Rev. Lett. 102, 085703 (2009).
  * Matena  _et al._ (2010) R. Matena, M. Dijkstra, and A. Patti, Phys. Rev. E 81, 021704 (2010).
  * Belli  _et al._ (2010) S. Belli, A. Patti, R. van Roij, and M. Dijkstra, The Journal of Chemical Physics 133, 154514 (2010), https://doi.org/10.1063/1.3505150 .
  * Falcón  _et al._ (2020) J. M. Falcón, C. Contreras-Aburto, M. Lara Peña, M. Heinen, C. Avendaño, A. Gil-Villegas, and R. Castañeda Priego, J. Chem. Phys. 153, 234901 (2020).