# Set Systems with Covering Properties and Low VC-Dimension

George Peterzil  Einstein Institute of Mathematics
The Hebrew University of Jerusalem
91904 Jerusalem
Israel<EMAIL_ADDRESS>and  Johanna Steinmeyer
<EMAIL_ADDRESS>
(Date: January 20, 2024)

###### Abstract.

Given natural numbers $k\leq s\leq n$, we ask: what is the minimal VC-
dimension of a family $\mathcal{F}$ of $s$-subsets of $[n]$ that covers all
$k$-subsets of $[n]$? We first show that for sufficiently large $n$ this
number is always $k$, and construct families which give a lower bound for the
actual growth of this stabilization point.

Consider a family of sets $\mathcal{F}$. We say that a set $A$ is shattered by
$\mathcal{F}$ if

$2^{A}=\left\\{A\cap S\mid S\in\mathcal{F}\right\\}.$

Define the VC-dimension of $\mathcal{F}$ to be the size of the largest finite
set shattered by $\mathcal{F}$, or $\infty$ if there is no such maximum.

Originiating from machine learning in [VC15], the notion of VC-dimension has
since been further studied also by combinatorialists (e.g. [FKKP22] and
[AMY15], for an overview see [Dud99]) as well as logicians (for a survey see
[Sim15]).

Consider the family of initial segments of the natural numbers. It
demonstrates that one can cover all finite subsets of $\mathbb{N}$ with a
family of finite sets of VC-dimension $1$. Bays, Ben-Neria, Kaplan and Simon
prove the following theorem, which can be seen as an uncountable analogue of
the above fact:

######  Theorem 0.1 ([BBKS23], Proposition 3.3 and Theorem 3.8).

Every family of finite sets which covers every finite subset of the first
uncountable cardinal has VC-dimension at least $2$, and there is such a family
with VC-dimension exactly $2$.

This result is motivated by model theory, namely by the existence of so-called
honest definitions for dependent local types (see [CS15], Definition 1.6).
This means that sets which are definable with a parameter _external_ to a
first-order structure by a formula of finite VC-dimension have their finite
subsets covered by a definable family (with an _internal_ parameter). It is
therefore natural to ask for effective results: suppose $k<s<n$ are natural
numbers. What is the minimal VC-dimension $D(k,s,n)$ of a family of subsets of
$[n]$, each of size $s$, which covers every $k$-element subset of $[n]$?

Finitary questions regarding VC-classes appear across the literature. A
celebrated one of these is the following, proved independently by Sauer and
Shelah.

######  Theorem 0.2 ([Sau72],[She72]).

Given $\mathcal{F}\subseteq 2^{[n]}$, if
$|\mathcal{F}|>\sum_{i=0}^{k-1}\binom{n}{i}$ then the VC-dimension of
$\mathcal{F}$ is at least $k$.

We can see our question above as an analogue of the Sauer-Shelah lemma,
replacing the assumption on the size of the family by the covering
assumptions. The following is our main result.

######  Theorem 0.3 (3.7).

For all $k\leq s$ and $n\in\mathbb{N}$ with $k^{2}\binom{s}{k}+k\leq n$ we
have $D(k,s,n)=k$.

In addition to this result, we construct $k$-covering families with VC-
dimension smaller than $k$ on smaller sets to demonstrate the requirement for
large $n$ in the above theorem. These results utilize particular examples of
families as well as canonical constructions of new families with prescribed
properties from existing ones.

Answering analogous questions for infinite cardinals is part of work in
progress of Omer Ben-Neria, Itay Kaplan and the first author.

## Acknowledgements

This work is part of G.P.’s Master’s thesis, done in the Hebrew University of
Jerusalem. He would like to thank his advisors Omer Ben-Neria and Itay Kaplan
for introducing him to this flavor of questions. G.P. is supported by the
Israel Science Foundation (Grant 1832/19). J.S. is supported by the Schaerf
fund of the Einstein Institute of Mathematics.

##  1\. Definitions and Examples

Given a natural number $n$, write $[n]=\\{1,...,n\\}$. For a set $A$ and
$s\in\mathbb{N}$, we write $\binom{A}{s}=\big{\\{}A_{0}\subseteq
A:|A_{0}|=s\big{\\}}$.

######  Definition 1.1.

Let $\mathcal{F}$ be a family of sets.

  * •

Say that a set $A$ is shattered by $\mathcal{F}$ if $2^{A}=\left\\{A\cap S\mid
S\in\mathcal{F}\right\\}$.

  * •

Define the VC-dimension of $\mathcal{F}$ to be the size of the largest finite
set shattered by $\mathcal{F}$, or $\infty$ if there is no such maximum.

######  Examples 1.2.

  * •

Given natural numbers $s\leq n$, the family $\binom{[n]}{s}$ has VC-dimension
$\min\\{s,n-s\\}$; the set $[\min\\{s,n-s\\}]\subseteq[n]$ is shattered, as
any $A\subseteq[\min\\{s,n-s\\}]$ of size $k$ has
$A=[s]\cap(A\cup\\{s+1,\dots,2s-k\\})$, and any larger set does not have the
empty set as the intersection of it with a member of the family by the
pigeonhole principle.

  * •

The family $\mathcal{F}$ of proper initial segments in a linearly ordered set
$X$ with at least two distinct elements $x<y$ has VC-dimension $1$. Indeed, it
is at least $1$ since there is $I\in\mathcal{F}$ containing $y$ and there is
$J\in\mathcal{F}$ which contains $x$ and not $y$, hence $\\{y\\}$ is
shattered. However, the set $\\{x,y\\}$ is not shattered – every
$I\in\mathcal{F}$ which contains $y$ contains $x$ as well.

######  Definition 1.3.

  * •

Let $X$ be a set, $\mathcal{F}\subseteq 2^{X}$, $k\in\mathbb{N}$ with
$k\leq|X|$. Say that $\mathcal{F}$ has the k-covering property if every
$A\in\binom{X}{k}$ is contained in a member of $\mathcal{F}$.

  * •

Given $k\leq s\leq n$, write $D(k,s,n)$ for the smallest VC-dimension of a
family $\mathcal{F}\subseteq\binom{[n]}{s}$ which has the $k$-covering
property.

Below are a few simple instances of calculations of this function.

######  Examples 1.4.

  * •

For all $k\leq n$ we have $D(k,n,n)=0$, as witnessed by the trivial family
$\\{[n]\\}$.

  * •

For all $k\leq s$ we have $D(k,s,n)\leq\min\\{s,n-s\\}$ by • ‣ 1.2.

  * •

Since the only $\mathcal{F}\subseteq\binom{[n]}{k}$ with the $k$-covering
property is $\binom{[n]}{k}$ itself, we have $D(k,k,n)=\min\\{k,n-k\\}$ by • ‣
1.2.

##  2\. General Bounds and Constructions

Very few families with covering properties can have VC-dimension $1$. The
following is a simple instance of this.

######  Proposition 2.1.

For $2\leq k\leq s$ such that $s<2k$ and $s+k<n$ we have $1<D(k,s,n)$.

###### Proof.

Take $S_{1}\in\mathcal{F}$. Using the $k$-covering property, pick $k$ points
in $[n]\setminus S_{1}$ and some $S_{2}\in\mathcal{F}$ which covers them, so
$|S_{1}\cap S_{2}|\leq s-k<k$. Take some $T\in\mathcal{F}$ which contains
$S_{1}\cap S_{2}$ and some point outside of $S_{1}\cup S_{2}$. Since all
members of $\mathcal{F}$ have size $s$, $T$ does not contain $S_{1}\setminus
S_{2}$ nor $S_{1}\setminus S_{2}$, hence there is some $p_{1}\in
S_{1}\setminus\big{(}S_{2}\cup T\big{)}$ and some $p_{2}\in
S_{2}\setminus\big{(}S_{1}\cup T\big{)}$. Then $\big{\\{}p_{1},p_{2}\big{\\}}$
is shattered by $S_{1},S_{2},T$ and any set which contains both $p_{1}$ and
$p_{2}$. ∎

######  Theorem 2.2.

Fix $k\leq s$. Then for sufficiently large $n\in\mathbb{N}$ we have $k\leq
D(k,s,n)$. In particular, one can take any $n\geq k^{2}\binom{s}{k}+k$.

###### Proof.

A short calculation shows that this is equivalent to
$k\binom{n}{k-1}<\frac{\binom{n}{k}}{\binom{s}{k}}$. Now since $2k\leq n$,
this gives us
$\sum_{i=0}^{k-1}\binom{n}{i}<\frac{\binom{n}{k}}{\binom{s}{k}}$. Let
$\mathcal{F}$ be a $k$-covering family of subsets of $[n]$ of size $s$, then
any $S\in\mathcal{F}$ can cover at most $\binom{s}{k}$ distinct subsets of
$[n]$ of size $k$. Since there are $\binom{n}{k}$ such subsets, the
cardinality of $\mathcal{F}$ is at least $\frac{\binom{n}{k}}{\binom{s}{k}}$,
hence by 0.2 we have $\operatorname{VC}(\mathcal{F})\geq k$. ∎

Our original approach to this bound was based on proving a finitary version of
Proposition 3.8 of [BBKS23], however it gave significantly worse estimates for
$n$ than the proof above.

The following two Propositions demonstrate the flexibility of pushing up $s$
and $n$, while keeping the VC-dimension and the $k$-covering property.

######  Proposition 2.3.

Suppose $\mathcal{F}\subseteq\binom{[n]}{s}$ has VC-dimension $m$ and the
$k$-covering property.

  1. (2.3.1)

There is a family $\mathcal{F}_{*}\subseteq\binom{[n+1]}{s+1}$ of VC-dimension
$m$ with the $k$-covering property (in fact, this family shatters the same
sets as $\mathcal{F})$. In particular, for every $k\leq s\leq n$ we have
$D(k,s+1,n+1)\leq D(k,s,n)$.

  2. (2.3.2)

For every natural number $\ell$, there is a family
$\mathcal{F}\times\ell\subseteq\binom{[n\cdot\ell]}{s\cdot\ell}$ of VC-
dimension $m$ with the $k$-covering property. In particular, for every $k\leq
s\leq n$ we have $D(k,s\cdot\ell,n\cdot\ell)\leq D(k,s,n)$.

###### Proof.

  1. (1)

Define $\mathcal{F}_{*}$ to be the family $\\{A\cup\\{n+1\\}\mid
A\in\mathcal{F}\\}\subseteq\binom{[n+1]}{s+1}$. We will claim that it has the
$k$-covering property. Indeed, take $A\subseteq[n+1]$ of size $k$. If it does
not contain $n+1$, take some $B\in\mathcal{F}$ containing it and cover it with
$B\cup\\{n+1\\}$. Otherwise, take some $x\in[n]\setminus A$ and some
$B\in\mathcal{F}$ containing $A\cup\\{x\\}\setminus\\{n+1\\}$, then
$A\cup\\{n+1\\}$ contains $A$. We will now claim that
$\operatorname{VC}\left(\mathcal{F}_{*}\right)\leq\operatorname{VC}(\mathcal{F})$.
Suppose $A\subseteq[n+1]$ is shattered by $\mathcal{F}_{*}$. If
$A\subseteq[n]$ then necessarily $A$ is shattered by $\mathcal{F}$ since for
every $S\in\mathcal{F}$ we have $S\cap A=S\cup\\{n+1\\}\cap A$. Otherwise we
have $n+1\in A$, which implies that every $S\in\mathcal{F}_{*}$ has $n+1\in
S\cap A$, a contradiction to $A$ being shattered by $\mathcal{F}_{*}$. We will
now show that every set shattered by $\mathcal{F}$ is shattered by
$\mathcal{F}_{*}$ as well. Indeed, suppose $A\subseteq[n+1]$ is shattered by
$\mathcal{F}$, so necessarily $A\subseteq[n]$, hence for every
$S\in\mathcal{F}$ we have $S\cup\\{n+1\\}\cap A=S\cap A$, so $A$ is shattered
by $\mathcal{F}_{*}$ as well.

  2. (2)

Define $\mathcal{F}\times\ell\subseteq\binom{[n]\times[\ell]}{s\cdot\ell}$ as
follows:

$\displaystyle\mathcal{F}\times\ell=\\{S\times[\ell]\mid S\in\mathcal{F}\\}$

We will first claim that $\mathcal{F}\times\ell$ has the $k$-covering
property. Fix $p_{i}=\left(v_{i},x\right)\in[n]\times[\ell]$ for $1\leq i\leq
n$. Take some $S\in\mathcal{F}$ which covers $\big{\\{}v_{i}\mid 1\leq i\leq
k\big{\\}}$, so $S\times[\ell]$ covers $\big{\\{}p_{1},\dots,p_{k}\big{\\}}$.
We will now claim that $\operatorname{VC}(\mathcal{F}\times\ell)\leq
m=\operatorname{VC}(\mathcal{F})$. Fix distinct points
$p_{i}=\left(v_{i},x_{i}\right)\in[n]\times[\ell]$ for $1\leq i\leq m+1$. If
there are $1\leq i_{1}<i_{2}\leq m+1$ with $v_{i_{1}}=v_{i_{2}}$, then for
every $S\in\mathcal{F}$ with $p_{i_{1}}\in S\times[\ell]$ we have
$p_{i_{2}}\in S\times[\ell]$, hence $\mathcal{F}\times\ell$ does not shatter
$\big{\\{}p_{i}\mid 1\leq i\leq m+1\big{\\}}$, so assume all $v_{i}$’s are
distinct. By the definition of $v$, there is some
$A\subseteq\big{\\{}v_{i}\mid 1\leq i\leq m+1\big{\\}}$ with $A\neq
S\cap\big{\\{}v_{i}\mid 1\leq i\leq m+1\big{\\}}$ for all $S\in\mathcal{F}$,
so $\big{\\{}p_{i}\mid v_{i}\in
A\big{\\}}\neq(S\times[\ell])\cap\big{\\{}p_{i}\mid 1\leq i\leq m+1\\}$ for
all $S\in\mathcal{F}$. To get equality in VC-dimension, note that for every
$A\subseteq[n]$ and $S\in\mathcal{F}$, we have $(A\cap
S)\times\\{1\\}=A\times\\{1\\}\cap S\times[\ell]$.

∎

##  3\. Specific Families

The following is a family of examples that shows that for small enough $n$, we
can expect low VC-dimension.

######  Proposition 3.1.

For all $k,m\in\mathbb{N}$ we have the following bound:

$D(k,k^{m},(k+1)^{m})\leq m\log_{2}(k+1)$

###### Proof.

Consider the abelian group $A=\big{(}\mathbb{Z}/(k+1)\big{)}^{m}$. Define the
family

$\mathcal{F}=\left\\{\prod_{i=1}^{m}X_{i}\mid
X_{i}\in\binom{\mathbb{Z}/(k+1)}{k}\right\\}\subseteq\binom{A}{k^{m}},$

where the members $\mathcal{F}$ are $m$-hypercubes of width $k$ in $A$. For
every $1\leq i\leq m$ and $p\in A$ define $\pi_{i}(p)$ to be the $i$’th
coordinate of $p\leavevmode\nobreak\ $. We will claim that $\mathcal{F}$ has
the $k$-covering property. Indeed, given points $p_{1},\dots,p_{k}\in A$, for
every $1\leq i\leq m$, there is at least one $1\leq t_{i}\leq k+1$ such that
$\pi_{i}(p_{j})\neq t_{i}$ for all $1\leq j\leq k$. This means that
$p_{1},...,p_{k}$ are all contained in the following hypercube:

$\prod_{i=1}^{m}\big{(}\big{(}\mathbb{Z}/(k+1)\big{)}\setminus\big{\\{}t_{i}\big{\\}}\big{)}$

This gives us the $k$-covering property. As for the VC-dimension, note that to
shatter a set of size $t$, the family must be of cardinality at least $2^{t}$.
Since $\mathcal{F}$ contains $|A|=(k+1)^{m}$ many elements, the result
follows. ∎

Figure 1. A member of $\mathcal{F}$ for $k=5$.

The following demonstrates that we can achieve the bound in 2.2 for all $n$.

######  Definition 3.2.

Fix $s\leq n$ and a family of sets $\mathcal{F}\subseteq\binom{[n]}{s}$. Say
that $\mathcal{F}$ has the _unique face property_ if for every
$S\in\mathcal{F}$ there is $K\subsetneq S$ such that the only member of
$\mathcal{F}$ containing $K$ is $S$.

######  Proposition 3.3.

If $\mathcal{F}\subseteq\binom{[n]}{s}$ has the unique face property, then it
has VC-dimension less than $s$.

###### Proof.

Fix $A\in\binom{[n]}{s}$. If $A\notin\mathcal{F}$ then $A$ itself is not of
the form $A\cap S$ for any $S\in\mathcal{F}$. Otherwise, by the unique face
property we have some $K\subseteq A$ such that the only member of
$\mathcal{F}$ containing it is $A$. In particular, $K$ cannot be separated
from $A\setminus K$ by $\mathcal{F}$. ∎

######  Construction 3.4.

Given $m\in\mathbb{N}$ and $k=1$, consider the following family:

$\mathcal{F}_{1}=\\{\\{2t-1,2t\\}\mid
t\leq\frac{m}{2}\\}\cup\\{\\{m-1,m\\}\\}\subseteq\binom{[m]}{2}$

Assuming we have defined $\mathcal{F}_{k}\subseteq\binom{[m+k-1]}{k}$, define
the family $\mathcal{F}_{k+1}\subseteq\binom{[m+k]}{k+1}$ as follows:

$\displaystyle\mathcal{F}_{k+1}=\bigcup_{i=1}^{m+k}\\{S\cup\\{i\\}\mid
S\in\mathcal{F}\cap 2^{[i-1]}\\}$

Figure 2. $\mathcal{F}_{2}$ when starting with $m=4$.

######  Proposition 3.5.

For all $m,k\in\mathbb{N}$, the family
$\mathcal{F}_{k}\subseteq\binom{[n]}{k+1}$ (for $n=m+k-1$) satisfies the
following properties:

  1. (1)

It has the $k$-covering property.

  2. (2)

It has the unique face property.

  3. (3)

For every ordered $\\{t_{1},\dots,t_{k+1}\\}\in\mathcal{F}_{k}$ and
$\hat{t}\in[n]$ with $t_{k}<\hat{t}<t_{k+1}$, we have
$\\{t_{1},\dots,t_{k},\hat{t}\\}\in\mathcal{F}_{k}$ as well.

  4. (4)

If $2k<n$ then $\big{\\{}n-k+1,\dots,n\big{\\}}$ is shattered by
$\mathcal{F}_{k}$.

###### Proof.

We will prove the result by induction on $k$. For $k=1$ the first two are
trivial, the third is void, and the fourth follows from
$\\{1,2\\}\in\mathcal{F}_{1}$ and $\\{n-1,n\\}\in\mathcal{F}_{1}$.

Assume now that these three properties hold for $\mathcal{F}_{k}$. Pick an
ordered $K=\big{\\{}t_{1},\dots,t_{k+1}\big{\\}}\subseteq[n+1]^{k+1}$ and let
$S\in\mathcal{F}_{k}$ cover $K\setminus\big{\\{}t_{k+1}\big{\\}}$. If $S$
covers $K$ then so does $S\cup\\{n+1\\}\in\mathcal{F}_{k+1}$. If $S$ does not
cover $K$ and $\max{S}<t_{k+1}$ then $S\cup\big{\\{}t_{k+1}\big{\\}}$.
Otherwise, there is some $r>0$ with $t_{k+1}+r=\max{S}$ (since $|S|=k+1)$.
Applying 3 to $S$ and $\hat{t}=t_{k+1}$ gives us
$S^{\prime}\in\mathcal{F}_{k}$ with $K\subseteq S$, reducing to the first
case.

For the unique face property, fix $S^{\prime}\in\mathcal{F}_{k+1}$, so some
$S\in\mathcal{F}_{k}$ has $S^{\prime}=S\cup\\{\max S^{\prime}\\}$. For $K_{0}$
a unique face of $S$, we will claim that $K=K_{0}\cup\\{\max S^{\prime}\\}$ is
a unique face of $S$. Suppose $T^{\prime}\in\mathcal{F}_{k+1}$ has $K\subseteq
T^{\prime}$, so some $T\in\mathcal{F}_{k}$ has $T^{\prime}=T\cup\\{\max
T^{\prime}\\}$. If $K\subseteq T$ then $K_{0}\subseteq T$, giving us $S=T$ and
$\max T^{\prime}=\max S^{\prime}$, hence $S=T$ and $S^{\prime}=T^{\prime}$.
Otherwise we have $\max T^{\prime}\in K$, hence $\max T^{\prime}=\max
S^{\prime}$ and so $S^{\prime}=T^{\prime}$ nonetheless.

Suppose $S^{\prime}=\big{\\{}s_{1},...,s_{k+2}\big{\\}}\in\mathcal{F}_{k+1}$
is ordered, so
$S=S^{\prime}\setminus\big{\\{}s_{k+2}\big{\\}}\in\mathcal{F}_{k}$. Fix
$s_{k+1}<\hat{t}<s_{k+2}$. By definition of $\mathcal{F}_{k+1}$, since $\max
S=s_{k+1}<\hat{t}$ we get $S\cup\\{\hat{t}\\}$ as well.

Assume that $2k+1<n$ and the claim holds up to $k$. By the induction
hypothesis, $\\{n-k+1,\dots,n\\}$ is shattered by $\mathcal{F}_{k}$. To
separate the empty set, note that $\\{1,\dots,k+1\\}\in\mathcal{F}_{k+1}$ and
$k+1<n-k+1$. For $A\subseteq\\{n-k+1,\dots,n+1\\}$ nonempty, let $a\in A$ be
maximal. By the induction hypothesis, there is some $S\in\mathcal{F}_{k}$ such
that $S\cap\\{n-k+1,\dots,n\\}=A\setminus\\{a\\}$. In particular we have
$S\subseteq 2^{[a-1]}$, so by definition we have
$S\cup\\{a\\}\in\mathcal{F}_{k+1}$, hence
$A=(S\cup\\{a\\})\cap\\{n-k+1,\dots,n+1\\}$. ∎

######  Corollary 3.6.

For all $k\leq s\leq n$ we have $D(k,s,n)\leq k$.

###### Proof.

For all $k\leq n$ have $D(k,k,n)\leq k$ by • ‣ 1.4. For $s=k+1$ we get the
result by 3.5 and 3.3. For general $s$, apply 2.3.1 $s-k-1$ many times to the
family $\mathcal{F}_{k}\subseteq\binom{[n-(s-k-1)]}{k+1}$. ∎

Joining 3.6 with 2.2 gives us our main theorem:

######  Theorem 3.7.

For all $k\leq s$ and $n\geq k^{2}\binom{k}{s}+k$ we have $D(k,s,n)=k$.

We end on a few natural questions arising from our analysis of the function
$D(k,s,n)$.

######  Question 3.8 (Stabilization estimates).

Theorem 3.7 gives us that for fixed $k\leq s$, the function $D(k,s,-)$
stabilizes eventually on the value $k$, so given $k,s$, we can define
$\operatorname{Stab}(k,s)$ to be the first $N$ such that for every $n\geq N$
we have $D(k,s,n)=k$. What is the actual growth rate of
$\operatorname{Stab}(k,s)$?

Another natural question is for the values that the function $D(k,s,-)$
attains before stabilizing.

######  Question 3.9 (Surjectivity).

Does the function $D(k,s,-)$ attain every natural number below $k$?

If we had considered family of sets of size _at most s_ rather than _exactly
s_ , the following question would have an immediate positive answer. However,
in this case it seems quite tricky, and might have a negative answer due to
divisibility issues.

######  Question 3.10 (Monotonicity).

When fixing two of the parameters and varying the third of $k,s,n$, is the
corresponding function monotone?

## References

  * [AMY15] Noga Alon, Shay Moran and Amir Yehudayoff  “Sign rank versus VC dimension”  In _Annual Conference Computational Learning Theory_ , 2015
  * [BBKS23] Martin Bays, Omer Ben-Neria, Itay Kaplan and Pierre Simon  “On large definable subsets in NIP”  In _Journal of the Institute of Mathematics of Jussieu_ Cambridge University Press, 2023, pp. 1–15  DOI: 10.1017/S1474748023000464
  * [CS15] Artem Chernikov and Pierre Simon  “Externally definable sets and dependent pairs II”  In _Transactions of the American Mathematical Society_ 367.7 American Mathematical Society, 2015, pp. 5217–5235  URL: http://www.jstor.org/stable/24513386
  * [Dud99] R.. Dudley  “Vapnik-Červonenkis Combinatorics”  In _Uniform Central Limit Theorems_ , Cambridge Studies in Advanced Mathematics  Cambridge University Press, 1999, pp. 134–169
  * [FKKP22] Nóra Frankl, Sergei Kiselev, Andrey Kupavskii and Balázs Patkós  “VC-saturated set systems”  In _European Journal of Combinatorics_ 104, 2022  DOI: https://doi.org/10.1016/j.ejc.2022.103528
  * [Sau72] Norbert Sauer  “On the density of families of sets”  In _Journal of Combinatorial Theory_ 13, A, 1972, pp. 145–147  DOI: 10.1016/0097-3165(72)90019-2
  * [She72] Saharon Shelah  “A combinatorial problem; stability and order for models and theories in infinitary languages”  In _Pacific Journal of Mathematics_ 41, 1972, pp. 247–261
  * [Sim15] Pierre Simon  “A Guide to NIP Theories”, Lecture Notes in Logic  Cambridge University Press, 2015
  * [VC15] Vladimir Vapnik and Alexey Chervonenkis  “On the Uniform Convergence of Relative Frequencies of Events to Their Probabilities”  In _Measures of Complexity: Festschrift for Alexey Chervonenkis_ Cham: Springer International Publishing, 2015, pp. 11–30  DOI: 10.1007/978-3-319-21852-6_3


# Harissa: Stochastic Simulation and Inference of Gene Regulatory Networks
Based on Transcriptional Bursting

Ulysse Herbach

(Université de Lorraine, CNRS, Inria, IECL, F-54000 Nancy, France
<EMAIL_ADDRESS>
Abstract. Gene regulatory networks, as a powerful abstraction for describing
complex biological interactions between genes through their expression
products within a cell, are often regarded as virtually deterministic
dynamical systems. However, this view is now being challenged by the
fundamentally stochastic, ‘bursty’ nature of gene expression revealed at the
single cell level. We present a Python package called Harissa which is
dedicated to simulation and inference of such networks, based upon an
underlying stochastic dynamical model driven by the transcriptional bursting
phenomenon. As part of this tool, network inference can be interpreted as a
calibration procedure for a mechanistic model: once calibrated, the model is
able to capture the typical variability of single-cell data without requiring
ad hoc external noise, unlike ordinary or even stochastic differential
equations frequently used in this context. Therefore, Harissa can be used both
as an inference tool, to reconstruct biologically relevant networks from time-
course scRNA-seq data, and as a simulation tool, to generate quantitative gene
expression profiles in a non-trivial way through gene interactions.

Keywords: Gene expression · Regulatory networks · Single-cell data

##  1 Introduction

Inferring graphs of interactions between genes has become a standard task for
high-dimensional statistics, while mechanistic models describing gene
expression at the molecular level have come into their own with the advent of
single-cell data. Linking these two approaches seems crucial today, but the
dialogue is far from obvious: statistical models often suffer from a lack of
biological interpretability, and mechanistic models are known to be difficult
to calibrate from real data.

Here, we present a Python package for both network simulation and inference
from single-cell gene expression data (typically scRNA-seq), called Harissa
(‘HARtree approximation for Inference along with a Stochastic Simulation
Algorithm’). It was implemented in the context of a mechanistic approach to
gene regulatory network inference from single-cell data [5] and is based upon
an underlying stochastic dynamical model driven by the transcriptional
bursting phenomenon. In this tool paper, we introduce briefly the main
concepts behind the package, and detail its usage through some application
examples.

##  2 Theory

Consider a network of $n$ genes. Our starting point is the well-known ‘two-
state model’ of gene expression [11], which corresponds to the following set
of elementary chemical reactions [5] for each gene $i\in\\{1,\dots,n\\}$:

$\begin{array}[]{ccc}G_{i}\xrightleftharpoons[k_{\text{off},i}]{k_{\text{on},i}}{G_{i}}^{*},&{G_{i}}^{*}\xrightarrow[]{s_{0,i}}{G_{i}}^{*}+X_{i},&X_{i}\xrightarrow[]{d_{0,i}}\varnothing,\\\
&X_{i}\xrightarrow[]{s_{1,i}}X_{i}+Z_{i},&Z_{i}\xrightarrow[]{d_{1,i}}\varnothing,\end{array}$
(1)

where $G_{i}$, ${G_{i}}^{*}$, $X_{i}$ and $Z_{i}$ respectively denote
‘inactive promoter’, ‘active promoter’, mRNA and protein copy numbers for gene
$i$. These reactions describe the main two stages of gene expression, namely
_transcription_ (rate $s_{0,i}$) and _translation_ (rate $s_{1,i}$), along
with degradation of mRNA (rate $d_{0,i}$) and protein (rate $d_{1,i}$)
molecules. Note that $[G_{i}]+[{G_{i}}^{*}]=1$ is a conserved quantity: the
particularity of the two-state model is that transcription of gene $i$ can
only occur when its promoter is active, corresponding to $[{G_{i}}^{*}]=1$.

More specifically, experimental data consistently suggest a particular regime
for this model [13]: $k_{\text{off},i}\gg k_{\text{on},i}$ and $s_{0,i}\gg
d_{0,i}$ with the ratio $k_{\text{on},i}s_{0,i}/(k_{\text{off},i}d_{0,i})$
remaining fixed, corresponding to short active periods during which many mRNA
molecules are produced. In this regime, mRNA is transcribed by ‘bursts’ of
tens to hundreds of molecules. Moreover, mRNA and protein copy numbers are not
conserved quantities in the model and can be reasonably described in a
continuous way using standard ‘mass action’ kinetics. This leads to a hybrid
dynamical model, consisting of ordinary differential equations subject to
random jumps (Fig. 1) where the related quantities are denoted by
$X(t)=(X_{1}(t),\dots,X_{n}(t))\in\mathbb{R}_{+}^{n}$ and
$Z(t)=(Z_{1}(t),\dots,Z_{n}(t))\in\mathbb{R}_{+}^{n}$.

So far the genes are not interacting: the main point of our approach is to
consider the burst frequency $k_{\text{on},i}$ as a function of proteins
levels $Z_{1},\dots,Z_{n}$ [5]. In the current version of Harissa, this
function takes the following form:

$k_{\text{on},j}(z)=\frac{k_{0,j}+k_{1,j}\exp(\beta_{j}+\sum_{i=1}^{n}\theta_{ij}z_{i})}{1+\exp(\beta_{j}+\sum_{i=1}^{n}\theta_{ij}z_{i})}\qquad\forall
j\in\\{1,\dots,n\\}$ (2)

with $k_{0,j}\ll k_{1,j}$, so that $(\theta_{ij})_{1\leq i,j\leq n}$ can be
easily interpreted as the _network interaction matrix_ while $\beta_{j}$
encodes basal activity of gene $j$ (see also Table 1). The rate
$k_{\text{off},i}$ (inverse of mean burst duration) is kept constant.

Figure 1: The stochastic dynamical model underlying Harissa. (A) Complete
model with mRNA levels $X(t)$ and protein levels $Z(t)$, used to perform
simulations. Transcription of gene $i$ occurs in bursts at random times with
rate $k_{\text{on},i}(Z(t))$ and the burst size follows an exponential
distribution $\mathcal{E}(b_{i})$ where $b_{i}$ corresponds to a scaling
factor. (B) Reduced model involving only protein levels $Z(t)$. Without loss
of generality, the parameter $c_{i}$ is set to an arbitrary value as protein
levels are not measured. This model serves as a basis for deriving a pseudo-
likelihood whose maximization is at the core of inference.

Mathematically speaking, the complete set of variables $(X(t),Z(t))_{t\geq
0}\in\mathbb{R}_{+}^{n}\times\mathbb{R}_{+}^{n}$ forms a _piecewise-
deterministic Markov process_ (PDMP) that is fully characterized by its
(continuous) master equation:

$\begin{array}[]{rll}\displaystyle\frac{\partial}{\partial
t}p(x,z,t)&\displaystyle=\sum_{i=1}^{n}\left[d_{0,i}\frac{\partial}{\partial
x_{i}}\\{x_{i}p(x,z,t)\\}+d_{1,i}\frac{\partial}{\partial
z_{i}}\\{(z_{i}-x_{i})p(x,z,t)\\}\right.&\vspace{2mm}\\\
&\qquad\displaystyle+\left.k_{\text{on},i}(z)\left(\int_{0}^{x_{i}}p(x-he_{i},z,t)b_{i}e^{-b_{i}h}\mathrm{d}h-p(x,z,t)\right)\right],&\end{array}$
(3)

where $p(x,z,t)$ is the probability distribution of $(X(t),Z(t))$. The master
equation describes the time evolution of this probability and is related to
the trajectory dynamics (see Fig. 1A). Notably, each $Z_{i}$ is rescaled so
that parameters $s_{1,i}$ do not appear anymore, and $s_{0,i}$ and
$k_{\text{off},i}$ aggregate in the bursty regime into the ‘burst size’
parameter $1/b_{i}=s_{0,i}/k_{\text{off},i}$ (Table 1).

Table 1: Parameters of the dynamical model underlying Harissa. Here we
consider an instance model = NetworkModel(n) where n is the number of genes in
the network.

_Degradation kinetics_
---
_Bursting kinetics_
_Network parameters_

Package variable | Notation |  |  Interpretation (gene $i$)
---|---|---|---
model.d[0][i] | $d_{0,i}$ |  | mRNA degradation rate
model.d[1][i] | $d_{1,i}$ |  | protein degradation rate
model.a[0][i] | $k_{0,i}$ |  | minimal burst frequency
model.a[1][i] | $k_{1,i}$ |  | maximal burst frequency
model.a[2][i] | $b_{i}$ |  | inverse of mean burst size
model.basal[i] | $\beta_{i}$ |  | basal activity
model.inter[i,j] | $\theta_{ij}$ |  |  interaction $i\to j$

This ‘complete’ model is simulated in Harissa using an efficient acceptance-
rejection method, similarly to the explicit construction given in [1]. A great
advantage of this method is that it is guaranteed to be exact without
requiring any numerical integration, contrary to the basic algorithm [6].

We also consider a ‘protein-only’ model (Fig. 1B), which can be seen as a
first-order approximation (Appendix A.1): this reduced model turns out to
provide a tractable inference procedure based on analytical results (Appendix
A.2).

##  3 Usage

The Harissa package has two main functionalities: network inference
interpreted as calibration of the stochastic dynamical model defined by (3),
and data simulation from the same model (e.g., scRNA-seq counts [10] but also
RT-qPCR levels [7, 12]). Besides, the package also allows for basic network
visualization (directed graphs with positive or negative edge weights) as well
as data binarization (using gene-specific thresholds derived from the data-
calibrated dynamical model).

The first step is to create an instance of the model:

model = NetworkModel(n)

where n is the number of genes. After setting the parameter values (see Table
1), the model can be simulated using model.simulate(time), where time is
either a single time or a list of time points.

Importantly, simulate() does not depend on time discretization and always
returns exact stochastic simulations: the resulting continuous-time
trajectories are simply extracted at user-specified time points. It is also
possible to consider a stimulus, represented by an additional protein $Z_{0}$
that receives no feedback and verifies $Z_{0}(t)=0$ for $t\leq 0$ and
$Z_{0}(t)=1$ for $t>0$. In order to reach a pre-stimulus steady state before
perturbation, an optional burnin parameter sets the time during which the
model is simulated with $Z_{0}(t)=0$.

#### Example: Repressilator Network.

As an example, we consider a ‘repressilator’ network made of 3 genes forming a
directed cycle of negative interactions (Fig. 2).

Figure 2: The repressilator network used in Fig. 3. This plot was made with
plot_network from the harissa.utils module.

Some critical parameters of the dynamical model are degradation rates
$d_{0,i}$ and $d_{1,i}$, which characterize the ‘responsiveness’ of mRNA and
protein levels. Here we set $d_{0,i}/d_{1,i}=10$, meaning that proteins are
ten times more stable than mRNAs. This is biologically realistic, but note
that this ratio is known to span a very wide range [9] so there is no single
choice.

An example of simulated single-cell trajectory for this network is shown in
Fig. 3. It is worth noticing that despite the strong level of stochasticity, a
robust periodic pattern is already emerging.

More stable proteins—with respect to mRNA—will lead to less ‘intrinsic noise’
in the system. We increase mRNA degradation rates and burst frequencies
instead, which is equivalent to a zoom-out regarding the time scale. Since
mRNA and protein scales are normalized, the overall levels do not depend on
degradation rates (but the _dynamics_ does).

The stochastic model converges as $d_{0}/d_{1}\to\infty$ to a slow-fast limit
[2] which turns out to be a nonlinear ODE system involving only proteins
(bottom plot of Fig. 3). It appears here that due to the deterministic
dynamics, the initial protein levels need to be perturbed so as not to stay in
a trivial manifold. This limit model can be simulated with
model.simulate_ode(time) and is generally useful to gain insight into the
system attractors (note however that it is a rough approximation of the
stochastic model, see remark 1).

In this slow-fast limit, which is part of the rationale behind the reduced
model (Appendix A.1), mRNA levels become _independent conditionally on protein
levels_ such that
$X_{i}(t)\sim\mathrm{Gamma}(k_{\mathrm{on},i}(Z(t))/d_{0,i},b_{i})$ for
$i=1,2,3$. This quasi-steady-state (QSS) behavior can be understood
intuitively from the top plot of Fig. 3. Regarding mRNA levels, simulate_ode
only returns the mean of the QSS distribution conditionally on protein levels
(the true limit model would consist in sampling from this distribution
independently for every $t>0$).

Figure 3: Simulation of a repressilator network using Harissa. The first two
plots show an example of single-cell trajectory of mRNA levels $X(t)$ and
protein levels $Z(t)$ from the full stochastic model (3). The bottom plot
shows a trajectory of $Z(t)$ from the related deterministic model
corresponding to the slow-fast limit (very stable proteins compared to mRNA,
i.e. $d_{0}/d_{1}\to\infty$). Importantly, the periodic pattern is emerging
and robust well before the limit is reached (here $d_{0}/d_{1}=10$). The burst
parameters (Table 1) are $k_{0,i}=0$, $k_{1,i}=2$ and $b_{i}=0.02$ for all
$i\in\\{1,2,3\\}$, which are current default values. The network parameters
are given in Appendix A.3.

#### Network Inference.

Here the main function is model.fit(), which takes as input a time-course
single-cell dataset (Fig. 4). Inference can be performed by creating a new
instance model = NetworkModel() without any size parameter, then loading a
dataset x and using model.fit(x). This will update all parameters of the model
(Table 1) except $d_{0,i}$ and $d_{1,i}$, which need to be provided by the
user from external data—or left to their default values, see remark 1—as they
cannot be inferred without seriously compromising identifiability of other
parameters. As an important feature, the updated model instance is ready for
simulation, to assess reproducibility of the original data or to predict the
outcome of network modifications. From a network inference viewpoint, the only
important parameter is model.inter which is a non-symmetric signed weight
matrix.

######  Remark 1.

Degradation rates $d_{0,i}$ and $d_{1,i}$ are in fact not required for the
_inference_ part, but they are important for the _simulation_ part. They are
currently set by default to $d_{0,i}=\ln(2)/9\approx 0.077\;\mathrm{h}^{-1}$
and $d_{1,i}=\ln(2)/46\approx 0.015\;\mathrm{h}^{-1}$ as median values from
thousands of genes [9]. This leads in particular to $d_{0,i}/d_{1,i}\approx
5.11$.

$\cdots$_cells_ _genes_ $t_{0}$$t_{1}$$t_{2}$$t_{3}$$t_{K}$

Figure 4: Typical structure of time-course scRNA-seq data, as required by the
inference module (here each row is a gene and each column is a cell). Each
group of cells collected at the same experimental time $t_{k}$ forms a
_snapshot_ of the biological heterogeneity at time $t_{k}$. Due to the
destructive nature of the measurement process, snapshots are made of different
cells. This data is therefore different from so-called ‘pseudotime’
trajectories, which attempt to reorder cells according to some smoothness
hypotheses.

##  4 Conclusions

While the development of Harissa started a few years ago, the tool has
recently become more mature and easy to use. Following early versions [5, 3],
an alternative method for the inference part called Cardamom was developed in
[14], which in turn influenced the current version of Harissa. The two
inference methods remain complementary at this stage and may be merged into
the same package in future development. They were evaluated in a recent
benchmark [16], giving encouraging results and showing, most importantly, that
the stochastic dynamical model (3) can indeed reproduce real time-course
scRNA-seq data through gene interactions. While the simulation part is quite
well optimized and can be considered as stable, the inference part is much
more challenging. Future directions include a better use of the dynamical
information contained in time-course transcriptional profiles, which has great
potential for performance improvement.

#### Code Availability.

The code of the package, a tutorial and some basic usage scripts are available
at https://github.com/ulysseherbach/harissa. In addition, Harissa is indexed
in the Python Package Index and can be installed via pip.

##  Appendix A Appendices

###  A.1 Reduced Model

The inference procedure is based on analytical results which are not available
for the two-stage ‘mRNA-protein’ model (3). On the other hand, such results
exist for a one-stage ‘protein-only’ model that is a valid approximation of
the former when proteins are more stable than mRNA (i.e. $d_{0,i}/d_{1,i}\gg
1$). The resulting process $(Z(t))_{t\geq 0}\in\mathbb{R}_{+}^{n}$ is also a
PDMP, whose master equation can be interpreted in terms of simplified
trajectories (Fig. 1B):

$\begin{array}[]{rll}\displaystyle\frac{\partial}{\partial
t}p(z,t)&\displaystyle=\sum_{i=1}^{n}\left[d_{1,i}\frac{\partial}{\partial
z_{i}}\\{z_{i}p(z,t)\\}\right.&\vspace{2mm}\\\
&\qquad\displaystyle+\left.\int_{0}^{z_{i}}k_{\text{on},i}(z-he_{i})p(z-he_{i},t)c_{i}e^{-c_{i}h}\mathrm{d}h-k_{\text{on},i}(z)p(z,t)\right].&\end{array}$
(4)

Given $Z(t)$, mRNA levels $X(t)$ are obtained by sampling independently for
every $i\in\\{1,\dots,n\\}$ and $t>0$ from
$X_{i}(t)\sim\mathrm{Gamma}(k_{\mathrm{on},i}(Z(t))/d_{0,i},b_{i}$), which is
the _quasi-steady-state_ (QSS) distribution of the complete model [5, 6].

###  A.2 Inference Algorithm

Now consider mRNA counts measured in $m$ cells, assumed independent, along a
time-course experiment following a stimulus. Each cell $k=1,\dots,m$ is
associated with an experimental time point $t_{k}$. We introduce the following
notation:

$\mathbf{x}_{k}=(x_{ki})\in\\{0,1,2,\dots\\}^{n}$ : mRNA counts (cell $k$,
gene $i$);

$\mathbf{z}_{k}=(z_{ki})\in(0,+\infty)^{n}$ : latent protein levels (cell $k$,
gene $i$);

$\alpha=(\alpha_{ij}(t_{k}))\in\mathbb{R}^{n\times n}$ : effective interaction
$i\to j$ at time $t_{k}$.

A stimulus is represented as gene $i=0$ and we therefore add parameters
$\alpha_{0j}(t_{k})$ for $j=1,\dots,n$ and $k=1,\dots,m$. We further set
$z_{k0}=0$ if $t_{k}\leq 0$ (before stimulus) and $z_{k0}=1$ if $t_{k}>0$
(after stimulus). Then, writing $a_{i}=k_{1,i}/d_{0,i}$, the underlying
statistical model of Harissa is defined by

$\displaystyle p(\mathbf{z}_{k})$
$\displaystyle=\prod_{i=1}^{n}{z_{ki}}^{c_{i}\sigma_{ki}-1}e^{-c_{i}z_{ki}}\frac{{c_{i}}^{c_{i}\sigma_{ki}}}{\mathrm{\Gamma}(c_{i}\sigma_{ki})}\;,$
(5) $\displaystyle p(\mathbf{x}_{k}|\mathbf{z}_{k})$
$\displaystyle=\prod_{i=1}^{n}\frac{1}{x_{ki}!}\frac{\mathrm{\Gamma}(a_{i}z_{ki}+x_{ki})}{\mathrm{\Gamma}(a_{i}z_{ki})}\frac{{b_{i}}^{a_{i}z_{ki}}}{(b_{i}+1)^{a_{i}z_{ki}+x_{ki}}}\;,$
(6)

with

$\sigma_{ki}={\left[1+\exp(-\\{\beta_{i}+\alpha_{0i}(t_{k})z_{k0}+\textstyle\sum_{j=1}^{n}\alpha_{ji}(t_{k})z_{kj}\\})\right]}^{-1}.$
(7)

Details of this derivation can be found in [5, 14, 4, 8]. Roughly, (5) comes
from a ‘Hartree’ approximation of (4), while (6) corresponds to a Poisson
distribution with random parameter sampled from the QSS distribution of $X(t)$
given $Z(t)$. Note that $p(\mathbf{z}_{k})$ is in general only a pseudo-
likelihood as $\sigma_{ki}$ depends on $\mathbf{z}_{k}$.

Since the preliminary version of Harissa [5], the global inference procedure
has been heavily improved using important identifiability results from [14,
15]. The final algorithm consists of three steps:

  1. 1.

_Model calibration:_ estimate $a_{i}$ and $b_{i}$ for each gene individually
from (6);

  2. 2.

_Bursting mode inference:_ estimate the frequency mode ($k_{0,i}$ or
$k_{1,i}$) for each gene in each cell (can be seen as a binarization step with
specific thresholds);

  3. 3.

_Network inference:_ consider $\mathbf{z}_{k}$ as observed from step 2 and
maximize (5) with respect to $\alpha$ after adding an appropriate penalization
term [14]. Each parameter $\theta_{ij}$ is then set to $\alpha_{ij}(t_{k})$
with $t_{k}$ that maximizes $|\alpha_{ij}(t_{k})|$.

###  A.3 Repressilator Network

Considering an instance model = NetworkModel(3), the repressilator network
simulated in Fig. 3 is defined as follows:

    model.d[0] = 1 # mRNA degradation rates
    model.d[1] = 0.1 # Protein degradation rates
    model.basal[1] = 5 # Basal activity of gene 1
    model.basal[2] = 5 # Basal activity of gene 2
    model.basal[3] = 5 # Basal activity of gene 3
    model.inter[1,2] = -10 # Interaction 1 -> 2
    model.inter[2,3] = -10 # Interaction 2 -> 3
    model.inter[3,1] = -10 # Interaction 3 -> 1

#### Acknowledgements.

The author is very grateful to Elias Ventre and Olivier Gandrillon for
fruitful discussions which led to improve the Harissa package.

## References

  * Benaïm et al. [2015] M. Benaïm, S. Le Borgne, F. Malrieu, and P.-A. Zitt.  Qualitative properties of certain piecewise deterministic Markov processes.  _Annales de l’Institut Henri Poincaré - Probabilités et Statistiques_ , 51(3):1040–1075, 2015.  https://doi.org/10.1214/14-AIHP619.
  * Faggionato et al. [2010] A. Faggionato, D. Gabrielli, and M. Crivellari.  Averaging and large deviation principles for fully-coupled piecewise deterministic Markov processes and applications to molecular motors.  _Markov Processes and Related Fields_ , 16(3):497–548, 2010.  https://doi.org/10.48550/arXiv.0808.1910.
  * Herbach [2018] U. Herbach.  _Modélisation stochastique de l’expression des gènes et inférence de réseaux de régulation_.  PhD thesis, Université de Lyon, 2018.
  * Herbach [2019] U. Herbach.  Stochastic gene expression with a multistate promoter: breaking down exact distributions.  _SIAM Journal on Applied Mathematics_ , 79(3):1007–1029, 2019.  https://doi.org/10.1137/18M1181006.
  * Herbach et al. [2017] U. Herbach, A. Bonnaffoux, T. Espinasse, and O. Gandrillon.  Inferring gene regulatory networks from single-cell data: a mechanistic approach.  _BMC Systems Biology_ , 11(1):105, 2017.  https://doi.org/10.1186/s12918-017-0487-0.
  * Malrieu [2015] F. Malrieu.  Some simple but challenging Markov processes.  _Annales de la Faculté de Sciences de Toulouse_ , 24(4):857–883, 2015.  https://doi.org/10.5802/afst.1468.
  * Richard et al. [2016] A. Richard, L. Boullu, U. Herbach, A. Bonnafoux, V. Morin, E. Vallin, A. Guillemin, N. Papili Gao, R. Gunawan, J. Cosette, O. Arnaud, J.-J. Kupiec, T. Espinasse, S. Gonin-Giraud, and O. Gandrillon.  Single-cell-based analysis highlights a surge in cell-to-cell molecular variability preceding irreversible commitment in a differentiation process.  _PLOS Biology_ , 14(12):e1002585, 2016.  https://doi.org/10.1371/journal.pbio.1002585.
  * Sarkar and Stephens [2021] A. Sarkar and M. Stephens.  Separating measurement and expression models clarifies confusion in single-cell RNA sequencing analysis.  _Nature Genetics_ , 53(6):770–777, 2021.  https://doi.org/10.1038/s41588-021-00873-4.
  * Schwanhäusser et al. [2011] B. Schwanhäusser, D. Busse, N. Li, G. Dittmar, J. Schuchhardt, J. Wolf, W. Chen, and M. Selbach.  Global quantification of mammalian gene expression control.  _Nature_ , 473(7347):337–342, 2011.  https://doi.org/10.1038/nature10098.
  * Semrau et al. [2017] S. Semrau, J. E. Goldmann, M. Soumillon, T. S. Mikkelsen, R. Jaenisch, and A. van Oudenaarden.  Dynamics of lineage commitment revealed by single-cell transcriptomics of differentiating embryonic stem cells.  _Nature Communications_ , 8(1):1096, 2017.  https://doi.org/10.1038/s41467-017-01076-4.
  * Shahrezaei and Swain [2008] V. Shahrezaei and P. S. Swain.  The stochastic nature of biochemical networks.  _Current Opinion in Biotechnology_ , 19(4):369–374, 2008.  https://doi.org/10.1016/j.copbio.2008.06.011.
  * Stumpf et al. [2017] P. S. Stumpf, R. C. G. Smith, M. Lenz, A. Schuppert, F.-J. Müller, A. Babtie, T. E. Chan, M. P. H. Stumpf, C. P. Please, S. D. Howison, F. Arai, and B. D. MacArthur.  Stem cell differentiation as a non-Markov stochastic process.  _Cell Systems_ , 5(3):268–282, 2017.  https://doi.org/10.1016/j.cels.2017.08.009.
  * Tunnacliffe and Chubb [2020] E. Tunnacliffe and J. R. Chubb.  What is a transcriptional burst?  _Trends in Genetics_ , 36(4):288–297, 2020.  https://doi.org/10.1016/j.tig.2020.01.003.
  * Ventre [2021] E. Ventre.  Reverse engineering of a mechanistic model of gene expression using metastability and temporal dynamics.  _In Silico Biology_ , 14(3-4):89–113, 2021.  https://doi.org/10.3233/ISB-210226.
  * Ventre et al. [2021] E. Ventre, T. Espinasse, C.-E. Bréhier, V. Calvez, T. Lepoutre, and O. Gandrillon.  Reduction of a stochastic model of gene expression: Lagrangian dynamics gives access to basins of attraction as cell types and metastabilty.  _Journal of Mathematical Biology_ , 83(5):59, 2021\.  https://doi.org/10.1007/s00285-021-01684-1.
  * Ventre et al. [2023] E. Ventre, U. Herbach, T. Espinasse, G. Benoit, and O. Gandrillon.  One model fits all: Combining inference and simulation of gene regulatory networks.  _PLOS Computational Biology_ , 19(3):e1010962, 2023.  https://doi.org/10.1371/journal.pcbi.1010962.


# Trigonal Symmetry Breaking and its Electronic Effects in Two-Dimensional
Dihalides and Trihalides

Alexandru B. Georgescu<EMAIL_ADDRESS>Department of
Materials Science and Engineering, Northwestern University, Evanston, Illinois
60208, USA  Andrew J. Millis<EMAIL_ADDRESS>Center for Computational
Quantum Physics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA
Department of Physics, Columbia University, 538 West 120th Street, New York,
New York 10027, USA  James M. Rondinelli<EMAIL_ADDRESS>Department of Materials Science and Engineering, Northwestern University,
Evanston, Illinois 60208, USA

###### Abstract

We study the consequences of the approximately trigonal ($D_{3d}$) point
symmetry of the transition metal (M) site in two-dimensional van der Waals MX2
dihalides and MX3 trihalides. The trigonal symmetry leads to a 2-2-1 orbital
splitting of the transition metal $d$ shell, which may be tuned by the
interlayer distance, and changes in the ligand-ligand bond lengths. Orbital
order coupled to various lower symmetry lattice modes may lift the remaining
orbital degeneracies, and we explain how these may support unique electronic
states using ZrI2 and CuCl2 as examples, and offer a brief overview of
possible electronic configurations in this class of materials. By building and
analysing Wannier models adapted to the appropriate symmetry we examine how
the interplay among trigonal symmetry, electronic correlation effects, and
$p$-$d$ orbital charge transfer leads to insulating, orbitally polarized
magnetic and/or orbital-selective Mott states. Our work establishes a rigorous
framework to understand, control, and tune the electronic states in low-
dimensional correlated halides. Our analysis shows that trigonal symmetry and
its breaking is a key feature of the 2D halides that needs to be accounted for
in search of novel electronic states in materials ranging from CrI3 to
$\alpha$-RuCl3.

##  I Introduction

Transition metal compounds exhibit electronic properties of high scientific
and technological interest, including ferroelectricity [1, 2, 3, 4], quantum
magnetism [5, 6, 7, 8, 9, 10], metal-insulator transitions [11, 12, 13, 14,
15, 16, 17, 18, 19, 20, 21] , and high transition-temperature
superconductivity [22, 23, 24, 25, 26]. Transition metal oxides derived from
the AMO3 perovskite structure have been a focus of particular attention
because any $3d$ or $4d$ transition metal can occupy the M site with
(typically) partially filled $d$ shells, while variation of the A-site ion can
tune the relative valence of the M site ion and the electronic bandwidth. The
perovskite structure is also highly polymorphic.

It allows for many variants of the basic structure that exhibit different
crystallographic symmetries, which activate interesting electronic states and
the pseudocubic structure allows for a wide variety of superlattices to be
built [19, 18]. Basic to the electronic physics of perovskites is the cubic
($O_{h}$) point symmetry of the M-site ion and its reduction to tetragonal
symmetry by even parity octahedral distortions.

Recently, two-dimensional (2D) van der Waals transition metal MX2 dihalides
and trihalides MX3 with X a halogen ligand have become of interest as they
exhibit layer-dependent ferromagnetism as in VI3 [27, 28, 29, 30, 31, 32],
possible Kitaev spin liquid behavior in RuCl3 [33, 34, 35, 36, 37, 38] and
other magnetic phenomena [39, 40, 41]. Additionally, the crystal and
electronic structures are highly two dimensional, so the materials can be
exfoliated, made in monolayer form, doped by gating and layering with other
compounds, and potentially twisted into Moiré materials [42]. Similar to the
perovskite transition metal oxides, the M site is six-fold coordinated by
ligands and can host essentially all $3d$ and $4d$ transition metals often
with partial $d$-orbital occupancy, which endows them with the aforementioned
physical properties. Dihalide and trihalide compounds with the same transition
metal will exhibit different nominal valence and octahedral coordination.
Table 1 lists structural and electronic features of several known 2D halides.

Although the sixfold coordination of the metal cations in the dihalides and
trihalides would suggest octahedral symmetry $O_{h}$, analogous to the
octahedral coordination in many perovskite oxides, the halides are more
appropriately described as trigonally coordinated. A local orbital basis
derived from $e_{g}$ and $t_{2g}$ representations of the $O_{h}$ group is not
the most useful description [43].

As Table 1 shows, the metal site symmetry in the bulk halides is either
trigonal ($D_{3d},D_{3},C_{3i}$) or a trigonal subgroup ($C_{2h}$) or even
lower [44] and as previously noted for MX2 compounds [45] a basis that
conforms to the symmetry enables a more straightforward treatment of the
physics.

Table 1:  Symmetry, atomic structure, $d$-electron configuration for the metal, and magnetic order of experimental transition metal dihalides and trihalides. AFM, FM, and HM indicate antiferromagnetic, ferromagnetic, and helimagnetic, respectively. NA and MV indicate not applicable and multiple values for the specified lengths, respectively. Structural data from Ref. [46] and magnetic data from Ref. [44]. Structure type | CdI2 | CdCl2 | BiI3 | CrCl3 | AlCl3
---|---|---|---|---|---
X packing | _hcp_ | _ccp_ | _hcp_ | _ccp_ |  $\sim$_ccp_
Space group | $P\bar{3}m1$ | $R\bar{3}m$ | $R\bar{3}$ | $P3_{2}12$ | $C2/m$
Point group | $D_{3d}$ | D3d | $C_{3i}$ | $D_{3}$ | $C_{2h}$
Compound | FeI2 | NiI2 | FeCl3 | CrCl3 | CrI3
M-$d^{n}$ config. | $d^{6}$ | $d^{8}$ | $d^{5}$ | $d^{3}$ | $d^{3}$
M-site sym. | $D_{3d}$ | $D_{3d}$ | $C_{3}$ | $C_{3}$ | $C_{2}$
Magnetism | AFM | HM | HM | AFM | FM
$L_{1}$ | 4.050 | 3.927 | 3.380 | 3.483 | MV
$L_{2}$ | 4.110 | 4.580 | 3.216 | 3.228 | MV
$L_{3}$ | N/A | N/A | 3.472 | 3.229 | MV
Compound | TiCl2 | MnCl2 | VCl3 |  $\alpha$-RuCl3 |  $\alpha$-RuCl3
M-$d^{n}$ config. | $d^{2}$ | $d^{5}$ | $d^{2}$ | $d^{5}$ | $d^{5}$
M-site sym. | $D_{3d}$ | $D_{3d}$ | $C_{3}$ | $C_{3}$ | $C_{2}$
Magnetism | AFM | AFM, HM | AFM | AFM | -
$L_{1}$ | 3.430 | 3.711 | 3.471 | 3.44 | MV
$L_{2}$ | 3.636 | 3.505 | 3.366 | 3.490 | MV
$L_{3}$ | NA | NA | 3.366 | 3.487 | MV

In this paper we present a general analysis of the local electronic structure
of the transition metal $d$ shells in dihalides and trihalides, considering
both ideal structures and consequences of cooperative atomic displacements. By
analogy to the more frequently studied perovskite compounds, we identify what
correlated electron behavior may arise. We focus on the implications of the
trigonal point symmetry and answer the following questions: What is the
appropriate orbital description for the correlated electrons? How to tune the
orbital structure by varying inter-layer distances and local MX6 geometry,
_i.e._ , ligand-ligand and metal-ligand distances? What are the consequences
of symmetry reductions for the correlated electron phases.

In addressing these questions, it is often helpful to construct a Wannier
basis that includes a representation of the local transition metal $d$
orbitals, which host the electron interactions associated with correlated
electron phenomena. However, standard Wannierization procedures for these
materials, do not easily produce a basis that transforms properly under
trigonal symmetry operations. We present a procedure for constructing an
appropriate symmetry-adapted basis, which can easily be adapted to other
materials. This allows us to discuss the correlated insulating states that may
be obtained in stoichiometric compounds, as well as new electronic states that
may appear as a result of further symmetry breaking, such as those present in
the semi-1D chains in ZrI2, and certain polymorphs of RuCl3, CrI3 and other
halides [44].

##  II Local structure: electronic states and symmetry considerations

The basic structural unit of 2D halides is a plane of transition metal ions,
with each transition metal ion coordinated with 6 halogen ions (Figure 1A).
Unlike in the layered perovskite-based transition metal compounds, the local
axes of the octahedron surrounding a transition metal ion in a halide are
rotated with respect to the axes that define the two-dimensional plane.
Specifically, choosing coordinates such that the monolayer plane is
perpendicular to [001] then the octahedral axes (M-X bond directions) are
[$\pm 1\pm 1\pm 1$] in the ideal halides (compare panels B and C of Figure 1).
In the ideal monolayer case, the point symmetry of the M ion is then $D_{3d}$.
The symmetry may be further reduced by spin-orbit coupling, additional
relaxations of the atomic positions, magnetic order and by different stackings
of the layered structures. We discuss some lattice symmetry reductions, as
well as magnetic order later.

Figure 1: A) Sequential illustration of symmetry breaking as a result of the
layered structure in dihalides and trihalides. Placing an ideal octahedron
within a monolayer will reduce its Oh symmetry to trigonal $D_{3d}$ symmetry:
3-fold rotations along the z axis are preserved, and so is mirror and
inversion symmetry - even without distortions to the octahedron itself. Sketch
of B) tetragonal and C) trigonal symmetry breaking of an octahedron. Heavy
magenta arrows indicate vector normal to monolayer planes. Layer structure of
MX6 octahedra for D) MX2 dihalides and E) MX3 trihalides. Definitions for
halide-halide lengths for octahedra in F) MX2 and G) MX3 halides. The red
edges with length $L_{1}$ surround the face parallel to the monolayer planes,
$L_{2}$ are edges that are not parallel to the material plane and that connect
nearest-neighbor MX6 octahedra. $L_{3}$ edges do not connect octahedra, are
not parallel to the plane, and only appear in MX3 compounds.

The key feature of the $D_{3d}$ symmetry is a three-fold rotation axis, in the
ideal case perpendicular to the monolayer plane of MX6 octahedra (Figure 1D,E)
and passing through a triangular X3 face (Figure 1F,G).

In $D_{3d}$ symmetry, the M-X bond lengths remain equal, however, the X-X
ligand distances $L_{i}$ ($i=1-3$) need not be equal. The inequivalence is
accommodated through changes in the X-M-X intra-octahedral angles. Similar
effects are known in edge sharing perovskites [47, 48] Because of the planar
structure formed by edge connectivity of the octahedra, the halogen-halogen
distances perpendicular to $(001)$ (denoted $L_{1}$ in panels F and G of
Figure 1) are inequivalent to X-X distances on other triangular faces.

This can be plainly observed in the MX2 compounds, where X-X edges not
parallel to the plane connect different octahedra, while those that are
parallel to the plane connect the same octahedra. This leads to the
inequivalence of the $L_{1}$ and $L_{2}$ X-X edges (Figure 1F).

One may parametrize the trigonal distortion away from perfect octahedral
coordination either by the $L_{2}/L_{1}$ ratio or by the X-M-X bond angles. We
find that the $L_{2}/L_{1}$ ratio is a more convenient parametrization in
$D_{3d}$ symmetry, as it plays an analogous role to the $c/a$ ratio for the
tetragonal symmetry and tunes the splitting between the doublet and singlet.

For the MX3 compounds only half of the X-X edges not parallel to the plane
connect MX6 octahedra, leading to a total of three inequivalent edges,
$L_{1}$, $L_{2}$, and $L_{3}$. Figure 1D illustrates that in both cases, there
are only two types of inequivalent faces: those parallel to the planes of MX6,
with three $L_{1}$ edges, and those oblique to the plane with edges $L_{1}$,
$L_{2}$ and $L_{3}$.

While the symmetry of an ideal halide monolayer is $D_{3d}$, different
arrangements of the planes can lead to symmetry reductions to $D_{3}$ or
$C_{3i}$ in the bulk compounds (see Table 1). These subgroups, however, lead
to the same trigonal orbital basis for the transition metal M $d$ orbitals.
Further, anisotropic semicovalent and van der Waals interactions or additional
electronically driven orderings can lead to further symmetry reductions (see
Table 1). In such cases the clearest way to understand the resulting orbital
structure is as an additional symmetry reduction beyond trigonal.

Table 2:  Relationships among the atomic $d$ orbitals for the standard cubic (tetragonal) basis with $z$ aligned to the (001) direction shown in Figure 1B and trigonal basis with $z$ axis aligned to (111) direction shown in Figure 1C. Note that in the absence of trigonal symmetry breaking, the $a_{1g}$ and e${}_{g}^{\sigma}$ orbitals become degenerate and combine to transform as the $t_{2g}$ representation.  Cubic Basis [001] | Trigonal Basis [111]
---|---
Symmetry | Orbitals | Symmetry | Orbitals
$t_{2g}$ | $d_{\bar{x}\bar{y}}$ | $a_{1g}$ | $d_{z^{2}}$
| $d_{\bar{y}\bar{z}}$ | $e_{g}^{\pi}$ | $\frac{2}{\sqrt{6}}d_{xy}+\frac{2}{\sqrt{3}}d_{yz}$
| $d_{\bar{x}\bar{z}}$ |  | $\frac{2}{\sqrt{6}}d_{x^{2}-y^{2}}-\frac{1}{\sqrt{3}}d_{xz}$
$e_{g}$ | $d_{\bar{x}^{2}-\bar{y}^{2}}$ | $e_{g}^{\sigma}$ | $\frac{1}{\sqrt{3}}d_{x^{2}-y^{2}}+\frac{2}{\sqrt{6}}d_{xz}$
| $d_{\bar{z}^{2}-r^{2}}$ |  | $\frac{1}{\sqrt{3}}d_{xy}-\frac{2}{\sqrt{6}}d_{yz}$

The trigonal symmetry has implications for the electronic structure. Under
$D_{3d}$ symmetry, the 5 $d$ orbitals transform as two doublets
($e_{g}^{\sigma}$ and $e_{g}^{\pi}$) and a singlet ($a_{1g}$). The wave
functions and level splittings are not constrained by symmetry and depend on
details of atomic-scale physics. The trigonal structure of the halides may be
viewed as a weak distortion of the $O_{h}$ symmetry familiar from cubic
perovskites, enabling a simpler interpretation of the basis functions. Under
$O_{h}$ symmetry, the atomic $d$ orbitals transform as an $e_{g}$ symmetry
doublet (wavefunctions conventionally chosen as $d_{\bar{x}^{2}-\bar{y}^{2}}$
and $d_{3\bar{z}^{2}-r^{2}}$) and a $t_{2g}$ triplet (wavefunctions
conventionally chosen as
$d_{\bar{x}\bar{y}},d_{\bar{x}\bar{z}},d_{\bar{y}\bar{z}}$) where
$\bar{x},\bar{y},\bar{z}$ are the three octahedral axes. Ligand fields arising
from hybridization lead to an energy separation between the doublet and
triplet states of the order of 2 eV. The further reduction of the symmetry
from $O_{h}$ to $D_{3d}$ does not additionally split the $e_{g}$ states, which
now form the $e_{g}^{\sigma}$ doublet representation of $D_{3d}$. It does,
however, split the $t_{2g}$ triplet into an $e_{g}^{\pi}$ doublet and an
$a_{1g}$ singlet with the $e_{g}^{\pi}-a_{1g}$ level splitting being typically
smaller than the energetic separation to the $e_{g}^{\sigma}$ doublet. Then to
first order in the trigonal distortion, the basis functions for the
$e_{g}^{\sigma}$ representation are linear combinations of the familiar cubic-
basis $e_{g}$ states while the basis functions for the $e_{g}^{\pi}$ and
$a_{1g}$ states are linear combinations of the familiar $t_{2g}$ states.
(Table 2) 111Note that for the cubic case, the conventional coordinate system
involves a $z$ axis chosen parallel to an M-X bond while for the trigonal case
the conventional coordinate system involves a $z$ axis passing through an
octahedral face (_i.e._ , along $[111]$ in conventional octahedral
coordinates). . Lowering the symmetry below trigonal, in particular by
breaking the $C_{3}$ rotational symmetry about the axis perpendicular to the
plane, will lift the degeneracies of the two $e_{g}$ doublets. Such
distortions may occur if the transition metal valence is such that one of the
orbital pairs forming a doublet are partially filled, enabling an electronic
symmetry-breaking transition.

##  III Physics of the Trigonal Distortion in TiCl2

Figure 2:  Electronic orbital splitting at the non-magnetic DFT level for
TiCl2 as a function of (A) Ti-Ti interlayer distance and (B) $L_{2}/L_{1}$
ratio, as quantified by crystal field splitting at the $\Gamma$ point and
orbital occupancy as obtained from the density matrix of the orbital
projections. In A, the octahedron is kept perfect, with a X-X distance of 3.43
Å. In B, the inter-layer distance is kept fixed to the experimental value,
while $L_{2}/L_{1}$ is varied. The experimental inter-layer distance and
$L_{2}/L_{1}$ ratio are indicated with broken magenta lines.

We now present density functional theory (DFT) calculations on the
representative dihalide TiCl2 222DFT calculations were performed using the
Quantum Espresso software package on a 3 atom formula unit, using ultrasoft
pseudopotentials using the PBE exchange-correlation functional, an energy
cutoff of 748eV and a 20x20x20 k-grid. Orbital occupations were obtained using
the orthonormal atomic projections, as implemented in Quantum Espresso; DFT+U
calculations were performed with the U applied to orthonormal atomic
projections as well, within the Dudarev simplified +U formalism., using the
experimental structure obtained from ICSD [46]. TiCl2 exhibits a $d^{2}$
electronic configuration and its primitive $P\bar{3}m1$ structure contains a
single formula unit. The Ti atoms are aligned along $z$, making for the
simplest possible stacking. The Ti site symmetry is $D_{3d}$. Details of its
crystal structure are presented in Table 1.

We begin by analyzing non-magnetic DFT calculations using the PBE functional,
in order to isolate the effect of the atomic structure on the on-site crystal
field splittings and orbital order before performing further analysis. We then
investigate correlation effects in the DFT+U approximation. This stepwise
approach allows us to disentangle the effects of the structure from that of
$d-$shell electron-electron interactions. We select different interlayer
distances, X-X distances, and $L_{2}/L_{1}$ ratios to probe the trigonal
symmetry effects on the orbital structure. We characterize the electronic
structure changes via the on-site energy of the orbitals at the $\Gamma$
point, and the orbital occupations, as obtained from the eigenvalues of the
density matrix defined from orbital projectors.

Experimentally, these perturbations to the inter-layer distance may be
realized via external pressure, while the $L_{2}/L_{1}$ ratio can be tuned via
epitaxial strain. The effects we find pertain to MX3 compounds as well, as the
MX3 structure can be obtained from the MX2 structure by removing 1/3 of the M
atoms and keeping the ligand octahedral structure intact.

###  III.1 Level Splittings and Orbital Occupancies

Figure 2 presents a summary of our results.

Starting with the experimental structure, we find that the energy splitting at
$\Gamma$ corresponds to the expected 2-2-1 orbital splitting as shown by the
vertical broken line in Figure 2A,B.

Upon visualizing isosurfaces of the wavefunctions in real space at the
$\Gamma$ point (not shown), we find the orbitals represent wavefunctions of
the trigonal basis.

These wavefunctions are similar to the Wannier functions presented in Figure
3.

The lowest energy Ti $d$-derived state at $\Gamma$ is the $a_{1g}$ state. The
wavefunction exhibits lobes directed through the faces parallel to the plane.
The next two higher energy states are the $e_{g}^{\pi}$ doublet with
wavefunction lobes directed at the faces of the octahedron transverse to the
plane. The highest energy states transform as the $e_{g}^{\sigma}$ doublet in
which the orbital lobes are directed along the octahedral axes, pointing
towards the Cl anions. This behavior is expected because they behave
equivalently to the $e_{g}$ states of the cubic basis. The
$e_{g}^{\sigma}$-$e_{g}^{\pi}$ splitting of approximately $1.5$ eV is somewhat
less than the approximately $2-3$ eV splitting typical in perovskite
transition metal oxides, reflecting the larger energy separation and weaker
hybridization of the metal $d$ states with the halogen ligand $p-$states.

Next, we examine the effects on the electronic structure of varying the
structural parameters (Figure 2). We begin by fixing the TiCl6 octahedra to
$L_{1}=L_{2}=3.430$ Å with all M-X bonds fixed to 2.425 Å and vary the inter-
layer distance. We find that the crystal-field splitting between the
$e_{g}^{\sigma}$ and the $e_{g}^{\pi}/a_{1g}$ complex increases with increased
inter-layer distance as does the corresponding difference in orbital occupancy
(Figure 2A). The $e_{g}^{\pi}-a_{1g}$ splitting is less dependent on the
interlayer distance. For inter-layer distances greater than or equal to the
experimental value measured at ambient pressure, the crystal field splitting
and orbital occupations converge to their asymptotic isolated layer values.

Next, we keep the inter-layer distance fixed to the ambient pressure
equilibrium experimental value of 6.1 Å and vary the $L_{2}/L_{1}$ ratio by
displacing the Cl atoms along $z$. This changes $L_{2}$ while keeping $L_{1}$
the same. We find that increasing $L_{2}$ increases the energy difference
between the highest two and lowest three orbital states states. In contrast to
the effect of increasing the inter-layer distance, changing $L_{2}/L_{1}$ also
decreases the $a_{1g}$ orbital energy relative to the $e_{g}^{\pi}$ energy.
The occupancies change in the corresponding manner. The relative sizes of the
energies and the orbital splittings, however, are not fixed. Importantly, we
find that the critical $L_{2}/L_{1}$ ratio at which the occupations are equal
occurs at a different critical $L_{2}/L_{1}$ ratio that gives a vanishing
energy difference for $e_{g}^{\pi}-a_{1g}$ orbitals. Both correspond to
$L_{2}/L_{1}\neq 0$, underlining that the symmetry of the representation of
the Ti $d$-shell is always at most trigonal, and never cubic.

###  III.2 Orbital Physics in the Trigonal Wannier Basis

Key to understanding the physical effects in a correlated a material is an
appropriate tight-binding model, which can be constructed via a Wannierization
procedure as implemented in standard electronic structure codes. A
Wannierization is in effect a choice of a basis that represents band states in
a certain energy range, along with a projection of the density functional
Hamiltonian onto this basis. It is desirable to choose a basis that is adapted
to the physics of the problem at hand; both because an appropriate choice
provides physical insight and because the form and magnitude of beyond-DFT
interactions depends on the basis chosen.

Standard application of the Wannier90 [51] code to di- and trihalide compounds
tends to lead to the $e_{g}$-$t_{2g}$ basis functions well adapted to $O_{h}$
symmetry. To obtain Wannier functions adapted to trigonal symmetry, we find it
is best to first use Wannier90 to obtain a basis, and then rotate the basis to
obtain a diagonal density matrix. Another option is to diagonalize the Wannier
Hamiltonian at the $\Gamma$-point. A third choice is to diagonalize the on-
site term in the real-space Wannier Hamiltonian. To obtain real-space
isosurfaces, we used the eigenvectors obtained by diagonalizing the density
matrix to build a linear superposition of the Wannier functions obtained
directly from Wannier90.

Figure 3 shows the density of states projected onto different combinations of
the Wannier states obtained for TiCl2 in a standard application of Wannier90
(left panel) and in the basis that diagonalizes the density matrix (center
panel) after the rotation. The total density of states is the same in both
cases. The $e_{g}^{\sigma}$-derived states, which are derived from the same
functions in both cases, are approximately the same. The trigonal functions
allow us to distinguish the higher-lying $a_{1g}$ states from the lower-lying
$e_{g}^{\pi}$ states. Figure 3 also shows the corresponding isosurfaces of the
resulting Wannier functions, which we obtain through linear mixing.

Figure 3: Projected density of states for the Wannier functions for TiCl2 in
the tetragonal basis (left) and trigonal basis (center). Eigenstates of the
density matrix, as obtained by diagonalizing the density matrix (right).
Figure 4:  A) $\vec{R}$ connecting transition metals M, for which we show the
hopping matrices of the Hamiltonian below in B). The diagonal elements are
shown in blue for $e_{g}^{\sigma}$, a1g in green, and $e_{g}^{\pi}$ red.

We now analyze the resulting Hamiltonian $H_{\vec{R}_{i}}$ (Figure 4), where
the rows (columns) correspond to the following $d$ orbitals: the first two
correspond to $e_{g}^{\sigma}$, the middle one to $a_{1g}$ and the last two to
$e_{g}^{\pi}$. For $\vec{R}=(0,0,0)$, we find that the $a_{1g}$ orbital is
orthogonal to the other sets of orbitals in $H_{\vec{R}_{0}}$, similar to what
is obtained in the atomic orbital projection matrix. While within the pairs of
degenerate doublets the mixing is practically 0, there is however mixing
between the two pairs. This $H_{\vec{R}_{0}}$ however possesses the 2-2-1
pattern of eigenvalues consistent with D3d symmetry. Hopping along $z$ between
the Ti atoms $\vec{R}=(0,0,c)$ is primarily driven by $a_{1g}$-$a_{1g}$
hopping, hopping between the two pairs of doublet states, as well as small but
non-zero hopping between equivalent orbitals as gleaned from $H_{\vec{R}_{z}}$
in Figure 4B. Finally, hopping in-plane between the atoms $\vec{R}=(a,0,0)$ is
driven by non-zero hoppings between all orbitals, however the lowest three
orbitals in $H_{\vec{R}_{x}}$ are the largest, as expected from the spatial
orientation of the orbital lobes discussed previously. We also note that one
of the eigenvalues of the hopping Hamiltonian is close to 0. The zero
eigenvalue pertains to hopping among the highest two energy ($e_{g}^{\sigma}$)
orbitals. As a result, there will always be a weakly dispersing flat band
among the higher $e_{g}$ states along high symmetry directions in the zone.

###  III.3 Lower structural symmetries

Some of the dihalides and trihalides exhibit symmetries lower than trigonal.
The point group $C_{2h}$, a subgroup of the trigonal point group $D_{3d}$,
appears in both MX2 and MX3 compounds, namely in ZrI2, as well as in some
polymorphs of RuCl3 and CrI3. This symmetry reduction includes a breaking of
the rotational symmetry, which splits the doublets, and may lead to minor
mixing of the orbital eigenvalues; generally, subgroups of $D_{3d}$ which do
not maintain the 3-fold rotational symmetry lead to a breaking of the
doublets’ degeneracy.

Figure 5: Single layer lattice structure for materials with point group
symmetry C2h, which breaks the degeneracy of the doublets: A) a layer of MX2,
with the structure from a layer of bulk ZrI2, black rectangle represents the
new unit cell B) a layer of MX3, with a similar pattern representing bulk
layers, for example RuCl3 and CrI3 with the AlCl3 structure type; all atoms in
the sketch are within one unit cell. Distances between the atoms are indicated
as follows: the red dotted lines are the longest, followed by purple, blue,
then mustard yellow. The distortions characteristic of ZrI2 are much stronger
than those of the MX3 materials. The M site symmetries are $C_{1h}$ for panel
A and $C_{2}$ for panel B. No d-orbital degeneracies are enforced by either
site symmetry.

Figure 5 shows that the single layers of ZrI2 and RuCl3 or CrI3 no longer have
three-fold rotational symmetry. The in-plane lattice parameters defining the
formula unit of a di- or tri-halide lattice are no longer equivalent, and are
replaced by perpendicular vectors describing the supercell that accomodates
for this lower symmetry as described in the caption for the figure. This
breaks the symmetry of the remaining two orbital doublet pairs. This symmetry
reduction is associated with inequivalence of both the X-X bonds and the M-X
bonds as indicated in Figure 5. In ZrI2 this is plainly discerned as the bond
distortions are sufficiently strong that the structures can be understood to
form 1D zig-zag chains [52]. A minor mixing between the $a_{1g}$ orbital and
the other orbitals appears as a result of the symmetry lowering, and
performing a Wannierization shows us that this minor mixing breaks its
rotational symmetry.

The ZrI2 structure also has additional forms of symmetry breaking beyond that
leading to the $C_{2h}$ point group; for example, the two layers that form the
minimum structural unit are not completely equivalent. The octahedral sizes
differ between the layers. The symmetry of ZrI2 is even lower in certain
experimentally reported structures with two non-equivalent octahedra forming
each layer leading to point group $C_{2v}$. Such further symmetry reductions
are beyond the scope of this paper.

##  IV Electronic Correlations and Local Structure Effects

Many transition metal compounds exhibit correlated insulating states. The
correlation physics that produces these states tends to favor high spin,
filled and empty state configurations; the exact order is determined by a
combination of correlation effects and electron-lattice coupling. Similar to
perovskite oxides, 2D di- and tri-halides can be in the Mott-Hubbard or in the
charge-transfer regime of the Zaanen–Sawatzky–Allen (ZSA) classification [53].
TiCl2 is a clear example of a Mott-Hubbard material as the $d$ and $p$ states
are clearly separated by 4 eV (Figure 6A). Similar to the transition metal
perovskite case, going right along the periodic table, _i.e._ , towards higher
orbital filling, the $p$-$d$ splitting decreases. In Figure 6B, we find the
$p$-$d$ splitting is low for the well-known material NiI2 [54, 55, 56] with
the $p$ and $d$ bands overlapping in energy, corresponding to this insulating
material being in the charge transfer regime.

Figure 6: DFT Band structure and Wannier ‘fat bands’ for a full $p-d$ model
for (left) TiCl2 and (right) NiI2, exemplifying materials that are in the
Mott-Hubbard (TiCl2) and charge-transfer (NiI2) regimes. Blue circles
correspond to transition metal M $d$ state projections, red squares to ligand
$p$ states.  Figure 7: Electronic band dispersions and orbital projected
density of states for ferromagnetic, high-spin TiCl2 at the DFT$+U=2$ eV
level: (A) experimental structure and (B) structure with $L_{2}$ extended to
4.02 Åand inter-layer distance kept constant. Left: Bands, majority spin in
magenta, minority spin in cyan. Right: projected density of states (PDOS),
minority spin rescaled by dividing by 5. Increasing $L_{2}/L_{1}$ closes the
gap by reducing the orbital polarization of the lower three orbitals. (C)
Dependence of the electronic band gap in TiCl2 with $L_{2}/L_{1}$ ratio for a
fixed inter-layer distance.

###  IV.1 Example of TiCl2

We performed $\mathrm{DFT}+U=2$ eV calculations on the dihalide TiCl2 in the
experimental structure and allowing for spin symmetry breaking to gain initial
insight into the role of correlations. The Hubbard on-site Coulomb interaction
allows the opening of a gap in the electronic spectrum ($\approx$ 0.375 eV),
creating a fully orbital and spin polarized state (Figure 7 A). Taking the
TiCl2 experimental structure and allowing for spin-symmetry breaking, we find
a high spin orbitally ordered state. The minority spin channel is completely
unoccupied. Two electrons reside in the majority spin channel of the two
$e_{g}^{\pi}$ orbitals, while the $e_{g}^{\sigma}$ and $a_{1g}$ occupations
are essentially zero.

Calculations with structures obtained by elongating $L_{2}$ while keeping
$L_{1}$ constant closes this gap, while shortening $L_{2}$ further opens the
gap, by shifting the relative energy levels of the lowest three states (Figure
7). Similar effects can be obtained by changing the inter-layer distance as
well. Importantly, we note that in the presence of electronic correlations,
the electronic configuration favors the state that leads to fully occupied and
empty states and the opening of a band gap, and a strong trigonal distortion
is needed to counteract this effect. We note that, in a tetragonal basis, as
discussed before, the lower three $t_{2g}$ states would remain equivalent, as
their lobes would point along equivalent directions, which would obscure the
identity of the active orbitals – the singlet and doublet – participating in
the band gap opening.

###  IV.2 Survey of Broken Symmetry Phases

Figure 8: Different possible electronic states for $d^{1}-d^{9}$ occupations
as discussed in the main text. Configurations that break the doublet
degeneracy would be associated with further lattice distortions, reducing
symmetry below trigonal. States with higher $a_{1g}/e_{g}^{\pi}$ occupation
ratios are likely to have a higher $L_{2}/L_{1}$ ligand bond ratio.

By varying the transition metal ion and considering both dihalides and
trihalides, one can obtain all $d^{n}$ configurations ranging from $d^{1}$
(TiI3) to $d^{9}$ (CuCl2, proposed below). As a result, a range of correlated
insulating states can be produced, as shown in Figure 8. Some of these states
will involve orbital order that breaks the trigonal symmetry, leading to the
appearance of one dimensional lattice structures - or, more generally,
structures that break the three-fold rotational symmetry.

In the $d^{1}$ configuration the natural Mott insulating state involves an
electron in the $a_{1g}$ orbital. This state preserves the trigonal symmetry.
However the multiple superexchange pathways involving also the $e_{g}^{\pi}$
orbitals suggest that unless $L_{2}/L_{1}$ can be made very small the ground
state will be ferromagnetic. However if the $L_{2}/L_{1}$ ratio can be made
large enough, a change in level ordering may occur and it is possible to have
an orbitally ordered Mott insulator with one electron in the $e_{g}^{\pi}$
orbitals and a corresponding orbital order and trigonal symmetry breaking
lattice distortion.

In the $d^{2}$ case, we may generically expect a high spin ground state with
two electrons in the $e_{g}^{\pi}$ oribitals, as found in the DFT$+U$
calculations for TiCl2. For a small enough $L_{2}/L_{1}$ ratio, however, a
state with one electron in the $a_{1g}$ orbital and one in $e_{g}^{\pi}$
doublet may lift the trigonal symmetry.

The relatively smaller ligand field splitting in the halides relative to the
octahedral perovskite oxides suggests that the $d^{4},~{}d^{5},d^{6}$ state
will all be high spin. The $d^{4}$ and $d^{6}$ configurations would also
exhibit orbitally ordered and trigonally broken states; $d^{6}$ can also
support an insulating diamagnetic state with the gap opened by the crystal
field splitting between the 2 higher energy and lower 3 orbitals.

The $d^{7}$ state is most likely to be a $S=1/2$ state with the lowest three
orbitals fully occupied, and one of the $e_{g}^{\sigma}$ orbitals half-filled
and spin-polarized, likely leading to a structurally broken state of the form
found in the materials with $C_{2h}$ point symmetry. The $d^{8}$ configuration
is most likely to be a $S=1$ state with each $e_{g}^{\sigma}$ orbital half-
filled and spin-polarized. In some scenarios, it may be orbitally-polarized
and non-magnetic, associated with a symmetry-breaking mode in the lattice,
most likely as previously discussed.

Such configurations may be susceptible to metal-insulator transitions (MIT),
which we assess using a recently devised machine-learning classification model
[57]. We found that the binary MIT-non MIT classifier tends to predict most of
the 2D halides are candidate MIT compounds, giving a positive MIT
classification for CrCl3, FeCl3, IrBr3, MnCl2, RuCl3, TiCl2, VCl3 and ZrI2,
and a negative classification for CrI3, FeI2 and NiI2. This is likely due to
the similarity of this class of materials to perovskite oxide MIT transition
compounds, both in the structural features exhibited, such as the average
transition metal-ligand distances which range between 2.4 Å and 3.01 Å, the
metal-metal distances of 3.3-4.2 Å, and electronic descriptors, e.g.,
estimated unscreened Hubbard $U$ values. However, due to their intrinsic
broken symmetry and reduced dimensionality, these materials generally will be
exclusively insulating. Similar to MIT compounds, the 2D halides will tend to
display coupled electronic and lattice transitions at low temperatures [30].

Figure 9: A) DFT band structure, and ’fat bands’ of the top two Wannier
orbitals on the hypothetical CuCl2 monolayer discussed in the main text. B)
Real-space isosurfaces of the highest energy d-orbitals, with occupations
0.912 and 0.93 per spin channel. Due to the pseudo-1D nature of the symmetry
broken lattice structure, there is relatively little dispersion along the
$\Gamma$-X direction.

For the $d^{9}$ and possibly $d^{4}$ configurations, the presence of one hole
in the $e_{g}^{\sigma}$ manifold implies a breaking of rotational symmetry and
likely, a one dimensional electronic structure. To examine this possibility,
we perform the following simulation. We take a single layer of the
experimental ZrI2 structure and replace Zr with Cu and I with Cl without
allowing the structure to relax. Owing to the strong charge-transfer character
of the resulting theoretical material, the Cu $d$ manifold is mostly filled.
Nonetheless, after performing the appropriate rotation of the $d$-shell basis,
we find that the two least-occupied orbitals are $e_{g}^{\sigma}$ like (Figure
9). These orbitals display non-zero orbital polarization with one orbital with
$d_{x^{2}-y^{2}}$ character pointing along Cu-Cl bonds along the 1D chain, and
one orbital of $d_{3z^{2}-r^{2}}$ character pointing towards Cu-Cl bonds that
connect the 1D chains. This type of symmetry breaking can likely be exploited,
provided the appropriate material strains, to form 1D conducting chains out of
a 2D structure. This behavior is similar to how the pseudo-2D electronic
structure of cuprate-like high-temperature superconductors forms from out of
3D conducting building blocks. Another orbitally polarized state can also be
obtained with a different orbital basis for the $d^{9}$ configuration [58].

##  V Conclusions

We showed that the highest possible metal orbital symmetry of 2D dihalides and
trihalides comprising edge-shared MX6 octahedra is trigonal, and analysed the
interplay of atomic lattice, orbital physics and correlation effects. Within
this trigonal basis, we showed that the amplitude of the electronic orbital
splittings can be tuned by both the interlayer distance, as well as through
changing the ligand-ligand bond ratio $L_{2}/L_{1}$, and that the effect of
correlations strongly favors ordered states.

In addition, the orbital occupancies can be sensitively tuned through changes
in these atomic structure features, which can be achieved experimentally via
pressure, strain, as well as possibly via optical excitations of the relevant
structural modes to enable control of magnetic configurations and other
electronic ordering. We showed that these materials can be analyzed from the
point of view of the ZSA classification.

We showed how one can build and analyze a Wannier model corresponding to this
reduced symmetry, symmetry which can lead to novel (correlated) electronic
states. Our work serves as a basis to understand correlated phenomena in this
class of materials, and their interplay with lattice symmetry modes, and
easily allows for models to disentangle their roles - similar to recently
built models on perovskites and their Ruddlesden-Popper phases [11, 17, 59,
60]. Our results may be particularly important in the search for candidate
spin-liquis. Spin liquid states are predicted to be found in the vicinity of a
Mott transition and in symmetric highly-frustrated structures. Our generic
finding of strong orbital ordering and associated trigonal symmetry breaking
will be important in appropriately modeling - and discovering - potential spin
liquids and other novel states.

###### Acknowledgements.

This research was supported in part by the National Science Foundation (NSF)
under DMREF Award DMR-1729303. The information, data, or work presented herein
was funded in part by the Advanced Research Projects Agency-Energy (ARPA-E),
U.S. Department of Energy, under Award Number DE-AR0001209.AJM is supported in
part by Programmable Quantum Materials, an Energy Frontier Research Center
funded by the U.S. Department of Energy (DOE), Office of Science, Basic Energy
Sciences (BES), under award DE-SC0019443. The views and opinions of authors
expressed herein do not necessarily state or reflect those of the United
States Government or any agency thereof. The Flatiron Institute is a division
of the Simons Foundation.

## References

  * Schlom  _et al._ [2007] D. G. Schlom, L. Q. Chen, C. B. Eom, K. M. Rabe, S. K. Streiffer, and J. M. Triscone, Annual Review of Materials Research 37, 589 (2007).
  * Garcia  _et al._ [2010] V. Garcia, M. Bibes, L. Bocher, S. Valencia, F. Kronast, A. Crassous, X. Moya, S. Enouz-Vedrenne, A. Gloter, D. Imhoff, C. Deranlot, N. Mathur, S. Fusil, K. Bouzehouane, and A. Barthelemy, Science 327, 1106 (2010).
  * Qi and Rabe [2021] Y. Qi and K. M. Rabe, http://arxiv.org/abs/2103.16466 (2021).
  * Chen  _et al._ [2014] H. Chen, Q. Qiao, M. S. J. Marshall, A. B. Georgescu, A. Gulec, P. J. Phillips, R. F. Klie, F. J. Walker, and C. H. Ahn, Nano. Lett 14, 4965 (2014).
  * Tokura [1998] Y. Tokura, Current Opinion in Solid State and Materials Science 3, 175 (1998).
  * Gibertini  _et al._ [2019] M. Gibertini, M. Koperski, A. F. Morpurgo, and K. S. Novoselov, Nature Nanotechnology 14, 408 (2019).
  * Koohfar  _et al._ [2019] S. Koohfar, A. B. Georgescu, A. N. Penn, J. M. LeBeau, E. Arenholz, and D. P. Kumah, npj Quantum Materials 4, 25 (2019).
  * Molegraaf  _et al._ [2009] H. J. A. Molegraaf, J. Hoffman, C. A. F. Vaz, S. Gariglio, D. Van Der Morel, C. H. Ahn, and J. M. Triscone, Advanced Materials 21, 1 (2009).
  * Mizuguchi  _et al._ [2008] Y. Mizuguchi, F. Tomioka, S. Tsuda, T. Yamaguchi, and Y. Takano, Applied Physics Letters 93, 152505 (2008).
  * Koohfar  _et al._ [2020] S. Koohfar, A. B. Georgescu, I. Hallsteinsen, R. Sachan, M. A. Roldan, E. Arenholz, and D. P. Kumah, Physical Review B 101, 064420 (2020).
  * Georgescu and Millis [2021] A. Georgescu and A. J. Millis, http://arxiv.org/abs/2105.02271 (2021).
  * Guzmán-Verri  _et al._ [2019] G. G. Guzmán-Verri, R. T. Brierley, and P. B. Littlewood, Nature 576, 429 (2019).
  * Shamblin  _et al._ [2018] J. Shamblin, M. Heres, H. Zhou, J. Sangoro, M. Lang, J. Neuefeind, J. A. Alonso, and S. Johnston, Nature Communications 9, 86 (2018).
  * Zubko  _et al._ [2011] P. Zubko, S. Gariglio, M. Gabay, P. Ghosez, and J.-m. Triscone, Annual Review of Condensed Matter Physics , 141 (2011).
  * Forst  _et al._ [2015] M. Forst, A. D. Caviglia, R. Scherwitzl, R. Mankowsky, P. Zubko, V. Khanna, H. Bromberger, S. B. Wilkins, Y. D. Chuang, W. S. Lee, W. F. Schlotter, J. J. Turner, G. L. Dakovski, M. P. Minitti, J. Robinson, S. R. Clark, D. Jaksch, J. M. Triscone, J. P. Hill, S. S. Dhesi, and A. Cavalleri, Nature Materials 14, 883 (2015).
  * Medarde  _et al._ [1998] M. Medarde, P. Lacorre, K. Conder, F. Fauth, and A. Furrer, Physical Review Letters 80, 2397 (1998).
  * Georgescu  _et al._ [2019] A. B. Georgescu, O. E. Peil, A. S. Disa, A. Georges, and A. J. Millis, Proceedings of the National Academy of Sciences of the United States of America 116, 14434 (2019).
  * Domínguez  _et al._ [2020] C. Domínguez, A. B. Georgescu, B. Mundet, Y. Zhang, J. Fowlie, A. Mercy, A. Waelchli, S. Catalano, D. T. Alexander, P. Ghosez, A. Georges, A. J. Millis, M. Gibert, and J. M. Triscone, Nature Materials 19, 1182 (2020).
  * Lee  _et al._ [2018] D. Lee, B. Chung, Y. Shi, N. Campbell, F. Xue, K. Song, J. P. Podkaminer, T. H. Kim, P. J. Ryan, T. R. Paudel, J. W. Spinuzzi, D. A. Tenne, E. Y. Tsymbal, M. S. Rzchowski, L. Q. Chen, J. Lee, and C. B. Eom, Science , 1037 (2018).
  * Caviglia  _et al._ [2013] A. D. Caviglia, M. Först, R. Scherwitzl, V. Khanna, H. Bromberger, R. Mankowsky, R. Singla, Y. D. Chuang, W. S. Lee, O. Krupin, W. F. Schlotter, J. J. Turner, G. L. Dakovski, M. P. Minitti, J. Robinson, V. Scagnoli, S. B. Wilkins, S. A. Cavill, M. Gibert, S. Gariglio, P. Zubko, J. M. Triscone, J. P. Hill, S. S. Dhesi, and A. Cavalleri, Physical Review B - Condensed Matter and Materials Physics 88, 220401(R) (2013).
  * Szymanski  _et al._ [2019] N. J. Szymanski, L. N. Walters, D. Puggioni, and J. M. Rondinelli, Physical Review Letters , 236402 (2019).
  * Chu  _et al._ [2012] J. H. Chu, H. Kue, J. Analytis, and I. Fisher, Science 337, 710 (2012).
  * Lee  _et al._ [2006] P. A. Lee, N. Nagaosa, and X. G. Wen, Reviews of Modern Physics 78, 17 (2006).
  * Norman [2020] M. R. Norman, Physics 13, 10.1103/physics.13.85 (2020).
  * Gariglio  _et al._ [2009] S. Gariglio, N. Reyren, A. D. Caviglia, and J. M. Triscone, Journal of Physics Condensed Matter 21, 164213 (2009).
  * Keimer  _et al._ [2015] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, Nature 518, 179 (2015).
  * Huang  _et al._ [2020] C. Huang, F. Wu, S. Yu, P. Jena, and E. Kan, Physical Chemistry Chemical Physics 22, 512 (2020).
  * Son  _et al._ [2019] S. Son, M. J. Coak, N. Lee, J. Kim, T. Y. Kim, H. Hamidov, H. Cho, C. Liu, D. M. Jarvis, P. A. Brown, J. H. Kim, C. H. Park, D. I. Khomskii, S. S. Saxena, and J. G. Park, Physical Review B 99, 041402(R) (2019).
  * Tian  _et al._ [2019] S. Tian, J. F. Zhang, C. Li, T. Ying, S. Li, X. Zhang, K. Liu, and H. Lei, Journal of the American Chemical Society 141, 5326 (2019).
  * Gati  _et al._ [2019] E. Gati, Y. Inagaki, T. Kong, R. J. Cava, Y. Furukawa, P. C. Canfield, and S. L. Bud’ko, Physical Review B 100, 94408 (2019).
  * Lyu  _et al._ [2020] B. Lyu, Y. Gao, Y. Zhang, L. Wang, X. Wu, Y. Chen, J. Zhang, G. Li, Q. Huang, N. Zhang, Y. Chen, J. Mei, H. Yan, Y. Zhao, L. Huang, and M. Huang, Nano Letters 20, 6024 (2020).
  * An  _et al._ [2019] M. An, Y. Zhang, J. Chen, H. M. Zhang, Y. Guo, and S. Dong, Journal of Physical Chemistry C 123, 30545 (2019).
  * Yadav  _et al._ [2016] R. Yadav, N. A. Bogdanov, V. M. Katukuri, S. Nishimoto, J. Van Den Brink, and L. Hozoi, Scientific Reports 6, 1 (2016).
  * Yokoi  _et al._ [2021] T. Yokoi, S. Ma, Y. Kasahara, S. Kasahara, T. Shibauchi, N. Kurita, H. Tanaka, J. Nasu, Y. Motome, C. Hickey, S. Trebst, and Y. Matsuda, Science 373, 568 (2021).
  * Suzuki  _et al._ [2021] H. Suzuki, H. Liu, J. Bertinshaw, K. Ueda, H. Kim, S. Laha, D. Weber, Z. Yang, L. Wang, H. Takahashi, K. Fürsich, M. Minola, B. V. Lotsch, B. J. Kim, H. Yavaş, M. Daghofer, J. Chaloupka, G. Khaliullin, H. Gretarsson, and B. Keimer, Nature Communications 12, 1 (2021).
  * Li  _et al._ [2021] H. Li, T. T. Zhang, A. Said, G. Fabbris, D. G. Mazzone, J. Q. Yan, D. Mandrus, G. B. Halász, S. Okamoto, S. Murakami, M. P. Dean, H. N. Lee, and H. Miao, Nature Communications 12, 4 (2021).
  * Takagi  _et al._ [2019] H. Takagi, T. Takayama, G. Jackeli, G. Khaliullin, and S. E. Nagler, Nature Reviews Physics 1, 264 (2019).
  * Hermanns  _et al._ [2018] M. Hermanns, I. Kimchi, and J. Knolle, Annual Review of Condensed Matter Physics 9, 17 (2018).
  * Blei  _et al._ [2021] M. Blei, J. L. Lado, Q. Song, D. Dey, O. Erten, V. Pardo, R. Comin, S. Tongay, and A. S. Botana, Applied Physics Reviews 8, 10.1063/5.0025658 (2021).
  * Wang and Sanyal [2021] D. Wang and B. Sanyal, Journal of Physical Chemistry C 125, 18467 (2021).
  * Kong  _et al._ [2019] T. Kong, S. Guo, D. Ni, and R. J. Cava, Physical Review Materials 3, 084419 (2019).
  * Chen  _et al._ [2020] X. Chen, X. Fan, L. Li, N. Zhang, Z. Niu, T. Guo, S. Xu, H. Xu, D. Wang, H. Zhang, A. S. McLeod, Z. Luo, Q. Lu, A. J. Millis, D. N. Basov, M. Liu, and C. Zeng, Nature Physics 16, 631 (2020).
  * Dionne [2009] G. F. Dionne, in  _Magnetic Oxides_ (Springer US, Boston, MA, 2009) pp. 37–106.
  * McGuire [2017] M. A. McGuire, Crystals 7, 10.3390/cryst7050121 (2017).
  * Botana and Norman [2019] A. S. Botana and M. R. Norman, Physical Review Materials 3, 044001 (2019).
  * Hellenbrandt [2004] M. Hellenbrandt, Crystallogr. Rev. 10, 17 (2004).
  * Streltsov and Khomskii [2017] S. V. Streltsov and D. I. Khomskii, Physics-Uspekhi 60, 1121 (2017).
  * Wagner  _et al._ [2019] N. Wagner, R. Seshadri, and J. M. Rondinelli, Phys. Rev. B 100, 064101 (2019).
  * Note [1] Note that for the cubic case, the conventional coordinate system involves a $z$ axis chosen parallel to an M-X bond while for the trigonal case the conventional coordinate system involves a $z$ axis passing through an octahedral face (_i.e._ , along $[111]$ in conventional octahedral coordinates).
  * Note [2] DFT calculations were performed using the Quantum Espresso software package on a 3 atom formula unit, using ultrasoft pseudopotentials using the PBE exchange-correlation functional, an energy cutoff of 748eV and a 20x20x20 k-grid. Orbital occupations were obtained using the orthonormal atomic projections, as implemented in Quantum Espresso; DFT+U calculations were performed with the U applied to orthonormal atomic projections as well, within the Dudarev simplified +U formalism.
  * Mostofi  _et al._ [2014] A. A. Mostofi, J. R. Yates, G. Pizzi, Y. S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Computer Physics Communications 185, 2309 (2014).
  * Guthrie and Corbett [1981] D. H. Guthrie and J. D. Corbett, Journal of Solid State Chemistry 37, 256 (1981).
  * Zaanen  _et al._ [1985] J. Zaanen, G. A. Sawatzky, and J. W. Allen, Physical Review Letters 55, 418 (1985).
  * Ronda  _et al._ [1987] C. R. Ronda, G. J. Arends, and C. Haas, Physical Review B 35, 4038 (1987).
  * Starnberg  _et al._ [1986] H. Starnberg, M. Johnson, and H. Hughes, Journal of Physics C: Solid State Physics 19, 2689 (1986).
  * Van Der Laan  _et al._ [1986] G. Van Der Laan, J. Zaanen, G. A. Sawatzky, R. Karnatak, and J. M. Esteva, Physical Review B 33, 4253 (1986).
  * B. Georgescu  _et al._ [2021] A. B. Georgescu, P. Ren, A. R. Toland, S. Zhang, K. D. Miller, D. W. Apley, E. A. Olivetti, N. Wagner, and J. M. Rondinelli, Chemistry of Materials 33, 5591 (2021).
  * Qin  _et al._ [2021] H. Qin, J. Chen, B. Sun, Y. Tang, Y. Ni, Z. Chen, H. Wang, and Y. Chen, Physical Chemistry Chemical Physics , 22078 (2021).
  * Peil  _et al._ [2019] O. E. Peil, A. Hampel, C. Ederer, and A. Georges, Physical Review B 99, 245127 (2019).
  * Han and Millis [2018] Q. Han and A. Millis, Physical Review Letters 121, 67601 (2018).