# Beyond pairwise network similarity: exploring Mediation and Suppression
between networks

Lucas Lacasa  School of Mathematical Sciences, Queen Mary University of
London, E14NS London (UK)  Institute for Cross-Disciplinary Physics and
Complex Systems (CSIC-UIB), Mallorca (Spain)  Sebastiano Stramaglia
Dipartimento Interateneo di Fisica, Università Degli Studi di Bari, Aldo Moro,
70126 Bari (Italy)  INFN, Sezione di Bari, 70126 Bari (Italy)  Daniele
Marinazzo  Data Analysis Department, Ghent University, 9000 Ghent (Belgium)
IRCCS Ospedale San Camillo s.r.l., 30126 Venice (Italy)

<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Network similarity measures quantify how and when two networks are
symmetrically related, including measures of statistical association such as
pairwise distance or other correlation measures between networks or between
the layers of a multiplex network, but neither can directly unveil whether
there are hidden confounding network factors nor can they estimate when such
correlation is underpinned by a causal relation. In this work we extend this
pairwise conceptual framework to triplets of networks and quantify how and
when a network is related to a second network directly or via the indirect
mediation or interaction with a third network. Accordingly, we develop a
simple and intuitive set-theoretic approach to quantify mediation and
suppression between networks. We validate our theory with synthetic models and
further apply it to triplets of real-world networks, unveiling mediation and
suppression effects which emerge when considering different modes of
interaction in online social networks and different routes of information
processing in the brain.

##  1 Introduction

Networks are usually seen as a parsimonious model to describe the backbone
architecture of complex systems [1]. Accordingly, comparing different systems
boils down to compare their architecture, leading to the notion of network
similarity measure [2, 3, 4, 5, 6]. In graph theory, two graphs are isomorphic
if there exist a vertex permutation which maps one network into the other,
naturally leading to a binary (and not very useful in real-world systems)
notion of similarity. More useful approaches proceed by projecting networks
into a suite of properties summarised in some vector ${\bf p}$ (e.g. degree
distribution, centrality vectors, eigenspectra, etc) and, subsequently
construct a similarity metric $\cal D$ by which two networks $A$ and $B$ are
closer in the space spanned by $\bf p$ if ${\cal D}(A,B)=||{\bf p}_{A}-{\bf
p}_{B}||$ is “small”. Other ideas include the formalisation of graph kernels
[2], comparing networks by comparing the statistics of random walks running
over them [7, 8], or using statistical approaches such as estimating
topological correlations between networks. While in all these approaches we
typically have ${\cal D}(A,B)={\cal D}(B,A)$, i.e. a symmetrical relation, in
many cases this undirected relation is hiding an actual direction (whether
causal or not). As an example, consider social networks. The different layers
of the social network of an individual are typically correlated: my friends
offline tend to be also friends in Facebook. However, such relation is
directional: when a new link –i.e, a new social relationship– is created, then
it is likely that such a link will be replicated within her online social
network too (Facebook, Instagram), but it is proportionally less likely that
the direction of influence is inverted. So the offline and the online social
network of a person are probably similar, but such similarity has a direction.
Furthermore, in many cases such influence is not direct (not causal).
Sometimes, there is a hidden network $\cal C$ that indeed confounds or
mediates the relation between $\cal A$ and $\cal B$. For instance, the
Facebook and Instagram networks of a certain individual are correlated not
because there is a direct, causal relationship between them, but because both
these networks are indeed related to the actual (offline) social network of
the person.

In this work we are interested in understanding and disambiguating when the
relation between two networks $\cal A$ and $\cal B$ (where for instance ${\cal
A}\ r\ {\cal B}$ if ${\cal D}({\cal A},{\cal B})<\epsilon$) is a direct one or
is underpinned by the hidden interaction with a third network $\cal C$. In
particular, $\cal C$ can be independent of $\cal A$ and $\cal B$ (leading to a
direct relation ${\cal A}\ r\ {\cal B}$). $\cal C$ can also act as a hidden
mediator or confounding factor (${\cal A}\ r\ {\cal C}$, ${\cal C}\ r\ {\cal
B}\Rightarrow{\cal A}\ r\ {\cal B})$. Finally, ${\cal C}$ can act as a
suppressor such that $[{\cal A}\oplus{\cal C}]\ r\ {\cal B}$, where $\oplus$
is here to be defined but conceptually means $\cal{A}$ and $\cal C$ interact
synergistically. The terms mediator, confounder and suppressor are inspired by
the information-theoretic framework described in [11].
In what follows we address these questions introducing a set-theoretical
approach where concepts such as network mediation or network suppression
emerge naturally. We benchmark our theory with simple generative models and
then apply it to a range of empirical networks, where we unveil and discuss
the concomitant roles of mediation and suppression.

##  2 Theory

Let $\mathcal{A,B,C}$ be three unweighted networks with adjacency matrices
$\bf A,B,C$, all with the same node set and respective edge sets $a,b$ and $c$
(i.e. they can also be identified with the layers of a multiplex network). Let
us define the network-Jaccard index of two networks
$\text{NJ}({\mathcal{A},B})$ as the Jaccard index over its edge sets

$\text{NJ}({\cal A,B}):=J(a,b)=\frac{|a\cap b|}{|a\cup b|}$ (1)

$\text{NJ}({\mathcal{A},B})$ is a similarity metric, and a distance can be
easily defined as $d({\mathcal{A},B})=1-\text{NJ}({\cal A,B})$. This quantity
alone can be used to initially establish if two networks are related.
Regardless the fact that such relation is effectively undirected or otherwise
is causal (influence), in order to explore whether such relation is
underpinned by a third network $\cal C$ we need to quantify the effect of
conditioning such relation on $\cal C$. Let us then define the partial
network-Jaccard index $\text{NJ}_{p}(\mathcal{A,B}|\mathcal{C})$ of two
networks $\mathcal{A,B}$ conditioned on a third one $\mathcal{C}$ as the
Jaccard index over the edge subsets of $\mathcal{A}$ and $\mathcal{B}$ formed
by those edges which are absent in $\mathcal{C}$:

$\text{NJ}_{p}(\mathcal{A,B}|\mathcal{C})=\frac{|(a\cap b)\setminus
c|}{|(a\cup b)\setminus c|}$ (2)

Let’s see intuitively the effect of conditioning with respect to $\mathcal{C}$
in this way. Suppose initially that $\mathcal{C}$ is totally independent from
$\mathcal{A}$ and $\mathcal{B}$. Then we may expect that the Jaccard index, on
average, will be the same if evaluated just on the links which are absent in
$\mathcal{C}$, so
$\text{NJ}_{p}(\mathcal{A,B}|\mathcal{C})\approx\text{NJ}({\cal A,B})$.
Suppose on the other hand that $\mathcal{A}$ is influencing $\mathcal{B}$
indirectly, with the mediation of $\mathcal{C}$. Then, intuitively, removing
the links of $\mathcal{C}$ would effectively push the partial Jaccard index to
zero. A similar scenario takes place if $\mathcal{A}$ and $\mathcal{B}$ are
undirectedly related through direct relation to a confounding factor
$\mathcal{C}$. Finally, $\mathcal{C}$ could be suppressing the influence of
$\mathcal{A}$ in $\mathcal{B}$. For example, imagine that $\mathcal{B}$
somewhat depends on whether $\mathcal{A}$ and $\mathcal{C}$ interact
synergistically, e.g. if links in $\mathcal{C}$ are more likely to occur if
they are in one network but not on the other (probabilistic XOR gate); then
removing the links of $\mathcal{C}$ will enhance the partial Jaccard index.
To distinguish these three scenarios, we define the Jaccard net difference

$\Delta[\mathcal{A,B;C}]:=\Delta=\text{NJ}_{p}(\mathcal{A,B}|\mathcal{C})-\text{NJ}({\cal
A,B}).$ (3)

Intuitively, if ${\cal C}$ is independent of the relation between ${\cal A}$
and ${\cal B}$ then $\Delta\approx 0$, if it mediates or confounds such
relation then $\Delta<0$, and if it acts as a suppressor then $\Delta>0$.
In what follows we construct simple generative models of independent, mediated
and suppressor interactions, detailed as algorithms, we prove that these
correctly generate these three types of trivariate relations, and depict
numerical simulations of the outcome for finite networks.

Figure 1: Pure models. $\text{NJ}_{p}(\mathcal{A,B}|\mathcal{C})$ vs
$\text{NJ}({\cal A,B})$, calculated on $1000$ realizations of triplets of
networks of $N=50$ nodes wired such that $\mathcal{C}$ plays no effect (green
crosses), plays a mediating effect (violet dots) or a suppression effect (red
crosses) in the relation between $\mathcal{A}$ and $\mathcal{B}$. These
interactions are constructed using the generative models described in
Algorithms 1, 2 and 3 ($p=0.5$ in every case, and $q=1$). For completeness, we
depict the histograms $P(\Delta)$ which certify that these algorithms generate
networks where $\cal C$ play an independent role ($\Delta\approx 0$), a
mediating role ($\Delta<0$) or a suppressing role ($\Delta>0$).

###  2.1 Independency

A simple generative model of independency is given by three independent,
Erdos-Renyi-type models, where in each of the networks each possible link
independently occurs with probability $p$, see Algorithm 1.

Algorithm 1 Uncorrelated()

Output: 3 Erdos-Renyi adjacency matrices $\bf A,B,C$ which are uncorrelated
($\Delta\approx 0)$

1:${\bf A}\leftarrow{\bf 0}$

2:${\bf B}\leftarrow{\bf 0}$

3:${\bf C}\leftarrow{\bf 0}$

4:for $i=1\ {\bf to}\ N$ do

5: for $j=i+1\ {\bf to}\ N$ do

6: if rand $<p$ then $A_{ij},A_{ji}\leftarrow 1$

7: if rand $<p$ then $B_{ij},B_{ji}\leftarrow 1$

8: if rand $<p$ then $C_{ij},C_{ji}\leftarrow 1$

9:return A,B,C

The following theorem can be easily proved.

######  Theorem 2.1.

Let $\bf A$, $\bf B$ and $\bf C$ be as in Algorithm 1. Then
$\mathbb{E}(\Delta)=0$ and the expected values of $\text{NJ}_{p}({\cal
A,B}|{\cal C})$ and $\text{NJ}({\cal A,B})$ are equal to $p/(2-p)$.

The proof of this theorem is given in the appendix. We conclude that, on
average, $\text{NJ}_{p}$ and NJ are nearly equal for uncorrelated networks
generated in this way, as partialization with respect to an independent
network $\mathcal{C}$ does not have any effect. In Fig.1 we illustrate this
case for finite networks with $N=50$ nodes and $p=0.5$, finding that indeed
$\Delta\approx 0$ and that $\text{NJ}_{p}({\cal A,B}|{\cal
C})\approx\text{NJ}({\cal A,B})\approx 1/3$, in good agreement with the
theorem.

Figure 2: Adjacency matrices cartoons displaying mediation and suppression.
(Top) Network $\cal C$ is mediating the relation between $\cal A$ and $\cal B$
(Algorithm 2). (Bottom) Network $\cal C$ acts as a suppressor between $\cal A$
and $\cal B$ (Algorithm 3).

###  2.2 Mediation

Suppose now that $\mathcal{A}$ and $\mathcal{B}$ are both dependent on
$\mathcal{C}$, i.e. where there is a link in $\mathcal{C}$, then there is a
link in $\mathcal{A}$ and $\mathcal{B}$ (see Fig.2 for an illustration of such
case, and Algorithm 2 for a formal recipe of this generative model).

Algorithm 2 Mediated()

Output: 3 Erdos-Renyi adjacency matrices $\bf A,B,C$ where $\bf C$ mediates
relation between $\bf A$ and $\bf B$ ($\Delta>0)$

1:${\bf A}\leftarrow{\bf 0}$

2:${\bf B}\leftarrow{\bf 0}$

3:${\bf C}\leftarrow{\bf 0}$

4:for $i=1\ {\bf to}\ N$ do

5: for $j=i+1\ {\bf to}\ N$ do

6: if rand $<p$ then $A_{ij},A_{ji}\leftarrow 1$

7: if rand $<p$ then $B_{ij},B_{ji}\leftarrow 1$

8: if rand $<p$ then $C_{ij},C_{ji},A_{ij},A_{ji},B_{ij},B_{ji}\leftarrow 1$

9:return A,B,C

This describes a situation where $\mathcal{C}$ mediates the relation between
$\mathcal{A}$ and $\mathcal{B}$ (or, alternatively, $\mathcal{C}$ is
confounding that relation). Partializing with respect to $\cal C$ removes the
dependence between ${\cal A}$ and ${\cal B}$ due to $\cal C$, which
intuitively leads to $\Delta<0$. The following theorem can be proved:

######  Theorem 2.2.

Let $\bf A$, $\bf B$ and $\bf C$ be as in Algorithm 2. If $\bf A$ and $\bf B$
share at least one edge besides the common edges shared with $\bf C$, then
$\Delta<0$.

The proof of this theorem is also put in an appendix. In Fig.1 we show
numerical results for finite networks with $N=50$, with $p=0.5$.

###  2.3 Suppression

Finally, let us consider the case where $\cal B$ depends on the interaction of
$\cal A$ and $\cal C$ such that, an edge occurs in $\cal B$ with a certain
probability if it appears in $\cal A$ but not in $\cal C$ or alternatively if
it appears in $\cal C$ but not in $\cal A$ (see Fig. 2 for an illustration).
This is akin to a probabilistic XOR gate. Then on average
$\text{NJ}_{p}(\mathcal{A,B}|\mathcal{C})>\text{NJ}(\mathcal{A,B})$, i.e.
partializing with respect to $\cal C$ in this case evidences suppression
effects.

Algorithm 3 Suppression()

Output: 3 Erdos-Renyi adjacency matrices $\bf A,B,C$ where $\bf C$ acts as a
suppressor between $\bf A$ and $\bf B$ ($\Delta>0)$

1:${\bf A}\leftarrow{\bf 0}$

2:${\bf B}\leftarrow{\bf 0}$

3:${\bf C}\leftarrow{\bf 0}$

4:for $i=1\ {\bf to}\ N$ do

5: for $j=i+1\ {\bf to}\ N$ do

6: if rand $<p$ then $A_{ij},A_{ji}\leftarrow 1$

7: if rand $<p$ then $B_{ij},B_{ji}\leftarrow 1$

8: if rand $<p$ then $C_{ij},C_{ji}\leftarrow 1$

9: if $A_{ij}*C_{ij}=0$ & $A_{ij}+C_{ij}>0$ & rand $<q$ then
$B_{ij},B_{ji}\leftarrow 1$

10:return A,B,C

Algorithm 3 encapsulates the generative model. The following theorem can be
proved:

######  Theorem 2.3.

Let $\bf A$, $\bf B$ and $\bf C$ be as in Algorithm 3. Then
$\mathbb{E}(\Delta)>0$.

The proof for this theorem is also put in an appendix. In Fig.1 we show
numerical results for finite networks with $N=50$, with $p=0.5$ and $q=1$,
which are in full agreement with the theorem.

Figure 3: (Left panel) Values of $\Delta$ obtained averaging 50 realizations
of an interpolating model giving a blend of mediation and suppression (see the
text), where each network has $N=300$ nodes and varying model parameter $p$,
as a function of the interpolation parameter $\mu$ (in each of the N(N-1)
steps we apply Algorithm 2 with probability $1-\mu$ and Algorithm 3 with
$\mu$). (Right panel) Revised version of the left panel, this time selectively
removing all suppression and mediation effects according to our theoretical
framework. The original curve, for an interpolating model generating a blend
of mediation and suppression, is depicted in black. The blue curve is a pure
randomisation, which generates $\Delta\approx 0$. The dotted green line
corresponds to a selective rewiring that removes all hidden suppression: in
that case the curve is kept always below zero (increasing $\mu$ increases the
amount of Algorithm 3, but then is selectively rewired, thus effectively
randomising the networks and pushing $\Delta$ to zero). The pink curve is the
result of a selectively rewiring that removes all hidden mediation: in that
case the curve is pushed to the regime $\Delta>0$. As $\mu$ increases, the
amount of Algorithm 3 (generating suppression) is increased, hence pushing
$\Delta$ to larger values. The dashed yellow and pink lines are the result of
selectively rewiring on the randomised networks, and only highlight the
residual values of suppression or mediation which occur by chance (as a finite
size effect) in randomised networks.

###  2.4 Coexistence of mediation and suppression effects

When we abandon ideal cases where only suppression or only mediation are
present, and we go towards a mixture of the two, it becomes evident that both
effects can be hidden and a single $\Delta$ cannot in principle tell us if the
system evidences only one out of two mechanisms. To investigate coexistence of
both mechanisms, we run a simulation in which algorithms 2 and 3 above are
combined: in each step with probability $1-\mu$ we take algorithm 2 and with
probability $\mu$ we take algorithm 3. The resulting model linearly
interpolates mediation and suppression, such that measurable
$\Delta=(1-\mu)\Delta_{\text{med}}+\mu\Delta_{\text{syn}}$, where
$\Delta_{\text{med}}$ and $\Delta_{\text{syn}}$ are hidden. The results are
depicted in the left panel of Fig.3, for different instances of parameter $p$
and $q=1$. We can have negative, null or positive values of $\Delta$
underpinned by a balance of both mediation and suppression mechanisms, and
actually for the concrete set of parameters, the effect of mediation in
$\Delta$ is slightly stronger than the effect of suppression (this unbalance
gets more pronounced for $q<1$). This simple interpolating model thus leads us
to conclude that, in real cases, we might for instance be naively measuring
$\Delta<0$ and misleadingly concluding that there is only mediation where in
fact both mediation and suppression could be at play. Accordingly, a measure
describing the effects of suppression and mediation is not enough to describe
and resolve the simultaneous presence of both.

Algorithm 4 Null()

Input: 3 adjacency matrices A, B, C and mode $X$ (mediation (M) or suppression
(S))

Output: Null model net difference $\Delta_{W,X}$

1:${\bf B}2\leftarrow\textsc{FullRewire}({\bf A,B,C})$

2:$\Delta_{X}\leftarrow\textsc{SelectiveRewire}({\bf A,B,C})$

3:$\Delta_{RX}\leftarrow\textsc{SelectiveRewire}({\bf A,B2,C})$

4:$\Delta_{W,X}\leftarrow[\Delta_{X}-\Delta_{RX}]/X_{\text{max}}$

5:return ${\Delta_{W,X}}$

1:function FullRewire($\bf G$)

2: for $\ell_{ij}\in{\bf G}$ do

3: if rand $<p$ then $\ell_{ij}\leftarrow\ell_{kl}$,
$\ell_{kl}\leftarrow\ell_{ij}$

4: return G

1:function SelectiveRewire(A, G, C, X)

2: if $X=S$ then

3: for $\ell^{\prime}_{ij}\in{\bf G}^{\prime}$ do

4: if $\ell^{\prime}_{ij}=1\ \&\ A_{ij}=1\ \&\ C_{ij}\neq 1$ then
$\ell^{\prime}_{kl}\leftarrow\ell^{\prime}_{ij};\ \ell^{\prime}_{ij}\leftarrow
0$

5: if $X=M$ then

6: for $\ell^{\prime}_{ij}\in{\bf G}^{\prime}$ do

7: if $\ell^{\prime}_{ij}=1\ \&\ A_{ij}=1\ \&\ C_{ij}=1$ then
$\ell^{\prime}_{kl}\leftarrow\ell^{\prime}_{ij};\ \ell^{\prime}_{ij}\leftarrow
0$

8: $\textsc{selective-rewire}\leftarrow\text{NJ}_{p}({\bf A,G}^{\prime}|{\bf
C})-\text{NJ}({\bf A,G}^{\prime})$

9: return SelectiveRewire

In order to disentangle both effects we now introduce Algorithm 4, which
applies both for constructing null models for mediation (M) and suppression
(S). To construct a surrogate where all suppression has been removed, starting
from $\cal A,B,C$, we perform a selective rewiring in $\cal B$, where only
those links in $\cal B$ which are also present in $\cal A$ but not in $\cal C$
(or that also appear in $\cal C$ but not in $\cal A$) are rewired randomly.
Similarly, to construct a surrogate where all mediation has been removed,
starting from $\cal A,B,C$, we perform a selective rewiring in $\cal B$, where
only those links in $\cal B$ which are also present in $\cal A$ and in $\cal
C$ are rewired randomly.
We then compute again the net Jaccard difference on the rewired versions,
which are are labelled $\Delta_{S}$ (applied to the case where suppression is
removed) and $\Delta_{M}$ (applied to the case where mediation is removed)
respectively. The heuristic is then simple: if there is e.g. hidden
suppression in the data (respectively mediation), then $\Delta_{S}<\Delta$
(respectively $\Delta_{M}>\Delta$), whereas if such mechanism is absent then
$\Delta_{S}\approx\Delta$ ($\Delta_{M}\approx\Delta$).
Now, we also need to take into account finite-size effects which irremediably
add spurious mediation and suppression effects (i.e. triplets of purely
random, uncorrelated networks will show small but non-zero mediation and
suppression due to chance). To counterbalance such effects, we also proceed to
selectively remove suppression and mediation from a completely randomly
rewired version of $\cal B$, which we call $\bf B2$, leading to two new
indices: $\Delta_{RS}$ and $\Delta_{RM}$. We can finally combine these to
produce normalised indices of mediation ($\bar{m}$) and suppression
($\bar{s}$) by normalising them dividing over the maximum possible value of
suppression (mediation) attainable by a generative model such as algorithm $3$
($2$) to the triplet of networks, i.e.

${\bar{m}}=\frac{\Delta_{S}-\Delta_{\text{R,S}}}{M_{max}},\ \
{\bar{s}}=\frac{\Delta_{M}-\Delta_{\text{R,M}}}{S_{max}}.$

Accordingly, the role that $\cal C$ plays in the relation of $\cal A$ and
$\cal B$ is described by the duple $(\bar{m},\bar{s})$.
Additionally, a significance value for these indices could be defined as:

$\sigma_{\text{\\{S,M\\}}}=\left\lvert\frac{\Delta_{\\{S,M\\}}-\Delta_{\text{R,\\{S,M\\}}}}{\text{std}(\Delta_{\text{R,\\{S,M\\}}})}\right\rvert.$

It is worth to stress that for large networks the finite size effects become
less common, and $\Delta_{\text{R,\\{S,M\\}}}$ will tend to zero. In the same
way any value of suppression and mediation will be highly significant.

Figure 4: Illustration of Algorithm 4.  We illustrate the procedure of how to
disentangle mediation and suppression in two examples. (a) Selective rewiring
applied on a triplet of independent networks generated by Algorithm 1, with
$N=300$ nodes and $p=0.3$. The original value of the net Jaccard difference
(close to zero) is denoted $\Delta_{0}$. Each selective rewiring is repeated
500 times, and histograms associated to each process are built. The selective
removal of suppression and mediation yields only a marginal change in
$\Delta$, similar to the one performed on a randomised version, what indicates
that the amount of suppression and mediation in this configuration is residual
and only due to finite size effects, as expected since $\cal C$ is independent
of $\cal A,B$ by construction. (b) Similar to panel (a) but applied on the C
Elegans triplet, where network $\cal C$ is assigned to the electrical layer.
The actual value $\Delta_{0}<0$, initially suggesting mediation. The selective
removal of suppression (mediation) significantly push the histograms towards
more negative (positive) values of $\Delta$ –much more than such selective
removal performed on a randomisation– suggesting that there exist a
significant amount of mediation and suppression. All histograms are built from
$500$ independent rewiring realisations.

For illustration, in the left panel of Fig.4 we show the effects of the
sequence of selective rewirings on $\Delta$ applied to a particular example of
three independent, Erdos-Renyi networks with $N=300$ nodes and wiring
probability $p=0.3$ (i.e., Algorithm 1). The original value of $\Delta$ is
very close to zero, as well as the ones obtained from a full randomisation of
$\cal B$. Since in this example the networks are independent, any mediation or
suppression is only a spurious residual due to finite size effects, thus this
residual is flagged out in similar terms by a selective rewiring on the actual
network $\cal B$ ($\Delta_{X}$) or on its full randomisation ($\Delta_{RX}$),
hence the violet and green histograms overlap, and similarly the pink and pale
blue ones also overlap. As an additional illustration, we applied the sequence
of selective rewirings on the results of the interpolating model. Results are
shown in the right panel of Fig.3.

##  3 Empirical networks

We now turn to real-world networks and consider four types of 3-layer
multiplex networks, including (i) different modes of social interaction in
Twitter during the 2014’s New York City Climate March (NYC), (ii) different
types of social interaction –proximity, phone call/text message, Facebook– as
collected in Denmark (Copenhagen), (iii) different interpersonal relations
inside a corporation (Lazega law firm) and (iv) different synaptic junctions
in a brain network (C Elegans), see appendix for details.
To begin with, in Fig.5 we confirm that, with the exception of the pair NYC
Retweets vs Replies, all other possible pairs of layers in the four examples
are indeed genuinely related –i.e., showing substantially more similarity than
a null model–. In each case we plot $\text{NJ}({\cal A,B})$ (blue bars) and as
a reference, as black lines we also plot the average result of
$\text{NJ}({\cal A,B})$ ($\pm$ one standard deviation) after $\cal A,B$ have
been appropriately randomised. We confirm that the similarity between each
pair of networks is not the result of a finite-size effect and thus exploring
the role of a third network ($\cal C$) is justified.
We then turn to analyse the role of $\cal C$. For illustration, the whole
selective rewiring process described in Algorithm 4 is depicted with detail
for a specific example (the case of the C Elegans multiplex where we explore
the role that the electrical synapses layer play in the relation between the
monadic and the polyadic layers) in the right panel of figure 4. We provide
the original value $\Delta_{0}$, and the distributions of the $\Delta$ values
obtained after each of the rewiring procedure, concluding that this network
indeed shows non-negligible mediation and suppression effects.

Figure 5: Similarity between pairs of real-world networks.  Values of
NJ$({\cal A,B})$ computed on the four empirical multiplex networks considered
(for each multiplex, we consider the three pair permutations). As a reference,
black horizontal lines display NJ${}_{\text{null}}({\cal A,B})$ which computes
the average over several randomisations of layers $\cal A,B$ (red lines
correspond to $\pm$ one standard deviation. We conclude that all pair of
layers are genuinely related with the possible exception of the pair NYC
replies vs retweet.

The normalized indices of mediation and suppression for the rest of
permutations in all the real-world multiplex networks are reported in figure
6. The first thing we can observe is that overall there is substantially more
mediation than suppression, although we also observe the latter mechanism. All
effects are statistically significant ($\sigma\gg 1$ in every case, data not
shown) except for the suppression in the Lazega advice layer and the NYC
Retweet layer, where $\sigma\approx 2$ in both cases.
In the case of the Copenhagen multiplex, the only layer which evidences a
significant role in the relation of the other two is the phone/sms layer,
which we show displays both mediation and suppression effects, although the
former is notably stronger. For the proximity network we considered averaged
values over the whole four weeks, and used an adjacency matrix of a density
comparable to the one of Facebook links, corresponding to the closest
proximity range. Also the phone network was built irrespective of the timing
of the interaction.
In the case of the Lazega law firm, all three layers display very high
mediation, but such effect is notably stronger for the co-working network,
i.e. within this firm dyadic friendships are related to the dyadic advisory
relations, and this is mediated by the fact that these are co-working. Only
the friendship layer displays a suppressor effect (the one played by the
advice layer is non-significant), i.e. pairs of individuals that are not co-
working can have an advisory relationship (or otherwise) because they are
friends, or pairs of co-working will also have an advisory relationship
without the needs of them being friends.
In the case of the Twitter triplet (NYC), only the Replies network shows a
mediating effect.
Finally, in the nervous system multiplex (C Elegans) we can see that all
layers display some amount of mediation and suppression. The electrical
synapses layer is the one displaying a stronger suppressing effect, whereas
the monadic chemical layer is the one that displays a larger mediation role.
The increased suppression role of the layer of electrical synapses reflects
the evidence that chemical and electrical synapses closely interact and serve
related functions [14], so that when the either of the chemical layers is
taken as $\cal{C}$ the presence of the other chemical layer accounts for a
reduced suppression/mediation.

Figure 6: Normalised indices of mediation and suppression for four real-world
multiplex networks. For each multiplex we permute the role of $\cal C$ across
the three layers. Dot clouds are the result of repeating the rewiring
procedures 500 realisations (for multiplexes with a large number of nodes, the
figure displays very little dispersion, and the cloud is only perceptible for
the Lazega triplet which is indeed the smallest multiplex). Figure 7:
Modulating mediation.  (Outer panel) $(\bar{m},\bar{s})$-plane for the role of
the proximity network with respect to the relation of the Facebook and the
phone/sms network (Copenhagen multiplex), for different proximity networks
reconstructed by changing the spatial scale over which links are defined. The
network displays a very clear mediating effect only when links represent close
physical proximity (yellow), for other spatial scales such effect is lost.
(Inset panel) Link average weight distribution (Received Signal Strength
(RSSI)) is depicted in black. According to this distribution, we build seven
non-overlapping RSSI windows. RSSI is inversely related to distance (the
larger the signal strength, the closer two nodes are) so these windows
represent different spatial scales. The right panel describes
$\bar{m}+\bar{s}$ for each spatial scale, further emphasising that only in the
closer spatial scale the proximity network is playing a strong mediating role.

As a final analysis, and in order to show how suppression and mediation can be
functionally modulated within a particular real-world example, we examine the
role played by the Proximity layer in the relation between the Facebook and
phone calls/sms layers when such layer is systematically varied. In this
multiplex, the proximity network is originally reconstructed using Bluethooth
signal strength between participants, by assigning a link between each pair of
nodes whose relative Bluetooth strength, averaged over the whole period of the
recording (four weeks) belongs to a given range. In order to build different
proximity networks (each of them accounting for a different spatial scale) and
at the same time keeping the edge density constant, we build non-overlapping
Bluetooth intensity ranges by taking into account the original Bluetooth
intensity distribution (see the inset of Fig.7). In this way, ranges are non-
uniform but the number of edges in each of this range is the same, hence the
resulting proximity networks have all the same edge density while describing
different scales of physical proximity. The intuition is that only the smaller
scale is a meaningful proximity network, and for larger spatial scales the
resulting proximity networks do not really imply any real interaction between
the nodes. Then, for each resulting proximity network, we estimate the role it
plays in the relation between the other two networks and plot it in the
$(\bar{m},\bar{s})$-plane. Results are shown in the outer panel of Fig.7. For
proximity networks describing large spatial scales, the network is essentially
independent of the other two, and it is only when the proximity network
captures smaller spatial scale (i.e. when links describe real physical
proximity) that an indirect effect (notably, mediation) gets amplified.

##  4 Discussion

In this paper we have proposed a simple strategy to assess the role that a
given network might play in shaping the relation between other two networks,
thus enlarging the paradigm of network similarity beyond the classical
pairwise comparison. This approach is aligned to a recent endeavour that aims
at going beyond dyadic interactions in the characterisation of complex systems
[10], and takes inspiration from the causal mediation literature [11, 12]. We
make use of a set-theoretic approach to define a similarity metric between a
pair of networks and to further explore if such relation is independent of,
mediated or suppressed by a third network which might be hidden. We introduce
simple generative models that, we prove, produce pure mediation and
suppression. We then explore the coexistence between mediation and suppression
and develop a procedure to disentangle both indirect effects. The whole
methodology is subsequently applied to a range of real-world, 3-layer
multiplex networks, and we unveil previously unnoticed mediation and
suppression effects in social and brain networks.
We hope this work sparks further research in several areas. First, the
simplicity and tractability of the approach makes it easily applicable across
the disciplines. Second, our approach can be readily generalised to consider
not just isolated triplets of networks. Indeed, one can sequentially apply
this protocol to a multiplex network of arbitrary number of layers, or to a
temporal network, and accordingly derive concepts of causality and
directionality in this context.

## Acknowledgements

LL acknowledges funding from EPSRC EP/P01660X/1. SS was supported by MIUR
project PRIN 2017WZFTZP “Stochastic forecasting in complex systems”. DM was
supported by the Belgian Embassy in the United Kingdom through the Belgian
Chair at the University of London and by the Flemish Fund for Research (FWO)
through a sabbatical bench fee.

## Appendix: Proofs of theorems

### Proof of theorem 1

###### Proof.

Since $\bf A$, $\bf B$ and $\bf C$ are ER(p), all links are independent. We
then use expected values, which should be representative as $N$ is large, to
directly compute the expected size of the different sets in eq. 1. With a
total of $\ell:=L(L-1)/2$ independent trials, it is easy to see that the
expected sizes are

$\mathbb{E}(|a|)=\mathbb{E}(|b|)=2\ell p,\ \mathbb{E}(|a\cap b|)=2\ell p^{2}.$

Since $|a\cup b|=|a|+|b|-|a\cap b|$ we also have

$\mathbb{E}(|a\cup b|)=2\ell(2p-p^{2}),$

such that

$\mathbb{E}(J(A,B))=\frac{p}{2-p}$

Consider now $J_{p}(A,B|C)$. On the one hand, we have $(a\cap b)\setminus
c=(a\cap b)\cap c^{\prime}$, where $c^{\prime}$ is the complement of $c$.
Since all networks are independent,

$\mathbb{E}(|(a\cap b)\setminus c|)=2\ell p^{2}(1-p)$

On the other hand, we have

$\mathbb{E}(|(a\cup b)\setminus c|)=2\ell\cdot(1-p)\cdot\underbrace{Prob(a\cup
b)}_{p(2-p)}=2\ell\cdot(1-p)\cdot(2p-p^{2}),$

so that altogether
$\mathbb{E}(J_{p}(A,B|C))=\mathbb{E}(J(A,B))=\frac{p}{2-p}$, and
$\mathbb{E}(\Delta)=\mathbb{E}(J_{p}(A,B|C))-\mathbb{E}(J(A,B))=0$. ∎

### Proof of theorem 2

###### Proof.

In this case the proof uses basic arguments of set theory. We aim to prove
that $J_{p}(A,B|C)<J(A,B)$, i.e.

$\frac{|(a\cap b)\setminus c|}{|(a\cup b)\setminus c|}<\frac{|a\cap b|}{|a\cup
b|}$

Let us define the residual sets $R_{a}=a\setminus c$, $R_{b}:=b\setminus c$,
and the residuals’ intersection $R_{i}=R_{a}\cap R_{b}$ and union
$R_{u}=R_{a}\cup R_{b}$. $R_{a}$ and $R_{b}$ are clearly non-empty according
to Algorithm 2. $R_{i}$ is not guaranteed to be non-empty as $A$ and $B$ are
probabilistic, so we need to assume in what follows that $R_{i}\neq\O$.
$R_{u}$ is trivially non-empty according to Algorithm 2.
Since for any three sets $x$, $y$ and $z$ we have that $(x\cap y)\setminus
z=(x\setminus z)\cap(y\setminus z)$ and $(x\cup y)\setminus z=(x\setminus
z)\cup(y\setminus z)$, we also have that $|(a\cap b)\setminus c|=|R_{i}|$ and
$|(a\cup b)\setminus c|=|R_{u}|$. Incidentally, our previous assumption can
now be easily interpreted: we need to assume that networks A and B share at
least one edge besides the common edges shared with C, that’s why this
assumption is indeed stated in the theorem.
Now, since in this case $c\subset(a\cap b)$ it is easy to see that $a\cap
b=c\cup R_{i}$ and $a\cup b=c\cup R_{u}$. Also by construction we have $c\cap
R_{i}=\O$ and therefore $|a\cap b|=|c|+|R_{i}|$. Similarly, also by
construction $c\cap R_{u}=\O$ and therefore $|a\cup b|=|c\cup
R_{u}|=|c|+|R_{u}|$. Therefore we aim to prove that

$\frac{|R_{i}|}{|R_{u}|}<\frac{|c|+|R_{i}|}{|c|+|R_{u}|}$

Since $R_{i}$ and $R_{u}$ are respectively the intersection and the union of
two sets, it trivially follows that $|R_{i}|\leq|R_{u}|$, where inequality
only saturates when $a=b$, and otherwise is strict. Let us indeed assume
$A\neq B$, thus enforcing the strict inequality. Rearranging terms:

$\displaystyle|R_{i}|<|R_{u}|$ $\displaystyle\iff$
$\displaystyle\frac{1}{|R_{i}|}>\frac{1}{|R_{u}|}$ $\displaystyle\iff$
$\displaystyle\frac{|c|}{|R_{i}|}>\frac{|c|}{|R_{u}|},\ (|c|>0)$
$\displaystyle\iff$ $\displaystyle
1+\frac{|c|}{|R_{i}|}>1+\frac{|c|}{|R_{u}|}$ $\displaystyle\iff$
$\displaystyle\frac{|R_{i}|+|c|}{|R_{i}|}>\frac{|R_{u}|+|c|}{|R_{u}|}$
$\displaystyle\iff$
$\displaystyle\frac{|R_{i}|}{|R_{u}|}<\frac{|c|+|R_{i}|}{|c|+|R_{u}|}$ (4)

∎

### Proof of theorem 3

###### Proof.

The first thing to observe is that in Algorithm 3 the condition
$c\subset(a\cap b)$ is not met, therefore theorem 2.2 does not hold in this
case. Let us now define the following sets:
Let $a_{p}\subset a$ be the subset of edges in $a$ such that $a_{p}\cap c=\O$.
The elements of this subset will be in $b$ with probability $q$, hence on
average $|a_{p}|q$ edges of $a_{p}$ will be in $b$.
Let $c_{p}\subset b$ be the subset of edges in $c$ such that $c_{p}\cap a=\O$.
The elements of this subset will be in $b$ with probability $q$, hence on
average $|c_{p}|q$ edges of $c_{p}$ will be in $b$.
Let $r\subset b$ be the subset of edges in $b$ which are neither in $a$ nor in
$c$, i.e. $r\cap(a\cup c)=\O$.
By symmetry, we have $\mathbb{E}(|a_{p}|)=\mathbb{E}(|c_{p}|)$. According to
Algorithm 3, we have

$\mathbb{E}(|(a\cap b)\setminus c|)=\mathbb{E}(|(a\cap
b))=\mathbb{E}(|a_{p}|)q,$ $\mathbb{E}(|(a\cup b)\setminus
c|)=\mathbb{E}(|a_{p}|)+\mathbb{E}(|r|),$ $\mathbb{E}(|(a\cup
b))=\mathbb{E}(|a_{p}|)+\mathbb{E}(|c_{p}|)q+\mathbb{E}(|r|),$

and thus

$\mathbb{E}(J_{p}(A,B|C))=\frac{\mathbb{E}(|a_{p}|)q}{\mathbb{E}(|a_{p}|)+\mathbb{E}(|r|)},$

whereas

$\mathbb{E}(J(A,B))=\frac{\mathbb{E}(|a_{p}|)q}{\mathbb{E}(|a_{p}|)+\mathbb{E}(|c_{p}|)q+\mathbb{E}(|r|)}.$

Therefore as long as $\mathbb{E}(|c_{p}|)q>0$, we have
$\mathbb{E}(J_{p}(A,B|C))>\mathbb{E}(J(A,B))$, yielding
$\mathbb{E}(\Delta)>0$. Since A and C are independent, $c_{p}$ is not empty
with large probability, so given a value of $q$, for sufficiently large $N$
this condition is always met. ∎

## Appendix: empirical networks

Summary of empirical networks
---
Triplet |  Type |  $N$ |  Network $\\#1$ |  Network $\\#2$ |  Network $\\#3$
C. Elegans |  brain |  279 |  monadic (1639 edges) |  polyadic (3193 edges) |  electrical (1031 edges)
NYC |  Twitter |  102439 |  retweet (213754 edges) |  mentions (131679 edges) |  replies (8062 edges)
Lazega law firm |  social (offline) |  71 |  cowork (892 edges) |  friendship (575 edges) |  advice (1104 edges)
Copenhagen |  social (offline/online) |  751 |  proximity (13020 edges) |  facebook (12847 edges) |  calls/sms (1760 edges)
Table 1: Summary of network specificities. The first three examples are
multiplex networks collected from comunelab.fbk.eu/data.php, whereas the
fourth one is collected from icon.colorado.edu. The C Elegans multiplex
describes the Caenorhabditis elegans connectome, where layers correspond to
different synaptic junctions: chemical monadic (”MonoSyn”), polyadic
(”PolySyn”) and electrical (”ElectrJ”) [13, 15]. The NYC multiplex describes
Twitter activity during an exceptional event, the NYC Climate March in 2014,
and layers correspond to retweet, mentions and replies [18]. The Lazega law
firm depicts three kinds of social relationships (Co-work, Friendship and
advice) between partners and associates of a corporate law partnership [16,
17]. Finally the Copenhagen multiplex describes social interaction in three
layers corresponding to phone calls and text messages (merged), Facebook
friendships, and proximity as measured with strength of Bluetooth signal [19].

## References

  * [1] V. Latora, V. Nicosia, G. Russo, Complex Networks: principles and applications (Cambridge University Press, 2017).
  * [2] S.V.N. Vishwanathan, N.N. Schraudolph, R. Kondor, R. and K.M. Borgwardt. Graph kernels. The Journal of Machine Learning Research, 11 (2010) pp.1201-1242.
  * [3] T.A. Schieber, L. Carpi, A. Diaz-Guilera, P.M. Pardalos, C. Masoller, and M.G. Ravetti, Quantification of network structural dissimilarities. Nature Communications 8 (1) (2017).
  * [4] M. De Domenico, and J. Biamonte, Spectral entropies as information-theoretic tools for complex network comparison. Physical Review X, 6(4), p.041062 (2016).
  * [5] A. Santoro, V. Nicosia, ”Algorithmic complexity of multiplex networks”, Phys. Rev. X 10, 021069 (2020).
  * [6] M. De Domenico, V. Nicosia, A. Arenas, V. Latora, Structural reducibility of multilayer networks, Nature Communications 6, 6864 (2015)
  * [7] V. Nicosia, M. De Domenico, V. Latora. Characteristic exponents of complex networks, EPL 106, 58005 (2014).
  * [8] L. Lacasa, J. Gomez-Gardeñes, Correlation dimension of complex networks, Physical Review Letters (2013).
  * [9] M. Lin, H.C. Lucas Jr, and G. Shmueli. Research Commentary—Too Big to Fail: Large Samples and the p-Value Problem. Information Systems Research 24:4, 906-917 (2017).
  * [10] F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J-G Young, and G. Petri. Networks beyond pairwise interactions: structure and dynamics. Physics Reports (2020).
  * [11] D.P. MacKinnon, J.L. Krull, and C.M. Lockwood. Equivalence of the Mediation, Confounding and Suppression Effect, Prev Sci. (2000).
  * [12] J.T. Lizier, N. Bertschinger,J. Jost, and M. Wibral. Information Decomposition of Target Effects from Multi-Source Interactions: Perspectives on Previous, Current and Future Work Entropy 20(4), 307 (2018).
  * [13] B.L. Chen, D.H. Hall, and D.B. Chklovskii. Wiring optimization can relate neuronal structure and function PNAS 103 (12) 4723–4728 (2006)
  * [14] A.E. Pereda. Electrical synapses and their functional interactions with chemical synapses. Nat Rev Neurosci. 5(4):250-63 (2014)
  * [15] M. De Domenico, M.A. Porter, and A. Arenas. MuxViz: A Tool for Multilayer Analysis and Visualization of Networks. Journal of Complex Networks 3 (2) 159-176 (2015)
  * [16] E. Lazega. The Collegial Phenomenon: The Social Mechanisms of Cooperation Among Peers in a Corporate Law Partnership. Oxford University Press (2001)
  * [17] T.A.B. Snijders, P.E. Pattison, G.L. Robins, and M.S. Handcock. New specifications for exponential random graph models.Sociological Methodology 99-153 (2006).
  * [18] E. Omodei, M. De Domenico, A. Arenas. Characterizing interactions in online social networks during exceptional events. Front. Phys. 3, 59 (2015)
  * [19] P. Sapiezynski, A. Stopczynski, D.D. Lassen, and S. Lehmann. Interaction data from the Copenhagen Networks Study. Scientific Data 6, 315 (2019)


# Quantum simulator of extended bipartite Hubbard model with broken sublattice
symmetry: magnetism, correlations, and phase transitions

Yasser Saleem  Department of Physics, University of Ottawa, Ottawa, Canada
Amintor Dusko  Department of Physics, University of Ottawa, Ottawa, Canada
Moritz Cygorek  Department of Physics, University of Ottawa, Ottawa, Canada
Marek Korkusinski  National Research Council of Canada, Ottawa, Canada
Department of Physics, University of Ottawa, Ottawa, Canada  Pawel Hawrylak
Department of Physics, University of Ottawa, Ottawa, Canada

###### Abstract

We describe here a quantum simulator of extended bipartite Hubbard model with
broken sublattice symmetry. The simulator consists of a structured lateral
gate confining two dimensional electrons in a quantum well into artificial
minima arranged in a hexagonal lattice. The sublattice symmetry breaking is
generated by forming an artificial triangular graphene quantum dot (ATGQD)
with zigzag edges. The resulting extended Hubbard model generates tunable
ratio of tunneling strength to electron-electron interactions andbof
sublattice symmetry with control over shape. The validity of the simulator is
confirmed for small systems using mean-field and exact diagonalization many-
body approaches which show that the ground state changes from a metallic to an
antiferromagnetic (AF) phase by varying the distance between sites or depth of
the confining potential. The one-electron spectrum of these triangular dots
contains a macroscopically degenerate shell at the Fermi level. The shell
persists at the mean-field level for weak interactions (metallic phase) but
disappears for strong interactions, in the AF phase. We determine the effects
of electron-electron interactions on the ground state, the total spin, and the
excitation spectrum as a function of filling of the ATGQD. We find that the
half-filled charge neutral shell leads to a partially spin polarized state in
both metallic and AF regimes in accordance with Lieb’s theorem. In both
regimes a relatively large gap separates the spin polarized ground state to
the first excited many-body state at half filling of the degenerate shell. By
adding or removing an electron, this gap drops dramatically, and alternate
total spin states emerge with energies nearly degenerate to a spin polarized
ground state.

††preprint: APS/123-QED

##  I Introduction

There is currently interest in understanding electronic properties of strongly
correlated quantum materials often modelled by an extended Hubbard model. It
is expected that progress in solutions to this intractable problem may be
achieved with quantum simulators. Here we describe a proposal of a quantum
simulator of extended bipartite Hubbard model with broken sublattice symmetry
inspired by graphene. Much progress in quantum simulators has been achieved
with cold atoms and trapped ions [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. Progress
in solid state and photonic based simulators [12, 13, 14, 15, 16, 17, 18, 19,
20, 18, 21, 22, 23, 24, 25, 26, 27, 28] is enabled by progress in new
materials, including quasi-two-dimensional electronic systems (2DES) in
semiconductor heterojunctions [29, 30, 31, 32] and graphene. The isolation of
a single carbon layer, graphene, introduced a new 2DES with unusual electronic
properties, including the zero energy band gap, relativistic nature of
quasiparticles, sublattice pseudospin and two non-equivalent valleys[33, 34,
35, 36, 37, 38, 39, 32]. Finite lateral size quantization of graphene opens up
an energy gap making graphene a 2D atomically thin semiconductor with a gap
tunable from THz to UV [40, 41, 42, 43, 44, 45, 46, 47, 48, 49].

Some of these known properties of graphene and their previous works inspire
our approach to create a simulator based on artificially structured gates on
top of a 2D electron gas in a field effect transistor. Artificial graphene
structures have been realized already using photonic lattices, nano-
patterning, modulation doping, and scanning probe methods for atomic
manipulation on metal surfaces.[2, 24, 23, 20, 18, 21, 50, 51, 26, 22, 27, 25,
52, 53]. Here we propose a quantum simulator of graphene inspired bipartite
extended Hubbard model with broken sublattice symmetry. The working of the
simulator is confirmed for small systems using mean-field and exact
diagonalization methods which show how magnetism, correlations, and phase
transitions emerge, as parameters of the simulator are varied.

We focus here on a very important property of hexagonal lattice, i.e., the
presence of two triangular sublattices. The sublattice symmetry can be broken
in finite, triangular quantum dots with zigzag edges, and a macroscopic band
of degenerate one electron states emerges at the Fermi level. Lieb predicted
that a bipartite Hubbard model with broken sublattice symmetry will have a
finite magnetic moment[54]. Indeed, for a half-filled system, Ezawa, Fernandes
Rossier, Kaxiras, Potasz, and others [55, 56, 57, 45, 42, 44, 43] found the
ground state to be partially spin polarized in agreement with Lieb’s theorem
[54] and with exchange interaction being responsible for aligning spins of
electrons on the zero-energy degenerate shell. This additional polarization
was found to be proportional to the imbalance of the two sublattices and to
degeneracy of the shell, with extra electron spins localized largely at the
edge of the triangular structure. Such small triangular structures were
realized experimentally and confirmed theoretical predictions at half-filling
[58, 59, 60]. The presence of a degenerate macroscopic shell at the Fermi
level, analogous to the lowest Landau level , was shown to lead to strong
correlations in the ground state. It was found that the addition of a single
additional electron to a half-filled zero-energy shell destroyed the spin
polarization[43, 45].

There are several advantages in transferring the physics of graphene to an
artificial graphene-like quantum simulator where carbon atoms are replaced by
gated quantum dots and the structure is embedded in a semiconductor host.
These advantages include tunable distance between the dots and depth of
confining potential, programmable lattice symmetry and termination (i.e. edge
type), tunable electron-electron interactions and interdot tunneling [25, 27,
18]. Further advantages of artificial graphene include the ability to control
values of $U/t$ in a bipartite Hubbard model, which is not possible with
graphene [61, 62]. Such control would allow to demonstrate different
electronic phases, including transition from semi-metal to an AF insulator
[63, 12] . Additionally, triangular graphene quantum dots are susceptible to
edge reconstruction as studied by Voznyy et. al. [64]. Edge reconstruction is
responsible for smearing out the distinction between sublattices, and reduces
the quantum dot symmetry. These combined features can destroy the magnetic
properties of the system. The difficulty of edge reconstruction is overcome in
an artificial system, where the edge is determined by the external gate.
Another important advantage of artificial graphene is that, unlike in
graphene, in AG a single electron can be placed in the system in order to
probe the single particle spectrum, directly demonstrating the existence of a
zero energy shell and relating it to many-electron properties.

In this paper we report on the study of the quantum simulator of an ATGQD with
zig-zag edges. The quantum dot is formed by a structured metallic gate
generating a lattice of potential minima arranged to form a hexagonal lattice
forming a triangle with zigzag edges in a semiconductor quantum well with 2D
electron gas. We compute tunneling matrix elements and Coulomb matrix elements
as a function of AG parameters. We confirm the existence of a macroscopic
shell of degenerate states as found in triangular graphene quantum dots. Then,
we compute the Hartree-Fock (HF) spectrum for a half-filled system as a
function of the strength of Coulomb interactions as measured by the ratio of
Coulomb on-site repulsion $U$ to nearest neighbor tunneling matrix element
$t$. We find a metallic phase with a zero-energy shell at the Fermi level for
weak Coulomb interactions and an insulating AF phase, without a shell, for
strong Coulomb interactions. We find the ground state of a half-filled system
to be partially spin polarized, in agreement with Lieb’s Theorem. With
explicit calculations of $U/t$ we determine how to transition between these
two phases. We use exact diagonalization techniques to show the spin
depolarization as a function of filling factors for metallic and AF phases.
These results are presented for a small size ATGQD where numerical
calculations are possible and test the viability of the Quantum simulator for
large, intractable, systems.

The paper is organized as follows: in Section II, we describe the model,
geometry and Hamiltonian, in section III we describe the one electron tight-
binding model. In section IV we discuss the many-body Hamiltonian and the
calculation of Coulomb matrix elements. Section V discusses the HF ground
state for weak,metallic, and strong, AF, Coulomb interactions. In Section VI
we add correlations within a TB-HF-configuration-interaction method, and
demonstrate that the inclusion of an extra electron collapses the energy gap
and depolarizes the electronic system.

Figure 1: Potential profile of ATGQD with zigzag edges, with $N=97$ sites
separated by $a=15$ nm. There is broken sublattice symmetry in this system in
which the number of A sites does not equal the number of B sites, as seen in
the difference in red and blue dots.

##  II MODEL OF TRIANGULAR ARTIFICIAL GRAPHENE QUANTUM DOT

We start with electrons confined to a quantum well described by a potential
$V(z)$ where z is the growth direction. A metallic gate is deposited on a
surface at a distance “$D$” from the quantum well. The potential on the gate
generates a potential $V(\vec{r}-\vec{R})$ in the plane of a quantum well
laterally confining an electron at $\vec{r}$ in the vicinity of position
$\vec{R}$. An artificial graphene structure is defined by structuring the
metallic gate resulting in an array of $N$ confining potentials, artificial
atoms, positioned on a hexagonal lattice of potential minima at $\vec{R_{i}}$
separated by a distance $"a"$. We next introduce a back gate from which
electrons are drawn into the artificial graphene structure leaving behind a
positive, compensating charge described by a gate potential $V_{g}$. Hence,
the artificial graphene Hamiltonian describing $N_{e}$ electrons in an array
of $N$ potential traps in a quantum well, in the presence of a compensating
gate potential and including electron-electron interactions is given by

$\scalebox{1.0}{$H=\sum\limits_{i=1}\limits^{N_{e}}\left[-\frac{\hbar^{2}}{2m^{*}}\nabla^{2}_{i}+\sum\limits_{j=1}\limits^{N}v(\vec{r_{i}}-\vec{R_{j}})+v(z_{i})\right.$}\\\
\scalebox{1.0}{$\left.+V_{g}\left(\vec{r_{i}},z_{i}\right)\right]+\sum\limits_{i<j}\limits^{N_{e}}\frac{e^{2}}{\kappa|\vec{r_{i}}-\vec{r_{j}}|},$}$
(1)

Here, we sum over $N_{e}$ electrons with effective mass $m^{*}$ in the field
produced by an array of $N$ sites, interacting with a gate, confined to a
quantum well by potential $v(z_{i})$ and interacting with each other via
Coulomb potential screened by a background dielectric function $\kappa$. The
potential minima $v(\vec{r_{i}}-\vec{R_{j}})$ are finite and smooth confining
potentials. Here we describe these potentials by a Gaussian potential

$v(\vec{r_{i}}-\vec{R_{j}})=-V_{0}e^{-\frac{|\vec{r_{i}}-\vec{R_{j}}|^{2}}{d^{2}}},$
(2)

with depth $V_{0}$, confinement length $d$, localized in the plane of a
quantum well at $\vec{R_{j}}$. The one electron potential
$V(r)=\sum_{j}v(\vec{r}-\vec{R_{j}})$ of the artificial graphene structure
studied here is shown in Fig. 1. Different structures with different size and
shape can be constructed analogously. Here, there are $N=97$ sites, with each
modelled as a Gaussian confining potential, with a depth $V_{0}=300$ meV,
confinement length $d=10$ nm and separation of $a=15$ nm. We see that the
confining potential forms a triangular quantum well, with minima arranged on a
hexagonal lattice, with visible benzene like rings and terminated by zigzag
edges. Such a structure is an example of a bipartite lattice with broken
sublattice symmetry, and as such Lieb’s theorem [54] will apply and play a
critical role in determining the nature of the ground state.

##  III ONE ELECTRON SPECTRUM IN THE TIGHT-BINDING MODEL

We now introduce one electron into the artificial graphene structure shown in
Fig. 1. Following the model described in Section II, the single particle
Hamiltonian is given by

$H_{0}=-\frac{\hbar^{2}}{2m^{*}}\nabla^{2}+\sum_{j}v(\vec{r}-\vec{R_{j}})+v(z),$
(3)

where the sum over $j$ extends over $N$ sites, and $v(\vec{r}-\vec{R_{j}})$ is
given in Eq. (2). Here $v(z)$ is a potential of an infinite quantum well with
width $0.1a$. Throughout this article, we assume strong confinement in the
$z$-direction, so that the eigenstates factorize into an in-plane part and
part in the z-direction, which we assume to be the lowest state of an infinite
quantum well $\xi\left(z\right)=\sqrt{\frac{2}{L}}\sin{\frac{\pi z}{L}}$. The
lateral confining potential is smooth and parabolic at low energies. Hence we
expand the in-plane part of the wavefunction $\varphi$ for a state $\nu$ in
terms of two-dimensional harmonic oscillator eigenfunctions $\alpha$ centered
on atom $j$:

$\varphi^{\nu}=\sum_{j,\alpha}A^{\nu}_{j\alpha}\phi^{0}_{j\alpha}.$ (4)

Acting with $H_{0}$ in Eq.(3) on the wavefunction Eq.(4), and noting that we
can decompose the Gaussian in terms of a parabolic confining potential plus a
correction as $v(\vec{r}-\vec{R_{j}})=v_{j}^{ho}+\delta V_{j}$ with
$v_{j}^{ho}=-V_{0}+\frac{\overline{\omega}^{2}}{4}|\vec{r}-\vec{R_{j}}|^{2}$
and $\overline{\omega}^{2}=\frac{4}{d^{2}}V_{0}$, we get matrix elements of
the Hamiltonian in Eq. (3) given by

$\scalebox{1.0}{$H_{m\beta,l\alpha}=\epsilon^{0}_{\beta}S_{m\beta,l\alpha}+\braket{\phi^{0}_{m\beta}}{\delta
V_{m}}{\phi^{0}_{l\alpha}}$}\\\ \scalebox{1.0}{$+\sum_{j\neq
m}\braket{\phi^{0}_{m\beta}}{v^{ho}_{j}+\delta V_{j}}{\phi^{0}_{l\alpha}}.$}$
(5)

Here $S_{m\beta,l\alpha}=\langle\phi^{0}_{m\beta}|\phi^{0}_{l\alpha}\rangle$
are overlap matrix elements for orbitals $\beta,\alpha$ localized on sites
$m,l$. Since the harmonic oscillator states on different sites defined in
Eq.(4) are not orthogonal, we orthogonalize the basis by solving the
generalized eigenvalue problem given by

$\textbf{\~{H}}\textbf{B}^{\nu}=\epsilon^{\nu}\textbf{B}^{\nu},$ (6)

where $\textbf{B}^{\nu}=\textbf{S}^{1/2}\textbf{A}^{\nu}$ and S is the overlap
matrix. The corresponding renormalized Hamiltonian given by

$\textbf{\~{H}}=\textbf{S}^{-\frac{1}{2}}\textbf{H}\textbf{S}^{-\frac{1}{2}}.$
(7)

The wavefunction is now expanded in terms of orthogonal states $\psi_{m\beta}$
localized on different sites

$\phi^{\nu}=\sum_{m,\beta}B^{\nu}_{m\beta}\psi_{m\beta}.$ (8)

The orthogonal orbitals localized on site $"m"$ are given explicitly by

$\psi_{m\beta}=\sum_{l\alpha}(S^{-\frac{1}{2}})_{l\alpha,m\beta}\phi^{0}_{l\beta}.$
(9)

Figure 2: TB spectrum for AG Triangular quantum dot with zigzag edges near the
Fermi level.

Fig. 2 shows the energy levels of a single electron obtained by diagonalizing
Eq.(7) for $N=97$ sites in the basis of 3 harmonic oscillator shells, $S$,
$P$, and $D$ per site, with separation of potential minima corresponding to
$a=12.5$ nm, and all other parameters given in section II. We observe a well-
defined shell of almost degenerate states at the Fermi level contained in the
1$S$-band. We note that the Fermi level corresponds to a half-filled electron
system. A large gap separates the 1$S$ band of levels from the $P$, and $D$
derived levels and the single particle spectrum resembles a spectrum obtained
by diagonalizing a tight binding model with nearest neighbor tunneling matrix
elements [43]. Thus we limit ourselves to only $S$ orbitals from here on in
order to simplify many-body calculations. We take the hopping parameter “$t$”
between nearest neighbour 1$S$ orbitals to be equal to half the bandwidth of
the 1$S$ band. We find $t=26.98$ meV for site separation of $a=12.5$ nm and
$t=3.18$ meV for site separation of $a=15$ nm. In Fig.2 we see a shell of
nearly degenerate zero-energy states at the Fermi level split by the
introduction of next-nearest-neighbour hopping generated by the itinerant
orbitals. Note that in artificial graphene, it is possible to have a single
electron moving in the system of potential minima, and optical experiments,
for example, could probe the existence of this zero energy shell.

##  IV Many-body Hamiltonian

With the orthogonalized orbitals $\psi_{m\beta}$ limited to the 1$S$ band, we
can write the many-body Hamiltonian, Eq. (1) in the second quantized form as

$\scalebox{1.0}{$H=\sum\limits_{i,\sigma}\epsilon_{i\sigma}c_{i\sigma}^{\dagger}c_{i\sigma}+\sum\limits_{i,j,\sigma}t_{ij}c_{i\sigma}^{\dagger}c_{j\sigma}+$}\\\
\scalebox{1.0}{$\frac{1}{2}\sum\limits_{i,j,k,l,\sigma,\sigma^{\prime}}\braket{ij}{V}{kl}c_{i\sigma}^{\dagger}c_{j\sigma^{\prime}}^{\dagger}c_{k\sigma^{\prime}}c_{l\sigma}+\sum\limits_{i,\sigma}v_{ii}^{g}c_{i\sigma}^{\dagger}c_{i\sigma},$}$
(10)

where $\epsilon_{i\sigma}$ are the onsite energies and $t_{ij\sigma}$ are the
hopping matrix elements computed above. $v_{ii}^{g}$ corresponds to the back
gate and $\braket{ij}{V}{kl}$ are Coulomb matrix elements given by

$\scalebox{1.0}{$\braket{ij}{V}{kl}=\int\int
d\vec{r}_{1}d\vec{r}_{2}\psi^{*}_{i}(\vec{r}_{1})\xi_{i}^{*}(z_{1})\psi^{*}_{j}(\vec{r}_{2})\xi_{j}^{*}(z_{2})$}\\\
\scalebox{1.0}{$\times\frac{2R_{y}}{|\vec{r}_{1}-\vec{r}_{2}+(z_{1}\hat{z}-z_{2}\hat{z})|}\psi_{k}(\vec{r}_{2})\xi_{k}(z_{2})\psi_{l}(\vec{r}_{1})\xi_{l}(z_{1}),$}$
(11)

where the functions $\psi_{i}$ are the 1$S$ localized and orthogonal orbitals
which are defined in Eq.(9). The functions
$\xi_{i}=\sqrt{\frac{2}{L}}sin(\frac{\pi z_{i}}{L})$ describe the lowest
energy sub-band of an infinite quantum well confining the electrons in the
$z$-direction. The back gate term is given by

$v_{ii}^{g}=\sum_{j}\frac{-\frac{2N_{p}}{N}}{\sqrt{(x_{i}-x_{j})^{2}+(y_{i}-y_{j})^{2}+d_{gate}^{2}}},$
(12)

where $N_{p}$ is the number of positive charges on the gate. The model assumes
that the number of positive charges on the gate $N_{p}$ is equal to the number
of electron charges in the system $N$ in order to enforce charge neutrality,
and is uniformly smeared on the gate.

Since the Hamiltonian in Eq.(10) cannot be solved exactly, we start by solving
the mean-field Hartree-Fock problem first.

##  V Mean-field Hartree-Fock

Figure 3: a) HF Spectrum for ground state of N=97 sites at half-filling with
lattice constant $a=12.5$ nm. We observe a shell of distinct degenerate states
near the Fermi level, filling up of levels with spin up and down electrons,
and a splitting of spin degeneracy of all levels due to an imbalance of up and
down spins found in HF solution. All states below the Fermi level ${E_{F}=0}$
(middle line) are occupied, all states above are unoccupied. b) HF spectrum
for lattice separation $a=15$ nm. Note the disappearance of a degenerate shell
at the Fermi level and emergence of a large gap proportional to Hubbard U
separating the valence and conduction band states, AF spin ordering and
partially spin polarized ground state. Figure 4: a) Spin density obtained for
the spin polarized HF ground state for $a=12.5$ nm. We observe the system to
be in the metallic phase, with the extra spins occupying the edge. b) Spin
density for the spin polarized ground state for $a=15$ nm . We notice one spin
is localized on one sublattice, while the other, is localized on the other
sublattice. Due to a sublattice imbalance this as well leads to spin polarized
electrons localized on the edge.

The magnetic properties of artificial graphene, and the nature of the ground
state for different electron numbers, can be studied by solving the Mean-field
HF Hamiltonian obtained from Eq.(10) and given by:

$\scalebox{1.0}{$H^{ATGQD}_{MF}=\sum\limits_{i,l,\sigma}\tau^{0}_{il\sigma}c_{i\sigma}^{\dagger}c_{l\sigma}$}\\\
\scalebox{1.0}{$+\sum\limits_{i,j,k,l,\sigma,\sigma^{\prime}}\left[\braket{ij}{V}{kl}-\braket{ij}{V}{lk}\delta_{\sigma\sigma^{\prime}}\right]$}\\\
\scalebox{1.0}{$\times\left(\rho_{jk\sigma^{\prime}}-\rho^{0}_{jk\sigma^{\prime}}\right)c_{i\sigma}^{\dagger}c_{l\sigma}+\sum\limits_{i,\sigma}v_{ii}^{g}c_{i\sigma}^{\dagger}c_{i\sigma},$}$
(13)

where $\rho_{jk\sigma^{\prime}}$ is the density matrix elements for the ATGQD,
and $\rho^{0}_{jk\sigma^{\prime}}$ is the density matrix for the bulk system.
The choice of $\rho^{0}_{jk\sigma^{\prime}}$ is discussed in appendix A, and
its purpose is to enhance convergence. $\tau^{0}_{il\sigma}$ is a tunneling
matrix element which describes the properties of bulk artificial graphene in
terms of the tunneling matrix element $t_{il}$ and bulk density matrix
$\rho^{0}_{jk\sigma}$. It is given by

$\tau^{0}_{il\sigma}=-t_{il\sigma}+\sum_{jk\sigma^{\prime}}\left[\braket{ij}{V}{kl}-\braket{ij}{V}{lk}\delta_{\sigma\sigma^{\prime}}\right]\rho^{0}_{jk\sigma^{\prime}},$
(14)

Here the superscript $0$ on $\tau^{0}_{il\sigma}$ is to denote that it is
computed with the bulk density matrix elements. In addition to the onsite
interaction terms, all direct terms are taken into account, as well as all
exchange terms up to next nearest neighbours. The bulk density matrix elements
are computed in the appendix for the AFM regime. They have been obtained by
Potasz et al. [45] for the metallic regime. We note the Hamiltonian in Eq.13
is symmetric with respect to spin, and thus can be diagonalized in separate
subspaces for spin up and spin down, but with the spin up Hamiltonian having a
dependence on the density of spin down electrons and vice versa. We focus on
$S_{z}\geq 0$ and proceed to diagonalize Eq.(13), with results shown in Fig.3
and Fig.4. Fig.3(a) shows the energy spectrum for spin up and down electrons
in the metallic regime , $a=12.5$ nm. We see a spin splitting of levels due to
a spin imbalance obtained in HF and consistent with Lieb’s theorem. We observe
a nearly degenerate shell at the Fermi level, with the blue spin up electrons
fully occupying a degenerate shell, leaving the red spin down levels
completely empty above the Fermi level. These extra spins are found to align
on the edge of the triangle as seen in Fig.4(a), with a uniform zero spin
density away from the edges indicating a semi-metallic regime. We now proceed
to strongly interacting regime by increasing the distance between lattice
sites to $a=15$ nm. We find the ground state to be again partially spin
polarized due to broken sublattice symmetry in agreement with Lieb’s theorem
[54]. We however lose a distinct degenerate shell at the Fermi level, instead
we find a large gap, proportional to Hubbard U, separating the valence and
conduction band, suggesting an insulating phase (Fig.3(b)). We observe an
antiferromagnetic spin ordering in the bulk as shown in Fig.4(b) as expected
in the large U/t regime [63] and ferromagnetic ordering at the edges of the
ATGQD.

Figure 5: The low-energy spectra in the metallic phase for a) 97-1 electrons
b) 97 electrons c) 97+1 electrons. For 97 electrons the ground state is
partially spin polarized, and we see a large gap that separates this state
from other spin states, but the introduction or removal of an electron
collapses this gap and many spins states exist at very close energies to the
spin polarized ground state for the $N=97-1$ and $N=97+1$ electron cases.

##  VI ELECTRONIC CORRELATIONS VIA CONFIGURATION INTERACTION

###  VI.1 ELECTRONIC CORRELATIONS IN THE METALLIC PHASE

We now turn to include electronic correlations. Let us begin with the semi-
metallic phase. In the example of ATGQD with $N=97$ sites, the breaking of
sublattice symmetry results in $N_{d}=7$ fold degenerate shell as shown in
Figure 3. The electronic correlations are most important for electrons
occupying the degenerate shell and we will treat them using configuration-
interaction (CI) method. The remaining background electrons are treated in HF.
Hence we proceed to solve the HF problem for $97-7=90$ electrons on 97 sites
with $N_{\uparrow}=N_{\downarrow}=45$. This leaves the degenerate shell at the
Fermi level empty. We then rotate the many body Hamiltonian, Eq.10, to the HF
basis[65, 45] for $N=90$ with final result given by

$\scalebox{1.0}{$H=\sum\limits_{p,\sigma}\epsilon^{HF}_{p\sigma}b_{p\sigma}^{\dagger}b_{p\sigma}-\sum\limits_{p,q,\sigma}t_{pq\sigma}b_{p\sigma}^{\dagger}b_{q\sigma}$}\\\
\scalebox{1.0}{$+\frac{1}{2}\sum\limits_{p,q,r,s,\sigma,\sigma^{\prime}}\braket{pq}{V}{rs}b_{p\sigma}^{\dagger}b_{q\sigma^{\prime}}^{\dagger}b_{r\sigma^{\prime}}b_{s\sigma}+\sum\limits_{p,q,\sigma}v^{add}_{pq}b_{p\sigma}^{\dagger}b_{q\sigma}$}$
(15)

where

$t_{pq\sigma}=\sum\limits_{i,l}\tau_{il\sigma}a^{*}_{il\sigma}a_{lq\sigma}$
(16)

and

$b_{q\sigma}=\sum\limits_{l}a_{lq\sigma}c_{l\sigma}$ (17)

with $\tau_{il\sigma}$ defined in Eqn.14 but computed with respect to ATGQD
density matrix elements. $a_{lq\sigma}$ are the eigenvectors obtained by
diagonalizing Eqn.13. We note that term $t_{pq\sigma}$ appears, in order to
lower the contribution of the quasiparticle-quasiparticle interaction term
$\braket{pq}{V}{rs}$. $\braket{pq}{V}{rs}$ are Coulomb matrix elements in the
basis of HF states, and are computed by rotating the real space matrix
elements in Eq.11, to the basis of HF states. They describe the remaining
interaction of HF quasiparticles beyond the mean-field. The last term in Eq.15
involving $v_{pq}^{add}$ describes additional HF quasiparticles added to the
degenerate shell. Since we solve the HF problem for $N_{e}=90$ electrons, when
adding HF quasiparticles at the CI level, we must compensate this charge with
additional positive charges on the gate to maintain charge neutrality. Since a
large gap separates the nearly degenerate shell from other states, it suffices
to take only the shell near the Fermi level for CI calculations and neglect
scattering from the valence band to the shell or from the shell to the
conduction band. Fig.5 shows the low energy spectra obtained by diagonalizing
Eq.15 for half-filled system with $N=97$ , with extra electron ($N=98$) and
with extra hole ($N=96)$. Focusing on the half-filled (N=97 electrons) case in
Fig.5(b), we see the ground state of the half-filled shell to be maximally
spin polarized, in agreement with Lieb’s theorem. The energy of this
configuration is well separated from other states with lower total spin S, in
other words the energy gap between our ground state, and first excited state
with a different total spin $S$ is large. This implies that the energy cost to
flip a spin is large. The removal or addition of a single electron, Fig.5(a)
and Fig.5(c),results in a ground state which is still maximally spin
polarized, but other low spin states lower their energy due to correlations,
with many total S states very close in energy. It costs practically zero
energy to flip a spin in this case. In contrast to regular graphene [45] where
the ground state corresponded to S=0, the ground state here has S=3, but we
observe a dramatic drop of the energy gap between different total S states, a
phenomenon seen in graphene as well [45, 43, 42]. Correlations in the lower
spin states cause a decrease of the energies, they become almost degenerate in
energy with the maximum spin state. It is worth noting that the spin of the
ground state for $N=96$ electrons is in agreement with [45], and the shrinking
of the gap between different total spins states is consistent as well.

###  VI.2 ELECTRONIC CORRELATIONS IN THE ANTI-FERROMAGNETIC PHASE

Due to the very different quasiparticle spectra for the AF phase and the semi-
metallic phase as seen in Fig.3, we require two different approaches to the
many body problem. Here we solve the many body problem in the AF phase in real
space by improving on the HF ground state in real space. In the AF phase there
is no degenerate shell at the Fermi level. Hence, we begin with the HF
solution in the AF phase. We obtain the HF solution for a fixed number of spin
up and spin down electrons by diagonalizing Eq.13, yielding a single Slater
determinant defined as

$\ket{HFGS}=\prod_{q=1}^{\lambda_{F}^{\uparrow}}b^{{\dagger}}_{q\uparrow}\prod_{q=1}^{\lambda_{F}^{\downarrow}}b^{{\dagger}}_{q\downarrow}\ket{0},$
(18)

where the Slater determinant is defined by filling up HF quasi particle levels
up to the Fermi level for each spin. We then rotate this HF state to the site
basis. In the site representation, we have a linear combination of Slater
determinants, and we select the largest contributing state. For example, at
half-filling for $S_{z}=\frac{7}{2}$, the ground state is given by

$\ket{GS}=\prod_{i\in A}c^{{\dagger}}_{i\uparrow}\prod_{i\in
B}c^{{\dagger}}_{i\downarrow}\ket{0},$ (19)

where we place spin up electrons on the A sublattice, and spin down electrons
on the B sublattice, representing a perfect antiferromagnetic phase. This is
the largest dominant real space configuration composing that HF groundstate
seen in Fig.4(b). We then divide the Hilbert space into segments for different
$S_{z}$ subspaces. This is done by starting with the ground state for
different $S_{z}$, as shown above and constructing configurations with the
same total $S_{z}$. The Hilbert space is divided into 5 sets of configurations
defined by

Figure 6: Energy Gap for different number of electrons as a function of
increasing Hilbert space size. In green, we are at half-filling, and the gap
between the ground state and the next spin state is very large, while for the
case where we have added or removed an electron from the system, the energy
gap collapses to almost zero. Figure 7: Energy gap vs filling factor, for AFM
(black) and metallic regime (red). In both cases a similar behaviour is
observed in which at half-filling the gap is large, and collapses away from
half-filling.

$\ket{O_{1,0}}=\sum\limits_{\begin{subarray}{c}i,j,\sigma\\\
{<0DO>}\end{subarray}}c_{j\sigma}^{\dagger}c_{i\sigma}\ket{GS}$ (20a)
$\ket{O_{1,1}}=\sum\limits_{\begin{subarray}{c}i,j,\sigma\\\
{<1DO>}\end{subarray}}c_{j\sigma}^{\dagger}c_{i\sigma}\ket{GS}$ (20b)
$\ket{O_{2,0}}=\sum\limits_{\begin{subarray}{c}i,j,k,l,\sigma,\sigma^{\prime}\\\
{<0DO>}\end{subarray}}c_{l\sigma}^{\dagger}c_{k\sigma^{\prime}}^{\dagger}c_{j\sigma^{\prime}}c_{i\sigma}\ket{GS}$
(20c)
$\ket{O_{2,1}}=\sum\limits_{\begin{subarray}{c}i,j,k,l,\sigma,\sigma^{\prime}\\\
{<1DO>}\end{subarray}}c_{l\sigma}^{\dagger}c_{k\sigma^{\prime}}^{\dagger}c_{j\sigma^{\prime}}c_{i\sigma}\ket{GS}$
(20d)
$\ket{O_{2,1}}=\sum\limits_{\begin{subarray}{c}i,j,k,l,\sigma,\sigma^{\prime},\\\
{<2DO>}\end{subarray}}c_{l\sigma}^{\dagger}c_{k\sigma^{\prime}}^{\dagger}c_{j\sigma^{\prime}}c_{i\sigma}\ket{GS},$
(20e)

where the brackets under the sum denote a restriction of the configurations in
which we include either 0,1 or 2 double occupancies (DO) measured from the
number of DO of the ground state. The subscripts of $\ket{O_{\mu,\nu}}$ are
defined with the first number denoting the number of electrons moved, and the
second number being the number of extra double occupancies. We also restrict
configurations by allowing only nearest neighbor scattering describing
correlations. We then proceed to diagonalize the Hamiltonian in Eq.10, as it
is written in the site basis, and not in the HF basis. We take the same
Coulomb matrix elements as described at the HF stage. After diagonalization,
we observe a maximally spin polarized ground state in the case of $N=96,97$
and $98$ electrons. In the case of adding or removing an electron, the
energies are very close to each other. In Fig.6 we observe that in the half-
filled case we have a large gap separating this state from any other spin
state. Meanwhile, when we add or remove an electron the energy gap between the
spin polarized ground state, and the spin state closest in energy collapses to
nearly zero. We now compare the metallic phase and AF phases. We first note
that in the metallic phase, we see an emergence of a shell at zero energy of
states which are practically degenerate. Here we are guided by the intuition
that exchange interaction will lower the energy of a spin polarized half-
filled system thus expecting Lieb’s theorem to be valid. In the AF regime,
without the presence of a degenerate shell, we still see a spin polarized
ground state. We then explore the regimes where we move away from half-
filling, where Lieb’s theorem doesn’t need to hold. In both regimes, we see a
collapse of the energy gap. This is consistent with previous work [45, 43].
Fig.7 shows schematically the collapse of the energy gap away from half-
filling. The gap peaks at half-filling and collapses when we add/remove an
electron. We note that in the metallic regime, we study this as a function of
filling factor, and see that the gap is maximum at half filling, where the
spin polarized ground state is most stable.

##  VII SUMMARY

In summary, we described here how to construct a quantum simulator of an
extended bipartite Hubbard model with broken sublattice symmetry using a
structured lateral gate confining two dimensional electrons in a quantum well
into artificial minima arranged in a hexagonal lattice. The sublattice
symmetry breaking was generated by forming an artificial triangular graphene
quantum dot (ATGQD) with zigzag edges. We demonstrated that in artificial
graphene quantum dots, by tuning $U/t$, we can reach two distinct regimes, a
semi-metallic and AFM. We showed for small systems that in both the metallic
and AF regimes, the system at half-filling is partially spin polarized in
agreement with Lieb’s theorem. The addition or removal of an electron in both
regimes collapses the energy gap and spin polarization. Such a simulator would
allow simulation of larger systems, verification of results presented here and
potential discovery of new phases resulting from strong electron-electron
interactions in hexagonal lattice systems inherent in graphene and transition
metal dichalcogenites.

###### Acknowledgements.

Y.S., A.D.,M.C., M.K. and P.H. acknowledges support from NSERC QC2DM Strategic
Grant No. STPG-521420, NSERC Discovery Grant No. RGPIN 2019-05714 and
University of Ottawa Research Chair in Quantum Theory of Quantum Materials,
Nanostructures, and Devices. Y.S. thanks J.Manalo, K.Sadecka, L.Szulakowska
for useful discussions.

*

## Appendix A

In this appendix we calculate the density matrix elements
$\rho_{jk\sigma^{\prime}}^{0}$ in the AFM (see [45] for derivation in metallic
phase). We start with the mean field Hamiltonian for bulk artificial graphene
is given by

$\scalebox{1.0}{$H^{0}_{MF}=\sum\limits_{i,l,\sigma}t_{il\sigma}c_{i\sigma}^{\dagger}c_{l\sigma}$}\\\
\scalebox{1.0}{$+\sum\limits_{i,j,k,l,\sigma,\sigma^{\prime}}\left[\braket{ij}{V}{kl}-\braket{ij}{V}{lk}\delta_{\sigma\sigma^{\prime}}\right]\rho^{0}_{jk\sigma^{\prime}}c_{i\sigma}^{\dagger}c_{l\sigma}$}$
(21)

Now, if we take only onsite density matrix elements, and we keep only terms
where $i=l$ (i.e. ignore small scattering elements) we can write the
Hamiltonian as

$H^{0}_{MF}=\sum\limits_{i,l,\sigma}t_{il}c_{i\sigma}^{\dagger}c_{l\sigma}+\sum\limits_{i,\sigma}\Delta_{i\sigma}c_{i\sigma}^{\dagger}c_{i\sigma}$
(22)

Where

$\Delta_{i\sigma}=\sum\limits_{j\sigma^{\prime}}\left[\braket{ij}{V}{ji}-\braket{ij}{V}{ij}\delta_{\sigma\sigma^{\prime}}\right]\rho^{0}_{jj\sigma^{\prime}}$
(23)

Postulating that we are in the AFM regime, so that

$\rho_{jj\uparrow}^{0}=\begin{cases}1&j\in\text{Subblatice A}\\\
0&j\in\text{Subblatice B}\end{cases}$ (24a)
$\rho_{jj\downarrow}^{0}=\begin{cases}0&j\in\text{Subblatice A}\\\
1&j\in\text{Subblatice B}\end{cases}$ (24b)

We arrive at the expressions

$\Delta_{A\uparrow}=\sum\limits_{j>1}\braket{1j}{V}{j1}$ (25a)
$\Delta_{A\downarrow}=U+\sum\limits_{j>1}\braket{1j}{V}{j1}$ (25b)
$\Delta_{B\uparrow}=U+\sum\limits_{j>1}\braket{1j}{V}{j1}$ (25c)
$\Delta_{B\downarrow}=\sum\limits_{j>1}\braket{1j}{V}{j1}$ (25d)

Where $\braket{1j}{V}{j1}$ corresponds to a direct interaction with a fixed
site we call $i=1$, and another site j which is independent of which
sublattice we fix our state $i=1$ to since all sites are identical. Now we can
subtract the constant $\frac{U}{2}+\sum\limits_{i>1}\braket{1i}{V}{i1}$ from
all terms since this leads to only a shift in energies we get

$\Delta_{A\uparrow}=-\frac{U}{2}$ (26a) $\Delta_{A\downarrow}=\frac{U}{2}$
(26b) $\Delta_{B\uparrow}=\frac{U}{2}$ (26c)
$\Delta_{B\downarrow}=-\frac{U}{2}$ (26d)

Taking nearest neighbours hopping only, the Hamiltonian in the basis of
sublattices $A\uparrow$,$B\uparrow$,$A\downarrow$,$B\downarrow$ is given by

$H\left(\vec{k}\right)=\begin{bmatrix}-\frac{U}{2}&-tf(\vec{k})&0&0\\\
-tf^{*}(\vec{k})&\frac{U}{2}&0&0\\\ 0&0&\frac{U}{2}&-tf(\vec{k})\\\
0&0&-tf^{*}(\vec{k})&-\frac{U}{2}\end{bmatrix}$ (27)

This Hamiltonian is block diagonal, with each block mimicking gapped graphene.
The valence band solution is given by

$E_{k}^{-}=-\sqrt{\left(\frac{U}{2}\right)^{2}+t^{2}|f(\vec{k})|^{2}}$ (28)
$\scalebox{1.0}{$\psi_{k}{{}_{-}}^{(\downarrow)}=-\sin\frac{\varphi_{k}}{2}e^{i\theta_{k}}\frac{1}{\sqrt{N_{U}}}\sum\limits_{\vec{R}_{A}}e^{i\vec{k}\cdot\vec{R}_{A}}\phi_{z}\left(\vec{r}-\vec{R}_{A}\right)$}\\\
+\scalebox{1.0}{$\cos\frac{\varphi_{k}}{2}\frac{1}{\sqrt{N_{U}}}\sum\limits_{\vec{R}_{B}}e^{i\vec{k}\cdot\vec{R}_{B}}\phi_{z}\left(\vec{r}-\vec{R}_{B}\right)$}$
(29a)
$\scalebox{1.0}{$\psi_{k-}^{(\uparrow)}=\cos\frac{\varphi_{k}}{2}\frac{1}{\sqrt{N_{U}}}\sum\limits_{\vec{R}_{A}}e^{i\vec{k}\cdot\vec{R}_{A}}\phi_{z}\left(\vec{r}-\vec{R}_{A}\right)$}\\\
\scalebox{1.0}{$-\sin\frac{\varphi_{k}}{2}\frac{1}{\sqrt{N_{U}}}\sum\limits_{\vec{R}_{B}}e^{i\vec{k}\cdot\vec{R}_{B}}\phi_{z}\left(\vec{r}-\vec{R}_{B}\right)$}$
(29b)

Where $\frac{U}{2}=|E^{-}_{k}|\cos{\varphi_{k}}$ defines $\varphi_{k}$ , and
$f(\vec{k})$ is the usual form factor of graphene. The density Matrix elements
are defined as

$\rho_{ij\sigma}^{0}=\sum_{\vec{k}}b^{*}_{\vec{R}_{i},\vec{k},\sigma}b_{\vec{R}_{j},\vec{k},\sigma}$
(30)

Where $b_{\vec{R}_{j},\vec{k},\sigma}$ are the coeficients of the wavefunction
defined in eqn.(29b). Computing the density matrix we have

$\rho_{BB\downarrow}^{0}=\rho_{AA\uparrow}^{0}=\sum_{\vec{k}}\cos^{2}{\frac{\varphi}{2}}\approx
1$ (31a)
$\rho_{BB\uparrow}^{0}=\rho_{AA\downarrow}^{0}=\sum_{\vec{k}}\left(1-\cos^{2}{\frac{\varphi}{2}}\right)\approx\frac{t}{U}$
(31b)
$\rho_{BB\uparrow}^{0}=\rho_{AA\downarrow}^{0}=\sum_{\vec{k}}\sin{\frac{\varphi}{2}}\cos{\frac{\varphi}{2}}e^{-i\theta_{k}}e^{i\vec{k}\cdot{b}}\approx
0$ (31c)

Where the approximate solution exploits the fact that we are in the AF phase
and so $\frac{t}{U}<<1$. This is in agreement with the numerical calculation
of the density matrix elements in a finite flake in the center of the
structure.

## References

  * Lloyd [1996] S. Lloyd, Science , 1073 (1996).
  * Bloch [2005] I. Bloch, Nature physics 1, 23 (2005).
  * Jaksch and Zoller [2005] D. Jaksch and P. Zoller, Annals of physics 315, 52 (2005).
  * Ortner  _et al._ [2009] M. Ortner, A. Micheli, G. Pupillo, and P. Zoller, New Journal of Physics 11, 055045 (2009).
  * Mazurenko  _et al._ [2017] A. Mazurenko, C. S. Chiu, G. Ji, M. F. Parsons, M. Kanász-Nagy, R. Schmidt, F. Grusdt, E. Demler, D. Greif, and M. Greiner, Nature 545, 462 (2017).
  * Kuhr [2016] S. Kuhr, National Science Review 3, 170 (2016).
  * Weimer  _et al._ [2010] H. Weimer, M. Müller, I. Lesanovsky, P. Zoller, and H. P. Büchler, Nature Physics 6, 382 (2010).
  * Barreiro  _et al._ [2011] J. T. Barreiro, M. Müller, P. Schindler, D. Nigg, T. Monz, M. Chwalla, M. Hennrich, C. F. Roos, P. Zoller, and R. Blatt, Nature 470, 486 (2011).
  * Hempel  _et al._ [2018] C. Hempel, C. Maier, J. Romero, J. McClean, T. Monz, H. Shen, P. Jurcevic, B. P. Lanyon, P. Love, R. Babbush, _et al._ , Physical Review X 8, 031022 (2018).
  * Islam  _et al._ [2011] R. Islam, E. Edwards, K. Kim, S. Korenblit, C. Noh, H. Carmichael, G.-D. Lin, L.-M. Duan, C.-C. J. Wang, J. Freericks, _et al._ , Nature communications 2, 1 (2011).
  * Leykam  _et al._ [2018] D. Leykam, A. Andreanov, and S. Flach, Advances in Physics: X 3, 1473052 (2018).
  * Salfi  _et al._ [2016] J. Salfi, J. Mol, R. Rahman, G. Klimeck, M. Simmons, L. Hollenberg, and S. Rogge, Nature communications 7, 1 (2016).
  * Buluta and Nori [2009] I. Buluta and F. Nori, Science 326, 108 (2009).
  * Aspuru-Guzik and Walther [2012] A. Aspuru-Guzik and P. Walther, Nature physics 8, 285 (2012).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, _et al._ , Nature 551, 579 (2017).
  * Cai  _et al._ [2013] J. Cai, A. Retzker, F. Jelezko, and M. B. Plenio, Nature Physics 9, 168 (2013).
  * Dusko  _et al._ [2018] A. Dusko, A. Delgado, A. Saraiva, and B. Koiller, npj Quantum Information 4, 1 (2018).
  * Singha  _et al._ [2011] A. Singha, M. Gibertini, B. Karmakar, S. Yuan, M. Polini, G. Vignale, M. Katsnelson, A. Pinczuk, L. Pfeiffer, K. West, _et al._ , Science 332, 1176 (2011).
  * Klembt  _et al._ [2018] S. Klembt, T. Harder, O. Egorov, R. G. K Winkler, M. Bandres, M. Emmerling, L. Worschech, T. Liew, M. Segev, C. Schneider, and S. Höfling, Nature 562, 552 (2018).
  * Park and Louie [2009] C.-H. Park and S. G. Louie, Nano letters 9, 1793 (2009).
  * Li  _et al._ [2021] Y. Li, S. Dietrich, C. Forsythe, T. Taniguchi, K. Watanabe, P. Moon, and C. R. Dean, Nature Nanotechnology 16, 525 (2021).
  * Shi  _et al._ [2019] L.-k. Shi, J. Ma, and J. C. Song, 2D Materials 7, 015028 (2019).
  * Forsythe  _et al._ [2018] C. Forsythe, X. Zhou, K. Watanabe, T. Taniguchi, A. Pasupathy, P. Moon, M. Koshino, P. Kim, and C. R. Dean, Nature nanotechnology 13, 566 (2018).
  * Gibertini  _et al._ [2009] M. Gibertini, A. Singha, V. Pellegrini, M. Polini, G. Vignale, A. Pinczuk, L. N. Pfeiffer, and K. W. West, Physical Review B 79, 241406 (2009).
  * Kylänpää  _et al._ [2016] I. Kylänpää, F. Berardi, E. Räsänen, P. García-González, C. A. Rozzi, and A. Rubio, New Journal of Physics 18, 083014 (2016).
  * Uehlinger  _et al._ [2013] T. Uehlinger, G. Jotzu, M. Messer, D. Greif, W. Hofstetter, U. Bissbort, and T. Esslinger, Physical review letters 111, 185307 (2013).
  * Räsänen  _et al._ [2012] E. Räsänen, C. Rozzi, S. Pittalis, and G. Vignale, Physical review letters 108, 246803 (2012).
  * Jacqmin  _et al._ [2014] T. Jacqmin, I. Carusotto, I. Sagnes, M. Abbarchi, D. Solnyshkov, G. Malpuech, E. Galopin, A. Lemaître, J. Bloch, and A. Amo, Physical Review Letters 112, 116402 (2014).
  * Editors  _et al._ [2000] Editors, EP2DS1999, P. Hawrylak, D. Lockwood, and A. Sachrajda, Physica E 6 (2000).
  * Editors  _et al._ [2010] Editors, EP2DS2009, S. Katsumoto, K. Kono, and S. Tarucha, Physica E 42, 673 (2010).
  * Piquero-Zulaica  _et al._ [2017] I. Piquero-Zulaica, J. Lobo-Checa, A. Sadeghi, Z. M. Abd El-Fattah, C. Mitsui, T. Okamoto, R. Pawlak, T. Meier, A. Arnau, J. E. Ortega, _et al._ , Nature communications 8, 1 (2017).
  * Wang  _et al._ [2014] S. Wang, L. Z. Tan, W. Wang, S. G. Louie, and N. Lin, Physical review letters 113, 196803 (2014).
  * Wallace [1947] P. R. Wallace, Physical review 71, 622 (1947).
  * Novoselov  _et al._ [2004] K. S. Novoselov, A. K. Geim, S. V. Morozov, D.-e. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, science 306, 666 (2004).
  * Zhang  _et al._ [2005] Y. Zhang, Y.-W. Tan, H. L. Stormer, and P. Kim, nature 438, 201 (2005).
  * Zhou  _et al._ [2006] S. Zhou, G.-H. Gweon, J. Graf, A. Fedorov, C. Spataru, R. Diehl, Y. Kopelevich, D.-H. Lee, S. G. Louie, and A. Lanzara, Nature physics 2, 595 (2006).
  * Neto  _et al._ [2009] A. C. Neto, F. Guinea, N. M. Peres, K. S. Novoselov, and A. K. Geim, Reviews of modern physics 81, 109 (2009).
  * Novoselov  _et al._ [2005] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. Grigorieva, S. Dubonos, and A. Firsov, Nature 438, 197 (2005).
  * Güttinger  _et al._ [2010] J. Güttinger, T. Frey, C. Stampfer, T. Ihn, and K. Ensslin, Physical review letters 105, 116801 (2010).
  * Saleem  _et al._ [2019] Y. Saleem, L. N. Baldo, A. Delgado, L. Szulakowska, and P. Hawrylak, Journal of Physics: Condensed Matter 31, 305503 (2019).
  * Güçlü  _et al._ [2014] A. D. Güçlü, P. Potasz, M. Korkusinski, and P. Hawrylak,  _Graphene Quantum Dots_ (Springer Verlag, 2014).
  * Güçlü  _et al._ [2011] A. Güçlü, P. Potasz, and P. Hawrylak, Physical Review B 84, 035425 (2011).
  * Güçlü  _et al._ [2009] A. Güçlü, P. Potasz, O. Voznyy, M. Korkusinski, and P. Hawrylak, Physical review letters 103, 246805 (2009).
  * Potasz  _et al._ [2010] P. Potasz, A. Güçlü, and P. Hawrylak, Physical Review B 81, 033403 (2010).
  * Potasz  _et al._ [2012] P. Potasz, A. Güçlü, A. Wójs, and P. Hawrylak, Physical Review B 85, 075431 (2012).
  * Ponomarenko  _et al._ [2008] L. A. Ponomarenko, F. Schedin, M. I. Katsnelson, R. Yang, E. W. Hill, K. S. Novoselov, and A. K. Geim, Science 320, 356 (2008).
  * Lu  _et al._ [2011] J. Lu, P. S. E. Yeo, C. K. Gan, P. Wu, and K. P. Loh, Nature nanotechnology 6, 247 (2011).
  * Ezawa [2010] M. Ezawa, Physical Review B 81, 201402 (2010).
  * Wunsch  _et al._ [2008] B. Wunsch, T. Stauber, and F. Guinea, Physical Review B 77, 035316 (2008).
  * Peleg  _et al._ [2007] O. Peleg, G. Bartal, B. Freedman, O. Manela, M. Segev, and D. N. Christodoulides, Physical review letters 98, 103901 (2007).
  * Gomes  _et al._ [2012] K. K. Gomes, W. Mar, W. Ko, F. Guinea, and H. C. Manoharan, Nature 483, 306 (2012).
  * Drost  _et al._ [2017] R. Drost, T. Ojanen, A. Harju, and P. Liljeroth, Nature Physics 13, 668 (2017).
  * Slot  _et al._ [2017] M. R. Slot, T. S. Gardenier, P. H. Jacobse, G. C. van Miert, S. N. Kempkes, S. J. Zevenhuizen, C. M. Smith, D. Vanmaekelbergh, and I. Swart, Nature physics 13, 672 (2017).
  * Lieb [1989] E. H. Lieb, Physical review letters 62, 1201 (1989).
  * Ezawa [2007] M. Ezawa, Physical Review B 76, 245415 (2007).
  * Fernández-Rossier and Palacios [2007] J. Fernández-Rossier and J. J. Palacios, Physical Review Letters 99, 177204 (2007).
  * Wang  _et al._ [2009] W. L. Wang, O. V. Yazyev, S. Meng, and E. Kaxiras, Physical review letters 102, 157201 (2009).
  * Pavliček  _et al._ [2017] N. Pavliček, A. Mistry, Z. Majzik, N. Moll, G. Meyer, D. J. Fox, and L. Gross, Nature nanotechnology 12, 308 (2017).
  * Mishra  _et al._ [2020] S. Mishra, D. Beyer, K. Eimre, S. Kezilebieke, R. Berger, O. Gröning, C. A. Pignedoli, K. Müllen, P. Liljeroth, P. Ruffieux, _et al._ , Nature nanotechnology 15, 22 (2020).
  * Su  _et al._ [2019] J. Su, M. Telychko, P. Hu, G. Macam, P. Mutombo, H. Zhang, Y. Bao, F. Cheng, Z.-Q. Huang, Z. Qiu, _et al._ , Science advances 5, eaav7717 (2019).
  * Wehling  _et al._ [2011] T. Wehling, E. Şaşıoğlu, C. Friedrich, A. Lichtenstein, M. Katsnelson, and S. Blügel, Physical Review Letters 106, 236805 (2011).
  * Wehling  _et al._ [2014] T. Wehling, A. M. Black-Schaffer, and A. V. Balatsky, Advances in Physics 63, 1 (2014).
  * Sorella and Tosatti [1992] S. Sorella and E. Tosatti, EPL (Europhysics Letters) 19, 699 (1992).
  * Voznyy  _et al._ [2011] O. Voznyy, A. D. Güçlü, P. Potasz, and P. Hawrylak, Physical Review B 83, 165417 (2011).
  * Ozfidan  _et al._ [2014] I. Ozfidan, M. Korkusinski, A. D. Güçlü, J. A. McGuire, and P. Hawrylak, Physical Review B 89, 085310 (2014).