# Quantum liquids with fractal branching patterns in Rydberg atom arrays

Tengzhou Zhang  Wilczek Quantum Center and Key Laboratory of Artificial
Structures and Quantum Control, School of Physics and Astronomy, Shanghai Jiao
Tong University, Shanghai 200240, China  Zi Cai<EMAIL_ADDRESS>Wilczek
Quantum Center and Key Laboratory of Artificial Structures and Quantum
Control, School of Physics and Astronomy, Shanghai Jiao Tong University,
Shanghai 200240, China  Shanghai Research Center for Quantum Sciences,
Shanghai 201315, China

###### Abstract

In this study, we propose an exotic quantum liquid state which does not order
at zero temperature in a Rydberg atom array with antiblockade mechanism. By
performing an unbiased large-scale quantum Monte Carlo simulation, we
investigate a minimal model with facilitated excitation in a clean system
without frustration, and find that the wavefunction of its ground state is
composed of configurations with fractal branching structures. This state
features a quasi-long-range order (correlation functions with algebraic decay)
as well as a heterogeneous structure of liquid and glass mixture. Thus it is
different from most well-established quantum phases of matter.

Introduction – The search for highly correlated quantum matters that do not
order even at zero temperature has been a focus across several disciplines of
quantum physics in the past decades. One of the best-studied examples in this
regards are the quantum spin liquids which feature long-range quantum
entanglement and topological excitationsWen (2017); Zhou et al. (2017).
Recently, Rydberg atom arrays has provided promising neutral-atom platforms
for quantum simulationWeimer et al. (2010); Semeghini et al. (2021); Ebadi et
al. (2021, 2022). A key element therein is the Rydberg blockade mechanism:
double excitations of adjacent atoms are strongly suppressedJaksch et al.
(2000); Lukin et al. (2001). The kinetic constraint, incorporating quantum
fluctuations, can give rise to intriguing quantum matters inSamajdar et al.
(2020, 2021); Verresen et al. (2021); Yue et al. (2021); Yan et al. (2023) and
outLesanovsky and Garrahan (2013); Lan et al. (2018); Bluvstein et al. (2021);
Wu et al. (2023) of equilibrium.

In this study, an exotic quantum liquid state other than the spin liquids is
proposed, which is based on a facilitation (antiblockade) mechanism with a
kinetic constraintAtes et al. (2007); Amthor et al. (2010): an atom can be
coherently excited if and only if just one of its neighbors has already been
in the Rydberg state. In another word, an excited atom facilitates the
excitation its neighbors. Unlike most work in this regards which mainly focus
on the dynamical properties Gärttner et al. (2013); Lesanovsky and Garrahan
(2014); Marcuzzi et al. (2017); Pérez-Espigares et al. (2017); Causer et al.
(2020); Helmrich et al. (2020); Ding et al. (2020); Liu et al. (2022), here we
study the ground state of a minimal model supporting such a kinetic constriant
using an unbiased quantum Monte Carlo (QMC) simulationSandvik (2003); Merali
et al. (2021). It is shown that even for such a simple disorder-free model
without frustration, the interplay between the kinetic constraint and quantum
fluctuation can give rise to an intriguing quantum liquid state, whose
wavefunction is composed of configurations with fractal branching structures.
This state does not order at zero temperature, and features an algebraic
decaying in its correlation function, distinguishing it from most well-
established quantum phases of matter.

Figure 1: (Color online) (a)Each atom has two internal states ($|g\rangle$ and
$|r\rangle$), $|g\rangle$ can be coherently excited to $|r\rangle$ by a laser
with Rabi frequency $\Omega$ and detuning $\Delta$. The transition becomes
resonant if the detuning compensates the NN interaction ($\Delta=V$);
(b)Examples of allowable (the lower two panels) and forbidden (the upper two
panels) transitions between different configurations under the kinetic
constraint. (c)A typical configuration generated by the importance sampling in
QMC simluation. (d)The area-radius relations ($\overline{M_{l}}\sim l$) for
three independently-sampled configurations in QMC simulation with $L=\beta=64$
and $\mu=0$.

Model and method – The system we studied is a Rydberg array where the atoms
are placed in a two-dimensional (2D) square lattice. Each atom possesses two
internal states: a ground state $|g\rangle$ can be coherently excited to a
Rydberg state $|r\rangle$ by a laser with Rabi frequency $\Omega$ and detuning
$\Delta$. The interaction strength between a pair of nearest-neighboring (NN)
excited atoms is $V$, and we ignore the longer-range interaction for
simplicity. The facilitated excitation emerges if the detuning happens to be
compensated by the NN interaction ($\Delta=V$), where the coherent transition
between $|g\rangle$ and $|r\rangle$ at a particular site becomes resonant only
if just one of its neighbors is already in state $|r\rangle$ (see Fig.1). This
anti-blockade mechanism has been realized in Rydberg aggregates, where the
strongly correlated growth induced by this facilitated excitation has been
explored experimentallySchempp et al. (2014); Urvoy et al. (2015).

Throughout this study, we ignore the incoherent transition and dissipation,
and the system can be considered as a closed system. The minimal model
describing the resulting quantum many-body system is rather simple, and its
Hamiltonian reads:

$H=\sum_{i}[-J\Xi_{\mathbf{i}}\sigma_{\mathbf{i}}^{x}-\mu n_{\mathbf{i}}],$
(1)

where $J$ and $\mu$ represent the amplitude of the coherent transition and the
effective detuning respectively.
$\hat{\sigma}^{x}_{\mathbf{i}}=|g\rangle_{\mathbf{i}}\langle
r|+|r\rangle_{\mathbf{i}}\langle g|$ and
$n_{\mathbf{i}}=|r\rangle_{\mathbf{i}}\langle r|$ is the density operator of
the Rydberg state on site $\mathbf{i}$. $\Xi_{\mathbf{i}}$ is the operator
projecting onto the subspace with only one of the four neighbors of site
$\mathbf{i}$ being in $|r\rangle$ state. $\Xi_{\mathbf{i}}$ can be expressed
in terms of the density operators of its neighbors as:

$\Xi_{i}=\sum_{\alpha}\big{[}n_{\mathbf{i}+\mathbf{e}_{\alpha}}\prod_{\gamma\neq\alpha}(1-n_{\mathbf{i}+\mathbf{e}_{\gamma}})\big{]}$
(2)

where $\alpha,\gamma=1\sim 4$, $\mathbf{e}_{1\sim 4}$ indicate the unit
vectors along $\pm x$ and $\pm y$ directions and $\mathbf{i}+\mathbf{e}_{1\sim
4}$ represent the four neighbors of site $\mathbf{i}$. It is easy to check
that $\Xi_{i}=1$ if and only if
$\sum_{\alpha}n_{\mathbf{i}+\mathbf{e}_{\alpha}}=1$, otherwise $\Xi_{i}=0$.

The kinetic constraint imposes additional conserved quantities on the system.
For each configuration, we define a cluster as a group of connected sites
occupied by $|r\rangle$ atoms, and the number of these clusters is conserved
since the spin flip terms under the kinetic constraint does not change the
connectivity of a cluster (e.g. they cannot break one cluster into two, see
Fig.1 b). This conserved quantity further divides the total Hilbert space into
different subspaces, each of which is characterized by its cluster number.
Here, we focus on the ground state of Hamiltonian.(1) within a particular
subspace with only one cluster, and general situations will be discussed in
the Supplementary materialSup .

We implement the QMC simulation with the stochastic series expansion (SSE)
algorithmSandvik (2003); Merali et al. (2021) on a system in a $L\times L$
square lattice with system size $N=L^{2}$. The periodical boundary condition
(PBC) is chosen and the inverse temperature is set to $\beta=L$ to scale to
the ground state in the thermodynamic limit. Notably, the SSE algorithm also
preserves the cluster number. In our QMC simulation, we focus on the single-
cluster subspace by choosing the initial state with only one “seed”: a
configuration with only one site is in $|r\rangle$ and all others are in
$|g\rangle$. This cluster will keep growing and ramifying in the QMC
simulation. It is straightforward to implement the kinetic constraint in the
SSE algorithm. Since this constraint as shown in Eq.(2) occurs in Fock basis
of the particle number representation, it does not cause negative sign problem
during the QMC sampling. As a result, the QMC simulation can give reliable
results for systems of sufficiently large size and low temperature.

Figure 2: (Color online) (a) For a given quantum many-body state
$|\Psi\rangle$, the relation between the averaged “area” of a cluster
($\langle\overline{M_{l}}\rangle_{\Psi}$) and it “radius ” ($l$) for systems
with a fixed $\mu=0$, but different size and inverse temperature ($L=\beta$).
(b)$\langle\overline{M_{l}}\rangle_{\Psi}$ v.s. $l$ with different $\mu$ and
$L=\beta=64$.

Fractal cluster dimension– We now focus on the properties of the ground state
within the single-cluster sector. The corresponding many-body wavefunction can
be expressed in terms of the Fock basis $|\sigma\rangle$ as:

$|\Psi\rangle=\sum_{\sigma}C_{\sigma}|\sigma\rangle$ (3)

Under the kinetic constraint, a “typical” configuration (a basis
$|\sigma\rangle$ with a relatively large coefficient) generated by the QMC
sampling exhibits loopless branching structure (e.g. Fig.1 c). Similar
structures widely exist in nature, such as lungs, coral, polymers, etc and
also appear in dynamical processes such as particle aggregation, crystal
growth and electrochemical decomposition, which have been shown to produce
fractal clustersVicsek (1992). The emergent dynamical fractal has also been
used to explain the anomalous noise observed in a clean spin ice systemHallen
et al. (2022).

There are different ways to measure the fractal dimension, here we adopt the
cluster-growing method, which has been used to characterize the fractal
dimension for the loopless branching structureVicsek (1992); Feder (1988). For
a given configuration $|\sigma\rangle$, all $|r\rangle$ sites form a connected
cluster, from which we randomly pick a site as origin. The chemical distance
between two sites within the cluster is defined as the length of the shortest
path connecting them in the clusterSup . Let $M_{l}$ be the total number of
the cluster sites whose chemical distance to the origin is not larger than $l$
($M_{l}$ being the “area” of the cluster surrounded by a “radius” $l$), then
the fractal dimension $d_{f}$ for such a single-cluster configuration $\sigma$
is defined via the area-radius relation as:

$\overline{M_{l}}\propto l^{d_{f}[\sigma]},\quad 1\ll l\ll L$ (4)

where $\overline{O}$ is the average over the results with random origin
choices. It is easy to check that $d_{f}[\sigma]$ takes integer values for
regular patternsSup , while for those typical fractal patterns,
$d_{f}[\sigma]$ is known to be non-integer (e.g. $d_{f}[\sigma]=1.67$ at a 2D
percolation critical pointBalankin (2015); Nakayama et al. (1994)).

For a quantum superposition state as in Eq.(3), all the configurations
$\sigma$ preserve the one-cluster structure. It is interesting to notice that
the $\overline{M_{l}}\sim l$ curves for different typical configurations in
the wavefunction.(3) are close to each other (see Fig.1 d), implying the
fractal cluster dimension $D_{f}[\Psi]$ is an intrinsic property of this
quantum state, and it can be defined similar as in Eq.(4):

$\langle\overline{M_{l}}\rangle_{\Psi}\propto l^{D_{f}[\Psi]},\quad 1\ll l\ll
L$ (5)

where $\langle O\rangle_{\Psi}=\sum_{\sigma}|C_{\sigma}|^{2}O[\sigma]$ is the
average over the quantum superposition state, which can be directly performed
via the importance sampling in QMC simulationSup . For the ground state within
the single-cluster sector, Fig.2 (a) indicates that the area-radius relation
in Eq.(5) can be numerically verified within an intermediate regime ($1\ll
l\ll L$), and this algebraic regime is extended to infinity in the
thermodynamic limit. In our simulation, we fit the data using Eq.(5) in the
regime $l\in[\frac{L}{4},\frac{3L}{4}]$. The fractal dimension can be
extracted as $D_{f}[\Psi]\approx 1.8$, which obviously deviates from $D=2$ as
shown in Fig.2 (a), indicating that it is indeed a fractal. The area-radius
relations for the ground states with different $\mu$ have also been plotted
Fig.2 (b), which seems to suggest that $D_{f}[\Psi]$ barely dependents on
$\mu$.

Figure 3: (Color online) (a) The correlation functions $C(\mathbf{r})$ along
the x-direction $\mathbf{r}=(r,0)$ for the system with different $\mu$
($L=\beta=64$), each of which exhibits algebraic decay with the same exponent
$\alpha$; (b) $C(\mathbf{r})$ in systems with a fixed $\mu=0$ but different
$L$ and $\beta$. Figure 4: (Color online) Structure factor along x-direction
$S(k_{x},0)$, which exhibits a peak at $k_{x}=k_{peak}$. The inset indicates
the dependence of the peak position on the average density of the $|r\rangle$
state. ($L=\beta=64$ and $\mu=0$). Figure 5: (Color online) The average
density of the $|r\rangle$ states $\rho$ as a function of the detuning $\mu$
with different system size and $L=\beta$.

Correlation functions with algebraic decay – One of the remarkable features of
fractal patterns is their scaling invariance, often exemplified by a power-law
behavior. To numerically validate this aspect, we calculate the correlation
function on the ground state:

$C(\mathbf{r})=\frac{1}{L^{2}}\sum_{\mathbf{i}}[\langle
n_{\mathbf{i}}n_{\mathbf{i}+\mathbf{r}}\rangle_{\Psi}-\langle
n_{\mathbf{i}}\rangle_{\Psi}\langle n_{j}\rangle_{\Psi}]$ (6)

as well as the corresponding structure factor:

$S(\mathbf{q})=\frac{1}{L^{2}}\sum_{\mathbf{r}}C(\mathbf{r})e^{i\mathbf{r}\cdot\mathbf{q}}$
(7)

In Fig.3, we plot $C(\mathbf{r})$ along x-direction, which exhibits an
algebraic decay accompanied by an oscillation. The exponent of the algebraic
decay ($\alpha=0.98(4)$) seems independent on $\mu$, and agrees with an
integer value $\alpha=1$ within the statistical error. Such a simple exponent
reminds us of the algebraic decay in the Coulomb phase of classical spin ice,
which also takes integer values ($\alpha=d$ with d being the system
dimension). However, our state is a pure quantum state thus is different from
the classical Coulomb phase in the spin ice.

A long-range order usually indicates a sharp peak in the structure factor,
with its height diverging in the thermodynamic limit. However, Fig. 4 reveals
no such sharp peaks in $S(\mathbf{k})$. Instead, $S(\mathbf{k})$ exhibits a
broad peak at the wavevector $k_{peak}$, which corresponds to the oscillation
period in $C(r)$. We also observe that $k_{peak}$ grows with the average
density $\rho=\frac{1}{L^{2}}\sum_{i}\langle n_{\mathbf{i}}\rangle_{\Psi}$
(see the inset of Fig. 4). One may wonder whether the emergence of such a
characteristic length $\frac{2\pi}{k_{peak}}$ contradicts the scaling
invariance feature of the fractal. However, as stated above, the scaling
invariance and fractal dimension manifest only on scales significantly larger
than the microscopic cut-off length (the lattice constant), since the emergent
period is also at a microscopic scale compatible to the lattice constant, it
does not contradicts the scaling invariance.

To distinguish this state from a Mott insulator, we calculate the average
density $\rho=\frac{1}{L^{2}}\sum_{i}\langle n_{\mathbf{i}}\rangle_{\Psi}$ as
a function of $\mu$. As illustrated in Fig. 5, there exists a critical value
at $\mu_{c}=J$, below which the energy cost of exciting an atom is so
substantial that the cluster cannot grow, thus $\rho\rightarrow 0$ in the
thermodynamic limit. Above $\mu_{c}$, $\rho$ increases monotonically with
$\mu$, and there is no plateau in the $\rho-\mu$ curve, indicating the absence
of a Mott insulator. The phase transition observed here is akin the absorbing
phase transition in non-equilibrium classical statistical modelsHenkel et al.
(2008).

Figure 6: (Color online) For a given QMC trajectory, the autocorrelation
functions $r_{i}(\tau)$ with the parameters $\mu=0$, $L=48$ and
$\beta=64J^{-1}$ (different curves indicate $r_{i}(\tau)$ on different sites.)

Autocorrelation function – The typical configurations exhibit branching
network structure as shown in Fig.1 (c), while the quantum fluctuations (the
spin flipped terms) create a superposition among them, leading to such an
exotic quantum state. For a given configuration, the kinetic constraint
dictates that the spin flipped terms can only operate on endpoints of the
cluster, while the interior points of the cluster are immune to quantum
fluctuations. This robustness against spin flipping indicates a glassy
behavior with slow relaxation dynamics for those interior points of the
network. For such a site, in spite of its robustness against single-spin flip,
it can slowly relax by transforming an interior point into a endpoint via a
sequence of spin flipping induced by quantum fluctuations.

To characterize the dynamical heterogeneity of these two types of sites
(endpoints and interior points), we calculate the on-site autocorrelation
function $r(\tau)$ in imaginary time $\tau$ on different system sites:

$r_{i}(\tau)=\frac{1}{R_{i}^{0}}\langle(n_{i}(0)-n)(n_{i}(\tau)-n)\rangle_{\tau}$
(8)

where the average $\langle O\rangle_{\tau}$ is performed over a single
trajectory of QMC simulation.
$R_{i}^{0}=\langle(n_{i}(0)-n)^{2}\rangle_{\tau}$, and this prefactor ensures
that for a given site where the quantum fluctuation is completely frozen,
$r_{i}(\tau\rightarrow\infty)\rightarrow 1$. Conversely, for those sites whose
spins are frequently flipped, $r_{i}(\tau)$ will decay to zero.

The autocorrelation functions on different sites are illustrated in Fig.6,
from which we can find that $r_{i}(\tau)$ decays with $\tau$ for the majority
of system sites, while for other sites which are randomly distribued, quantum
fluctuations remain frozen up to the longest imaginary time scale
($\tau=\beta/2$) in our simulation, indicating a glassy behavior. Different
from the spin glass, our model is free from disorder, thus the glassy
mechanism in such a translational invariant system can be attribued to the
kinetic constraints, akin to those observed in genuine glass systems. Those
flippable sites enhance the quantum coherence that is responsible for the
distinctive quasi-long-range order (algebraic instead of exponential decay of
$C(\mathbf{r})$), which makes it qualitatively differ from the conventional
glass phases with only short-range order. Such a heterogeneous structure of
liquid and glass mixture is reminiscent of an intriguing classical spin liquid
state dubbed “spin slush”, which was first proposed by Rau et.al in a
classical spin ice systemRau and Gingras (2016).

Discussion – Throughout this paper, we focus on a specific subspace of Hilbert
space with only one cluster. Although our model enforces a stringent kinetic
constraint, real experimental systems implement this constraint through an
energy penalty mechanism, which allows for potential violations in principle.
Such a constraint-violating perturbation may induce quantum tunneling between
subspaces characterized by different cluster numbers. In another word, this
perturbation transforms the previously stable ground state with a single
cluster, as discussed earlier, into a metastable state. The lifetime of this
metastable state is inversely proportional to the strength of the
perturbation. Consequently, as long as the lifetime of this proposed state
greatly exceeds the typical time scale of a practical Rydberg atom system, it
remains observable, permitting experimental investigation into its properties.

Conclusion and outlook – In conclusion, we proposed an exotic fractal quantum
states, which is different from most well-established quantum phases of
matter. It is neither an ordered phase, nor a Mott insulator. The presence of
quasi-long-range order distinguishes it from conventional spinMezard et al.
(1986) and Bose glassesGiamarchi and Schulz (1988); Fisher et al. (1989). This
phase also stands apart from the U(1) quantum spin liquidHermele et al. (2004)
or the Coulomb phase in the classical spin iceHenley (2010); C.Castelnovo et
al. (2012) due to the absence of emergent gauge field. Compared to other
states with fractal structures (e.g. the critical states in
percolationNakayama et al. (1994); Balankin (2015) or quantum many-body
systemsZhou et al. (2021)), the proposed state is not intrinsically tied to
criticality, thus it exists in a wide regime instead of at a particular point
in the parameter space.

One of the unsolved problems worthy of further study is the physical origin of
the algebraic decay and the associated simple power exponent. In the Coulomb
phase of spin ice, the algebraic decay can be understood as a consequence of
the emergent $U(1)$ gauge field due to the constraint imposed by the ice
ruleHenley (2010), one may wonder whether there exists a similar mechanism
associated with kinetic constraint in our model. Another unexplored avenue is
the effect of thermal fluctuation, for instance, whether this quasi-long-range
order can persist at low temperature analogous to the two-dimensional XY
modelKosterlitz and Thouless (1973), or is not robust against any thermal
fluctuation like the Coulomb phase in the spin ice. A related question is
what’s the nature of the elementary excitations for such a state, do they
manifest as collective excitations or an individual ones?

Acknowledgments.— This work is supported by the National Key Research and
Development Program of China (Grant No. 2020YFA0309000), NSFC of China (Grant
No.12174251), Natural Science Foundation of Shanghai (Grant No.22ZR142830),
Shanghai Municipal Science and Technology Major Project (Grant
No.2019SHZDZX01).

## References

  * Wen (2017) X.-G. Wen, Rev. Mod. Phys. 89, 041004 (2017).
  * Zhou et al. (2017) Y. Zhou, K. Kanoda, and T.-K. Ng, Rev. Mod. Phys. 89, 025003 (2017).
  * Weimer et al. (2010) H. Weimer, M. Muller, I. Lesanovsky, P. Zoller, and H. P. Buchler, Nature Phys. 6, 382 (2010).
  * Semeghini et al. (2021) G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, et al., Science 374, 1242 (2021).
  * Ebadi et al. (2021) S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, et al., Nature 595, 227 (2021).
  * Ebadi et al. (2022) S. Ebadi, A. Keesling, M. Cain, T. T. Wang, H. Levine, D. Bluvstein, G. Semeghini, A. Omran, J. Liu, R. Samajdar, et al., Science 376, 1209 (2022).
  * Jaksch et al. (2000) D. Jaksch, J. I. Cirac, P. Zoller, S. L. Rolston, R. Côté, and M. D. Lukin, Phys. Rev. Lett. 85, 2208 (2000).
  * Lukin et al. (2001) M. D. Lukin, M. Fleischhauer, R. Cote, L. M. Duan, D. Jaksch, J. I. Cirac, and P. Zoller, Phys. Rev. Lett. 87, 037901 (2001).
  * Samajdar et al. (2020) R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev, Phys. Rev. Lett. 124, 103601 (2020).
  * Samajdar et al. (2021) R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev, Proc. Nat. Acad. Sci. 118, e2015785118 (2021).
  * Verresen et al. (2021) R. Verresen, M. D. Lukin, and A. Vishwanath, Phys. Rev. X 11, 031005 (2021).
  * Yue et al. (2021) M. Yue, Z. Wang, B. Mukherjee, and Z. Cai, Phys. Rev. B 103, L201113 (2021).
  * Yan et al. (2023) Z. Yan, Y.-C. Wang, R. Samajdar, S. Sachdev, and Z. Y. Meng, Phys. Rev. Lett. 130, 206501 (2023).
  * Lesanovsky and Garrahan (2013) I. Lesanovsky and J. P. Garrahan, Phys. Rev. Lett. 111, 215305 (2013).
  * Lan et al. (2018) Z. Lan, M. van Horssen, S. Powell, and J. P. Garrahan, Phys. Rev. Lett. 121, 040603 (2018).
  * Bluvstein et al. (2021) D. Bluvstein, A. Omran, H. Levine, A. Keesling, G. Semeghini, S. Ebadi, T. T. Wang, A. A. Michailidis, N. Maskara, W. W. Ho, et al., Science 371, 1355 (2021).
  * Wu et al. (2023) X. Wu, Z. Wang, F. Yang, R. Gao, C. Liang, M. Khoon Tey, X. Li, T. Pohl, and L. You, arXiv e-prints arXiv:2305.20070 (2023), eprint 2305.20070.
  * Ates et al. (2007) C. Ates, T. Pohl, T. Pattard, and J. M. Rost, Phys. Rev. Lett. 98, 023002 (2007).
  * Amthor et al. (2010) T. Amthor, C. Giese, C. S. Hofmann, and M. Weidemüller, Phys. Rev. Lett. 104, 013001 (2010).
  * Gärttner et al. (2013) M. Gärttner, K. P. Heeg, T. Gasenzer, and J. Evers, Phys. Rev. A 88, 043410 (2013).
  * Lesanovsky and Garrahan (2014) I. Lesanovsky and J. P. Garrahan, Phys. Rev. A 90, 011603 (2014).
  * Marcuzzi et al. (2017) M. Marcuzzi, J. c. v. Minář, D. Barredo, S. de Léséleuc, H. Labuhn, T. Lahaye, A. Browaeys, E. Levi, and I. Lesanovsky, Phys. Rev. Lett. 118, 063606 (2017).
  * Pérez-Espigares et al. (2017) C. Pérez-Espigares, M. Marcuzzi, R. Gutiérrez, and I. Lesanovsky, Phys. Rev. Lett. 119, 140401 (2017).
  * Causer et al. (2020) L. Causer, I. Lesanovsky, M. C. Bañuls, and J. P. Garrahan, Phys. Rev. E 102, 052132 (2020).
  * Helmrich et al. (2020) S. Helmrich, A. Arias, G. Lochead, T. M. Wintermantel, M. Buchhold, S. Diehl, and S. Whitlock, Nature 577, 481 (2020).
  * Ding et al. (2020) D.-S. Ding, H. Busche, B.-S. Shi, G.-C. Guo, and C. S. Adams, Phys. Rev. X 10, 021023 (2020).
  * Liu et al. (2022) F. Liu, Z.-C. Yang, P. Bienias, T. Iadecola, and A. V. Gorshkov, Phys. Rev. Lett. 128, 013603 (2022).
  * Sandvik (2003) A. W. Sandvik, Phys. Rev. E 68, 056701 (2003).
  * Merali et al. (2021) E. Merali, I. J. S. De Vlugt, and R. G. Melko, arXiv e-prints arXiv:2107.00766 (2021), eprint 2107.00766.
  * Schempp et al. (2014) H. Schempp, G. Günter, M. Robert-de Saint-Vincent, C. S. Hofmann, D. Breyel, A. Komnik, D. W. Schönleber, M. Gärttner, J. Evers, S. Whitlock, et al., Phys. Rev. Lett. 112, 013002 (2014).
  * Urvoy et al. (2015) A. Urvoy, F. Ripka, I. Lesanovsky, D. Booth, J. P. Shaffer, T. Pfau, and R. Löw, Phys. Rev. Lett. 114, 203002 (2015).
  * (32) See the supplementary material for an elaboration of the fractal cluster dimension, as well as a comparision between the ground states within the single-cluster subspace and subspaces with other cluster numbers.
  * Vicsek (1992) T. Vicsek, _Fractal Growth Phenomena, Part IV_ (World Scientific, Singapore, 1992).
  * Hallen et al. (2022) J. N. Hallen, S. A. Grigera, D. A. Tennant, C. Castelnovo, and R. Moessner, Science 378, 1218 (2022).
  * Feder (1988) J. Feder, _Fractals_ (Springer Science, New York, 1988).
  * Balankin (2015) A. S. Balankin, Phys. Rev. E 92, 062146 (2015).
  * Nakayama et al. (1994) T. Nakayama, K. Yakubo, and R. L. Orbach, Rev. Mod. Phys. 66, 381 (1994).
  * Henkel et al. (2008) M. Henkel, H. Hinrichsen, and S. Lubeck, _Non-Equilibrium Phase Transitions, Volume 1: Absorbing Phase Transitions_ ( Springer Dordrecht, 2008).
  * Rau and Gingras (2016) J. G. Rau and M. J. P. Gingras, Nat. Comms 7, 12234 (2016).
  * Mezard et al. (1986) M. Mezard, G. Parisi, and M. Virasoro, _Spin Glass Theory and Beyond: An Introduction to the Replica Method and Its Applications_ ( World Scientific, 1986).
  * Giamarchi and Schulz (1988) T. Giamarchi and H. J. Schulz, Phys. Rev. B 37, 325 (1988).
  * Fisher et al. (1989) M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Phys. Rev. B 40, 546 (1989).
  * Hermele et al. (2004) M. Hermele, T. Senthil, M. P. A. Fisher, P. A. Lee, N. Nagaosa, and X.-G. Wen, Phys. Rev. B 70, 214437 (2004).
  * Henley (2010) C. Henley, Annu. Rev. Condens. Mat. Phys. 1, 179 (2010).
  * C.Castelnovo et al. (2012) C.Castelnovo, R.Moessner, and S. Sondhi, Annu. Rev. Condens. Mat. Phys. 3, 35 (2012).
  * Zhou et al. (2021) Z. Zhou, X.-F. Zhang, F. Pollmann, and Y. You, arXiv e-prints arXiv:2105.05851 (2021), eprint 2105.05851.
  * Kosterlitz and Thouless (1973) J. M. Kosterlitz and D. J. Thouless, Journal of Physics C: Solid State Physics 6, 1181 (1973).

Supplemental material

Figure 7: (Color online) Definition of the chemical distance for a regular
pattern (left) and branching cluster pattern (right) Figure 8: (Color online)
(a)-(c) Three typical configurations independently generated by the QMC
sampling. (d) The area-radius relation for these three independent
configurations. The parameters are chosen as $L=\beta=64$, $\mu=0$.

##  I Fractal cluster dimension

We first elaborate the definition of the fractal cluster dimension via the
area-radius relation, then explain the details of generating the typical
configurations via the importance sampling in QMC simulations.

###  I.1 Definition

For a given configuration $|\sigma\rangle$ with a single connected cluster, we
randomly choose a site as origin, then the chemical distance $l$ between the
origin and another site is defined as the length of the shortest path in the
cluster connecting these two sites (see Fig.7 for an example).

Let $M_{l}$ be the total number of $|r\rangle$ sites whose chemical distance
to the origin is not larger than $l$, then the fractal dimension $d_{f}$ of in
its cluster configuration $\sigma$ is defined based on the relationship
between $M_{l}$ and $l$ as:

$\overline{M_{l}}\propto l^{d_{f}[\sigma]},\quad 1\ll l\ll L$ (9)

where $\overline{O}$ is the average over the results staring from different
choices of the origin within the same cluster.

Notice that the area-radius relation in Eq.(9) only holds in an intermediate
regime of the length scale, which is much larger than the microscopic cut-off
length (the lattice constant), but significantly smaller than the lattice
length L. The short distance behavior depends on the microscopic details of
the system, thus is not universal in general. On the other hand, the area of a
cluster starts to be saturated once its radius reaches the boundary of the
system, as shown in Fig. 8 (d). In our simulation, we extract the fractal
dimension by fitting the area-radius curves using Eq.(9) within the regime
$\frac{L}{4}\leq l\leq\frac{3L}{4}$.

For a regular patterns with Euclid geometry, it is easy to check that
$d_{f}[\sigma]$ takes integer values (e.g. $d_{f}[\sigma]=2$ for the
configuration as shown in Fig.7 a). While for those typical fractal patterns,
$d_{f}[\sigma]$ is known to be non-integer. A typical example is the
configurations at a 2D percolation critical point, where
$d_{f}[\sigma]=1.67$Balankin (2015); Nakayama et al. (1994).

###  I.2 Typical configurations generated by importance sampling and their
fractal cluster dimensions

In general, a quantum many-body wavefunction can be expanded in terms of the
Fock basis $|\sigma\rangle$ in the particle number representation as:

$\Psi=\sum_{\sigma}C_{\sigma}|\sigma\rangle$ (10)

for an operator $\hat{O}$ which is diagonal in terms of Fock basis, its
average value can be written as:

$\langle\hat{O}\rangle_{\Psi}=\sum_{\sigma}|C_{\sigma}|^{2}\langle\sigma|\hat{O}|\sigma\rangle$
(11)

Usually, the Hilbert space dimension of a quantum many-body system grows
exponentially with the system size, which makes it impossible to perform the
summation over all the basis. The spirt of Monte Carlos method lies in the
importance sampling: instead of averaging over all configurations, one only
needs to perform average over those “important” configurations, whose number
is much smaller than the total Hilbert space dimension. In Monte Carlo
simulation, the average value of a diagonal operator can be approximated as:

$\langle\hat{O}\rangle_{\Psi}\approx\overline{\langle\sigma|\hat{O}|\sigma\rangle}$
(12)

where the average $\overline{O}$ is over those representative configurations,
which are randomly selected according to their probabilities (weights) in the
wavefunction or partition function.

For certain quantum models (e.g. the model studied in the main text), a
$d$-dimensional quantum thermal ensemble is equivalently described by the path
integral representing the partition function in $d+1$ Euclidean space-time
with the time being imaginary. In the quantum Monte Carlo simulation, the
space-time configurations are sampled according to their weight in the
partition function. To calculate the average value of a diagonal operator, one
can choose the spacial configurations at different (imaginary) time slices,
calculate the corresponding expectation values of the operator, then perform
average over all of them.

The fractal cluster dimension is actually a diagonal operator since it only
involves density operators, thus in principle one can first calculate the
fractal cluster dimension for different configurations generated by QMC
sampling and perform average over them. In Fig.8, we plot three typical
configurations independently generated by the QMC sampling, and the relations
between $M_{l}$ and $l$ for these three configuration have been shown in Fig.8
(d), from which we can see that for a sufficiently large system, these curves
agree with each other within the statistical error bar. This result indicates
that the fractal dimension is an intrinsic property for the ground state
wavefunction, instead of a consequence of an average over different integer
values.

Figure 9: (Color online) The correlation functions along the x-direction for
the ground states within subspaces with different cluster numbers
$\mathcal{N}$ with parameter $L=\beta=64$ and $\mu=0$.

##  II Ground states in subsections of Hilbert space with cluster numbers
other than one

The number of the cluster is a conserved quantity due to the kinetic
constraint. Throughout the main text, we only focus on the ground state in
single-cluster subspace. Here we study the cases with other cluster numbers to
check how the results depend on this conserved quantity.

###  II.1 The ground states with a few clusters

We first focus on the cases with a few clusters (the cluster number
$\mathcal{N}$ is an integer that is independent on the system size). The SSE
update algorithm preserves the cluster number during the QMC simulation, so we
choose the initial state as an almost empty state with only $\mathcal{N}$
occupied $|r\rangle$ sites ($\mathcal{N}$ seeds) whose positions are randomly
distributed. Staring from these initial states, we calculate the correlation
function $C(r)$ and compare them to the correlation function of the ground
state within the single-cluster sector. We plot $C(r)$ for the cases with
$\mathcal{N}=1,2,3$ in Fig.9, from which we can see that in all the three
cases, $C(r)$ in the corresponding ground state decays algebraically with
almost the same power exponent within the statistical error bar. These results
indicates that in the thermodynamic limit, the ground states with a few
clusters share the same nature.

Figure 10: (Color online) (a) The correlation function $C(\mathbf{r})$ in the
ground state in the overall Hilbert space with parameter $L=\beta=64$ and
$\mu=0$. (b)The same correlation function along the diagonal ($x=y=r$) and
anti-diagonal ($x=-y=r$) directions.

###  II.2 The ground state in the overall Hilbert space

Even though the ground states in the subspaces with a few clusters are similar
to each other, further increasing the cluster number $\mathcal{N}$ (e.g.
$L\rightarrow\infty$, while keep $\frac{\mathcal{N}}{L^{2}}$ being a constant)
might qualitatively change the nature of the ground state. To verify this
point, we consider the ground state within the overall Hilbert space.

Since $\mathcal{N}$ is a conserved quantity, in principle, to find the ground
state in the overall Hilbert space, one needs to compare the energies of the
ground states in all the subspaces, which is difficult due to the typical
number of $\mathcal{N}\sim O(L^{2})$ is large thus it’s inefficient to explore
all the subspaces. To avoid this complexity, we introduce a small perturbation
into the original Hamiltonian, which slightly violates the kinetic constraint,
thus leads to the tunneling between the subspaces characterized by different
$\mathcal{N}$. In another word, $\mathcal{N}$ is no longer a conserve quantity
in the QMC simulation, and the system will equilibrate towards the ground
state of the overall Hilbert space for a sufficiently low temperature in our
QMC simulation.

We focus on the case with $\mu=0$ for simplicity, and the perturbed
Hamiltonian reads:

$H=J\sum_{i}[(1-\epsilon)\Xi_{\mathbf{i}}\sigma_{\mathbf{i}}^{x}+\epsilon(\mathbb{I}-\Xi_{\mathbf{i}})\sigma_{\mathbf{i}}^{x}$
(13)

where $\mathbb{I}$ is the identity operator, and the last term in Eq.(13)
indicates the constraint-violating perturbation. The strength of this
perturbation $\epsilon$ is supposed to be sufficiently small that it does not
qualitatively change the nature of the ground state of the original $H_{0}$,
but not too small so we can search the overall ground state within a
reasonable time scale in our QMC simulation.

We first fix $\epsilon=0.01J$, and perform the QMC simulation for a
sufficiently long time until the system falls into the subspace with the
lowest ground state energy, then we slowly quench $\epsilon$ to zero, to make
sure that the system finally converge to the ground state in this subspace
without perturbation. We find that this subspace satisfies
$\mathcal{N}/L^{2}=1/8$ for the system with $\mu=0$. We focus on the
correlation function $C(\mathbf{r})$ of the overall ground state. The 2D
correlation function is plotted in Fig.10, from which we can find that the
ground state in this case exhibits a charge-density-wave (CDW) ordering. It is
interesting to notice that this CDW phase does not only spontaneously breaks
the translational symmetry as it is supposed to be, but also breaks the
rotational symmetry from a $Z_{4}$ symmetry group into a $Z_{2}$ one. This
spontaneous $Z_{2}$ rotational symmetry breaking can be seen more clearly by
comparing the $C(\mathbf{r})$ along the diagonal $(x=y)$ and anti-diagonal
$(x=-y)$ directions. As shown in Fig.10, $C(\mathbf{r})$ along both directions
exhibit persistent oscillation, while their amplitude are different.


# Tailoring coercive fields and the Curie temperature via proximity coupling
in WSe2/Fe3GeTe2 van der Waals heterostructures

Guodong Ma1,4, Renjun Du1,4, Fuzhuo Lian1,2, Song Bao1, Zijing Guo1, Xiaofan
Cai1, Jingkuan Xiao1, Yaqing Han1, Di Zhang1, Siqi Jiang1, Jiabei Huang1,
Xinglong Wu1, Alexander S. Mayorov1,∗, Jinsheng Wen1,3, Lei Wang1,3,∗ and
Geliang Yu1,3,∗ 1 National Laboratory of Solid State Microstructures, School
of Physics, Nanjing University, Nanjing 210093, China  2 Chongqing 2D
Materials Institute, Liangjiang New Area, Chongqing, 400714, China  3
Collaborative Innovation Center of Advanced Microstructures, Nanjing
University, Nanjing 210093, China  4 These authors contributed equally to this
work.  ∗ Authors to whom any correspondence should be addressed.
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Hybrid structures consisting of two-dimensional (2D) magnets and
semiconductors have exhibited extensive functionalities in spintronics and
opto-spintronics. In this work, we have fabricated WSe2/Fe3GeTe2 van der Waals
(vdW) heterostructures and investigated the proximity effects on 2D magnetism.
Through reflective magnetic circular dichroism (RMCD), we have observed a
temperature-dependent modulation of magnetic order in the heterostructure. For
temperatures above $40\mathord{\thinspace\rm K}$, WSe2-covered Fe3GeTe2
exhibits a larger coercive field than that observed in bare Fe3GeTe2,
accompanied by a noticeable enhancement of the Curie temperature by
$21\mathord{\thinspace\rm K}$. This strengthening suggests an increase in
magnetic anisotropy in the interfacial Fe3GeTe2 layer, which can be attributed
to the spin-orbit coupling (SOC) proximity effect induced by the adjacent WSe2
layers. However, at much lower temperatures ($T<20\mathord{\thinspace\rm K}$),
a non-monotonic modification of the coercive field is observed, showing both
reduction and enhancement, which depends on the thickness of the WSe2 and
Fe3GeTe2 layers. Moreover, an unconventional two-step magnetization process
emerges in the heterostructure, indicating the short-range nature of SOC
proximity effects. Our findings revealing proximity effects on 2D magnetism
may shed light on the design of future spintronic and memory devices based on
2D magnetic heterostructures.

††:

, ,

  *

  * Keywords: vdW heterostructure, 2D magnetism, spin-orbit coupling, proximity effect

##  1 Introduction

The integration of magnetic materials with semiconductors has been extensively
explored in spintronics through active manipulations of the spin degrees of
freedom[1, 2]. This exploration has led the investigation of various
prototypical devices, such as spin field-effect transistors[3] and magnetic
tunnel junctions[4, 5], to unlock fascinating functionalities within these
hybrid structures. The interfaces between magnets and semiconductors play a
critical role in determining these novel properties. Recently, the utilization
of two-dimensional materials to form van der Waals (vdW) magnetic
heterostructures has provided a new platform for the development of ultrathin
spintronic devices[6, 7, 8, 9], owing to their high-quality crystalline
structures and bond-free interfaces[10, 11, 12].

Incorporating 2D magnets with 2D semiconductors in a vdW heterostructure
enables a wide range of intriguing phenomena through magnetic proximity
effects. For instance, WSe2/CrI3 heterostructures demonstrate a spontaneous
emergence of both valley Zeeman splitting and polarization[13, 14]. The former
is highly sensitive to the overall magnetization while the latter is dominated
by the interfacial layer of CrI3 and can be controlled through electrical
gating[15, 16]. Additionally, the magnetic proximity effect in MoSe2/CrBr3
exhibits a dependency on the charge state[17]. The examples mentioned above
focus on the impact of magnetic layers on nonmagnetic semiconductors.
Conversely, there is growing interest in modulating 2D magnets through
adjacent nonmagnetic layers[18, 19, 20, 21]. In 2D magnets, spin-orbit
coupling (SOC) from heavy elements is crucial for magnetic anisotropy, which
sustains long-range magnetic order even at the monolayer limit[22, 23]. It has
been proposed that combining 2D magnets with materials that exhibit strong SOC
can effectively enhance the magnetic anisotropy[24, 25, 26, 27, 21], resulting
in a modification of the magnetic order. However, experimental studies are
few, which require extensive efforts.

Here, we investigate the influence of proximity effects on a magnetic layer in
vdW heterostructures, consisting of a 2D itinerant magnet Fe3GeTe2 (FGT) and a
2D semiconductor WSe2 with strong SOC[28, 7]. Through reflective magnetic
circular dichroism (RMCD) studies, we observe a pronounced enhancement of the
coercive field and Curie temperature in the WSe2/Fe3GeTe2 heterostructure
compared with bare FGT, arising from SOC occurring at the interface. However,
a non-monotonic tuning of coercive fields, including both reduction and
enhancement, emerges in WSe2/Fe3GeTe2 regions at much lower temperatures
($<20\mathord{\thinspace\rm K}$), depending on the thickness of WSe2 and
Fe3GeTe2 layers. Moreover, an unconventional magnetic stratification is
observed in the heterostructure with monolayer WSe2, leading to an obvious
two-step magnetization reversal, which reveals the short-range nature of SOC
proximity effects.

##  2 Results and discussions

Figure 1a shows the schematic of the WSe2/Fe3GeTe2 heterostructure, assembled
by an all-dry viscoelastic stamping method[29]. In our experiments, we
exfoliated Fe3GeTe2 flakes onto a Si substrate with $285\mathord{\thinspace\rm
nm}$ SiO2. The thin WSe2 flakes were exfoliated on polydimethylsiloxane
(PDMS). Then we transferred the WSe2 onto the Fe3GeTe2 flake via a high-
precision transfer setup to construct the final heterostructure. The
fabrication was carried out inside a glovebox with an argon atmosphere to
avoid the reaction with oxygen and moisture. Figure 1b shows an optical image
of the heterostructure H1 with few-layer FGT and bilayer WSe2. The thickness
of the FGT flake is $6.2\mathord{\thinspace\rm nm}$ (about 7 layers)
determined by atomic force microscopy (AFM) (see the inset of Fig. 1b). The
room-temperature photoluminescence (PL) spectra of the WSe2 flakes (Fig. 1c)
exhibit a distinct dependence on thickness, in which a strong single peak
centered at $1.67\mathord{\thinspace\rm eV}$ emerges in the monolayer
WSe2[30]. This allows us to identify the layer numbers of the corresponding
WSe2 flakes in heterostructures.

Figure 1: (a) Schematic of the RMCD measurement configuration and crystal
structure of the WSe2/Fe3GeTe2 heterostructure. (b) Optical image of bilayer
WSe2 on Fe3GeTe2. The inset shows the profile extracted along the red line in
(b). (c) PL spectra of the WSe2 flakes with different thicknesses. (d,e)
Comparison of normalized RMCD signal for bare FGT and WSe2-covered FGT in H1
at (d) $110\mathord{\thinspace\rm K}$ and (e) $4\mathord{\thinspace\rm K}$,
respectively. The grey dashed lines mark the positions of coercive fields of
bare FGT.

To investigate the WSe2-induced proximity effect on magnetic order, we
performed a series of RMCD measurements by an optical cryostat. A
$633\mathord{\thinspace\rm nm}$ HeNe laser with excitation power of
$10\mathord{\thinspace\rm\upmu W}$ was used and focused onto the sample (Fig.
1a) with a spot size of about $1\mathord{\thinspace\rm\upmu m}$. Figures 1d,e
show a comparison of the RMCD signals of WSe2-covered FGT (WSe2+FGT) region
and the bare FGT region from the same heterostructure H1 (bilayer WSe2/6.2 nm
FGT), at the temperature of $110\mathord{\thinspace\rm K}$ and
$4\mathord{\thinspace\rm K}$, respectively. The typical rectangular hysteresis
loops emerge in both regions. Notably, we observe an opposite effect on the
coercive field $H_{\mathrm{c}}$ at different temperatures. At
$110\mathord{\thinspace\rm K}$, the loop in the WSe2+FGT region (red curve at
the bottom) is broader than that measured in bare FGT (blue curve at the top),
as shown in Fig. 1d, while $H_{\mathrm{c}}$ exhibits a reduction in the
WSe2+FGT region at $4\mathord{\thinspace\rm K}$.

Figures 2a,b present the temperature dependence of magnetic hysteresis loops
in sample H1, compared between the bare FGT and WSe2+FGT regions. As
temperature increases, the hysteresis loop shrinks rapidly in both regions due
to the stronger thermal fluctuation. It disappears at
$135\mathord{\thinspace\rm K}$ in bare FGT, implying a transition from the
ferromagnetic phase to paramagnetic phase, but retains a rectangular loop in
the WSe2+FGT region under the same conditions. Figure 2c shows the remanent
RMCD signal as a function of temperature, in which we also compare the results
for bare FGT (blue squares) and WSe2+FGT (red circles). The corresponding
dashed curves are plotted by the critical power-law form
$\alpha(1-T/T_{\mathrm{c}})^{\beta}$, by which we extract the Curie
temperature $T_{\mathrm{c}}$ and the critical exponent $\beta$. The
$T_{\mathrm{c}}$ of WSe2-covered FGT is increased by $21\mathord{\thinspace\rm
K}$ due to the WSe2 proximity effect ($135\mathord{\thinspace\rm K}$ in bare
FGT and $156\mathord{\thinspace\rm K}$ for WSe2+FGT). The value of $\beta$ we
obtain is $0.19\pm 0.01$, consistent with the previous report on a few layers
of FGT[23]. Figure 2d displays the corresponding extracted $H_{\mathrm{c}}$ as
a function of temperature for both two regions. The enhancement of
$H_{\mathrm{c}}$ can be observed in the same temperature range. The
significant increase in $T_{\mathrm{c}}$ and the simultaneous enhancement of
$H_{\mathrm{c}}$ confirm the improved magnetism in FGT, which can be
attributed to the proximity effect from the adjacent WSe2.

In contrast to the enhancement of coercive fields at higher temperatures
($T>40\mathord{\thinspace\rm K}$), as shown in Fig. 2, WSe2-covered FGT
demonstrates a distinct reduction of $H_{\mathrm{c}}$ ($\sim
330\mathord{\thinspace\rm Oe}$) at $4\mathord{\thinspace\rm K}$ (see Fig. 3a).
We measured the layer dependence of WSe2 on $H_{\mathrm{c}}$, as shown in Fig.
3a, demonstrating a discernible reduction in $H_{\mathrm{c}}$ with the
increase in layer numbers. Figures 3b,c depict the difference of coercive
fields ($\Delta
H_{\mathrm{c}}=H_{\mathrm{c}}^{\mathrm{W}}-H_{\mathrm{c}}^{\mathrm{FGT}}$) as
a function of temperature in both 2L-WSe2/FGT and 5L-WSe2/FGT regions, where
$H_{\mathrm{c}}^{\mathrm{W}}$ and $H_{\mathrm{c}}^{\mathrm{FGT}}$ represent
the extracted coercive field in WSe2-covered FGT and bare FGT, respectively.
As shown in Fig. 3b, the negative $\Delta H_{\mathrm{c}}$ associated with the
suppression of magnetic anisotropy is evident in both 2L- and 5L-WSe2 covered
regions at low temperatures, persisting as the temperature increases until
$20\mathord{\thinspace\rm K}$. $H_{\mathrm{c}}$ in 5L-WSe2+FGT is smaller than
that in 2L-WSe2+FGT, suggesting a stronger suppression of magnetic anisotropy.
Additionally, as the temperature rises above $12\mathord{\thinspace\rm K}$,
$H_{\mathrm{c}}$ of the 2-layer (red circles) and the 5-layer (yellow
triangles) WSe2+FGT are almost equal within the margin of error. Figure 3c
presents the results for the overall temperature range below $T_{\mathrm{c}}$.
The sign of $\Delta H_{\mathrm{c}}$ changes from negative (blue panel) to
positive (red panel) at approximately $20\mathord{\thinspace\rm K}$, and the
regions with different thicknesses of WSe2 show no pronounced difference at
elevated temperatures.

Figure 2: Proximity-induced enhancing PMA in H1. (a,b) Comparison of
temperature-dependent RMCD signal for bare FGT (a) and WSe2-covered FGT (b).
(c) Remanent RMCD signal as a function of temperature for different regions in
H1. The dashed fitting curves are plotted by the critical power-law form
$\alpha(1-T/T_{\mathrm{c}})^{\beta}$. (d) Extracted coercive field
$H_{\mathrm{c}}$ as a function of temperature for different regions in H1.
Figure 3: (a) Normalized RMCD signal measured at different regions: blue line,
bare FGT; red line, 2L WSe2/FGT; yellow line, 5L WSe2/FGT. (b)
$H_{\mathrm{c}}$ difference, $\Delta
H_{\mathrm{c}}=H_{\mathrm{c}}^{\mathrm{W}}-H_{\mathrm{c}}^{\mathrm{FGT}}$, as
a function of temperature, where $H_{\mathrm{c}}^{\mathrm{W}}$ and
$H_{\mathrm{c}}^{\mathrm{FGT}}$ represent the coercive field in regions of
WSe2-covered FGT and bare FGT, respectively. (c) $\Delta H_{\mathrm{c}}$ as a
function of temperature with the full temperature range. Two distinct
background colors represent the sign of $H_{\mathrm{c}}$ (blue, negative; red,
positive, respectively).

Thus far, our observations have revealed the coexistence of the proximity-
induced enhancement and suppression of magnetic anisotropy in FGT, suggesting
the likelihood of multiple underlying mechanisms. In 2D magnets, magnetic
anisotropy arising from spin-orbit coupling plays a crucial role in
stabilizing magnetic order. It has been proposed that attaching a nonmagnetic
layer with strong SOC to 2D magnets can effectively enhance the magnetic
anisotropy[27]. Recent experiments involving analogous heterostructures[18,
21] reveal the stabilization of 2D magnetism in thick FGT flakes, which is
induced by the SOC proximity effect. Therefore, it is reasonable to ascribe
the enhancing magnetic anisotropy to the SOC proximity effect. On the other
hand, the magnetic anisotropy of FGT can also be modulated by strain
engineering[31, 32]. Given that the lattice constant in WSe2 is smaller than
that in FGT, a slightly compressive strain may occur in the $ab$-plane of
FGT[24], potentially contributing to the suppression of magnetic anisotropy at
low temperatures.

Figure 4: (a) Comparison of normalized RMCD signal in different
heterostructures with FGT thicknesses of $6.2\mathord{\thinspace\rm nm}$ and
$8.6\mathord{\thinspace\rm nm}$ at $110\mathord{\thinspace\rm K}$. (b)
Relative differences of the coercive field $\Delta
H_{\mathrm{c}}/H_{\mathrm{c}}^{\mathrm{W}}$ as a function of temperature. (c)
Thickness dependence of $T_{\mathrm{c}}$ in the bare FGT and the WSe2-covered
FGT.

Figure 4a shows the RMCD measurements in two WSe2+FGT heterostructures with
varying thicknesses of FGT but the same thickness of WSe2 flakes (2 layers).
The heterostructure with a thinner FGT exhibits a more pronounced enhancement
of $H_{\mathrm{c}}$ at $110\mathord{\thinspace\rm K}$. We present the relative
difference $\Delta H_{\mathrm{c}}/H_{\mathrm{c}}^{\mathrm{W}}$ as a function
of temperature in Fig. 4b. It is noteworthy that the reduction of
$H_{\mathrm{c}}$ at low temperatures is only observable in the thinnest FGT
($6.2\mathord{\thinspace\rm nm}$), but is absent in the heterostructures with
thicker FGT, which implies that the potential strain effect rapidly diminishes
as the FGT thickness increases. Figure 4c illustrates the dependence of
$T_{\mathrm{c}}$ on FGT thicknesses. At the thickness of
$8.6\mathord{\thinspace\rm nm}$, $T_{\mathrm{c}}$ of WSe2/FGT and bare FGT are
nearly identical; however, a significant difference is observed at the
thickness of 6.2 nm, showing that the enhancement induced by SOC proximity is
strengthened with decreasing thickness of FGT.

To explore more details of this interfacial SOC proximity effect, we
fabricated another heterostructure H2 for RMCD measurements, consisting of a
monolayer WSe2 and a $7.8\mathord{\thinspace\rm nm}$-thick FGT layer. A thin
h-BN was used to covered the heterostructure as a comparison. Figures 5a,b
show the normalized RMCD loops at several fixed temperatures (solid curves for
WSe2-covered FGT and dashed curves for bare FGT). As the temperature
increases, H2 generates four distinct regions marked by Roman numerals. In
region I, the enhancement of $H_{\mathrm{c}}$ is present but diminishes above
$20\mathord{\thinspace\rm K}$ in region II. No observable difference of
$H_{\mathrm{c}}$ persists until $85\mathord{\thinspace\rm K}$. In region III,
the enhanced $H_{\mathrm{c}}$ recovers at $90\mathord{\thinspace\rm K}$,
accompanied by the emergence of a two-step magnetization reversal (as marked
by red arrows in Fig. 5b). Unlike the rectangle loop in the bare FGT, the
WSe2-covered FGT exhibits an evident two-step hysteresis loop at about
$95\mathord{\thinspace\rm K}$, with two distinct coercive fields
$H_{\mathrm{c1}}$ and $H_{\mathrm{c2}}$ ($H_{\mathrm{c1}}<H_{\mathrm{c2}}$)
during the magnetization process. The two-step loop is unstable at elevated
temperatures, vanishing at $110\mathord{\thinspace\rm K}$ but reappearing in
region IV. We depict the modulation as a function of temperature in Fig. 5c,
defining the modulation strength as the ratio of $H_{\mathrm{c}}$ in the
WSe2-covered FGT to that in the bare FGT. The two-step loop is only observable
with separate coercive fields at higher temperatures. In addition, the
enhancement of $H_{\mathrm{c}}$ in H2 with hBN encapsulation confirms that the
altered magnetic order results from the interfacial proximity effect rather
than the degeneration (see Figure S3).

The two-step magnetization reversal in the ferromagnetic layer can be achieved
by spin exchange coupling with an adjacent antiferromagnetic layer[33, 34].
However, this is not the scenario in this work as WSe2 is a nonmagnetic
semiconductor. Figure 5d shows the two-step hysteresis loop manifesting four
stable magnetic states, labeled by red numbers from 1 to 4. The additional
intermediate states (state 2 and state 4) could be explained by the short-
range nature of the SOC proximity effect. Assuming that the FGT layers fall
into two parts due to the proximity of WSe2 (see the schematics in Fig. 5e).
For the FGT layers at the interface, the SOC proximity effect is dominant,
leading to enhanced magnetic anisotropy and a larger coercive field. While the
underneath FGT layers are minimally affected by the interface, resulting in a
smaller magnetic field required for flipping. Thus the arrangement of magnetic
moments can be illustrated by schematics in Fig. 5e. States 1 and 3 exhibit
fully spin-polarized magnetization with all-spin-up and all-spin-down
configurations, respectively. While the magnetizations in the layers adjacent
to and far away from the interface point in opposite directions, resulting in
intermediate states 2 and 4. A similar two-step hysteresis loop induced by SOC
proximity was recently observed in the MoS2/FGT heterostructure[19].

Figure 5: Temperature dependence of the SOC proximity effect in H2. (a,b)
Comparison of normalized RMCD signal for bare FGT (dashed curves) and
WSe2-covered FGT (solid curves) at several temperatures. (c) The ratio of
$H_{\mathrm{c}}$ in WSe2-covered FGT to that in bare FGT as a function of
temperature. Separated parts of two curves correspond to the emergence of the
magnetic intermediate states. (d) RMCD signal for WSe2-covered FGT at
$135\mathord{\thinspace\rm K}$. (e) Schematics of moment configurations for
the corresponding magnetic states in (d). The layer with solid green circles
represents the monolayer WSe2, and the solid circles with arrows represent the
magnetization orientations of different FGT layers. The gradient blue
background depicts the influence induced by SOC proximity.

At low temperatures, characterized by reduced thermal fluctuation, the
interlayer exchange coupling is sufficiently strong to maintain a uniform
$H_{\mathrm{c}}$ in FGT layers. However, the interlayer exchange coupling
diminishes with strong thermal fluctuation at elevated temperatures, resulting
in non-uniform coercive fields for the interfacial and underneath FGT layers,
i.e., the two-step magnetization reversal process. At much lower temperatures
($<20\mathord{\thinspace\rm K}$), the proximity effect exhibits a non-
monotonic modification depending on the thickness of the WSe2 and the FGT
layers. The enhancement of $H_{\mathrm{c}}$ emerges in monolayer-WSe2/FGT, but
disappears in others. While the reduced $H_{\mathrm{c}}$ is observed with the
thinnest FGT ($6.2\mathord{\thinspace\rm nm}$). This may suggest a complex
competition between different mechanisms including interlayer exchange
coupling, SOC proximity, interfacial strain, thermal fluctuation and so on,
which requires further investigation in the future.

In conclusion, we systematically investigated the proximity effect on 2D
magnetism in WSe2/FGT vdW heterostructures. Both $H_{\mathrm{c}}$ and
$T_{\mathrm{c}}$ can be noticeably enhanced in the WSe2-covered region, which
is more pronounced as the FGT thickness decreases. This enhancement is
attributed to the SOC proximity effect induced from the adjacent WSe2. At low
temperatures, we observed a reduction of the coercive field in the
heterostructure with a $6.2\mathord{\thinspace\rm nm}$ thick FGT flake,
potentially due to the strain effect from the interface. In addition, the
presence of unconventional magnetic stratification demonstrates the short-
range nature of the SOC proximity effect. Our findings offer a practical
method for modifying 2D magnetic order by combining a 2D nonmagnetic
semiconductor with strong SOC, thereby promising for advancing the field of
spintronics based on vdW magnetic heterostructures.

Financial support from the National Key R&D Program of China (Nos.
2022YFA120470, 2021YFA1400400), the National Natural Science Foundation of
China (Nos. 12004173, 11974169), and the Fundamental Research Funds for the
Central Universities (Nos. 020414380087, 020414913201) are gratefully
acknowledged.

## References

  * [1] Prinz G A 1990 Science 250 1092–1097
  * [2] $\mathrm{\check{Z}}$uti$\mathrm{\acute{c}}$ I, Fabian J and Das Sarma S 2004 Rev. Mod. Phys. 76 323–410
  * [3] Datta S and Das B 1990 Appl. Phys. Lett. 56 665–667
  * [4] Moodera J S, Kinder L R, Wong T M and Meservey R 1995 Phys. Rev. Lett. 74 3273–3276
  * [5] Miyazaki T and Tezuka N 1995 J. Magn. Magn. Mater. 139 L231–L234
  * [6] Ahn E C 2020 npj 2D Mater. Appl. 4 17
  * [7] Sierra J F, Fabian J, Kawakami R K, Roche S and Valenzuela S O 2021 Nat. Nanotech. 16 856–868
  * [8] Zhang L, Zhou J, Li H, Shen L and Feng Y P 2021 Appl. Phys. Rev. 8 021308
  * [9] Ghising P, Biswas C and Lee Y H 2023 Adv. Mater. 35 2209137
  * [10] Novoselov K S, Mishchenko A, Carvalho A and Castro Neto A H 2016 Science 353 aac9439
  * [11] Liu Y, Weiss N O, Duan X, Cheng H C, Huang Y and Duan X 2016 Nat. Rev. Mater. 1 16042
  * [12] Liu Y, Huang Y and Duan X 2019 Nature 567 323–333
  * [13] Zhong D, Seyler K L, Linpeng X, Cheng R, Sivadas N, Huang B, Schmidgall E, Taniguchi T, Watanabe K, McGuire M A, Yao W, Xiao D, Fu K M C and Xu X 2017 Sci. Adv. 3 e1603113
  * [14] Zhong D, Seyler K L, Linpeng X, Wilson N P, Taniguchi T, Watanabe K, McGuire M A, Fu K M C, Xiao D, Yao W and Xu X 2020 Nat. Nanotech. 15 187–191
  * [15] Zollner K, Faria Junior P E and Fabian J 2019 Phys. Rev. B 100 085128
  * [16] Li L, Jiang S, Wang Z, Watanabe K, Taniguchi T, Shan J and Mak K F 2020 Phys. Rev. Mater. 4 104005
  * [17] Lyons T P, Gillard D, Molina-Sánchez A, Misra A, Withers F, Keatley P S, Kozikov A, Taniguchi T, Watanabe K, Novoselov K S, Fernández-Rossier J and Tartakovskii A I 2020 Nat. Commun. 11 6021
  * [18] Kim S J, Choi D, Kim K W, Lee K Y, Kim D H, Hong S, Suh J, Lee C, Kim S K, Park T E and Koo H C 2021 ACS Nano 15 16395–16403
  * [19] Tu Z, Zhou T, Ersevim T, Arachchige H S, Hanbicki A T, Friedman A L, Mandrus D, Ouyang M, Žutić I and Gong C 2022 Appl. Phys. Lett. 120 043102
  * [20] Gong C, Zhang P, Norden T, Li Q, Guo Z, Chaturvedi A, Najafi A, Lan S, Liu X, Wang Y, Gong S J, Zeng H, Zhang H, Petrou A and Zhang X 2023 Nat. Commun. 14 3839
  * [21] Choi E M, Kim T, Cho B W and Lee Y H 2023 ACS Nano 17 15656–15665
  * [22] Huang B, Clark G, Navarro-Moratalla E, Klein D R, Cheng R, Seyler K L, Zhong D, Schmidgall E, McGuire M A, Cobden D H, Yao W, Xiao D, Jarillo-Herrero P and Xu X 2017 Nature 546 270–273
  * [23] Fei Z, Huang B, Malinowski P, Wang W, Song T, Sanchez J, Yao W, Xiao D, Zhu X, May A F, Wu W, Cobden D H, Chu J H and Xu X 2018 Nat. Mater. 17 778–782
  * [24] Wang H, Liu Y, Wu P, Hou W, Jiang Y, Li X, Pandey C, Chen D, Yang Q, Wang H, Wei D, Lei N, Kang W, Wen L, Nie T, Zhao W and Wang K L 2020 ACS Nano 14 10045–10053
  * [25] Katmis F, Lauter V, Nogueira F S, Assaf B A, Jamer M E, Wei P, Satpati B, Freeland J W, Eremin I, Heiman D, Jarillo-Herrero P and Moodera J S 2016 Nature 533 513–516
  * [26] Zhu W, Song C, Han L, Bai H, Wang Q, Yin S, Huang L, Chen T and Pan F 2022 Adv. Funct. Mater. 32 2108953
  * [27] Dong X J, You J Y, Zhang Z, Gu B and Su G 2020 Phys. Rev. B 102 144443
  * [28] Le D, Barinov A, Preciado E, Isarraraz M, Tanabe I, Komesu T, Troha C, Bartels L, Rahman T S and Dowben P A 2015 J. Phys.: Condens. Matter 27 182201
  * [29] Castellanos-Gomez A, Buscema M, Molenaar R, Singh V, Janssen L, van der Zant H S J and Steele G A 2014 2D Mater. 1 011002
  * [30] Zhao W, Ghorannevis Z, Chu L, Toh M, Kloc C, Tan P H and Eda G 2013 ACS Nano 7 791–797
  * [31] Zhuang H L, Kent P R C and Hennig R G 2016 Phys. Rev. B 93 134407
  * [32] Wang Y, Wang C, Liang S J, Ma Z, Xu K, Liu X, Zhang L, Admasu A S, Cheong S W, Wang L, Chen M, Liu Z, Cheng B, Ji W and Miao F 2020 Adv. Mater. 32 2004533
  * [33] Chen P, Huang Z, Li C, Zhang B, Bao N, Yang P, Yu X, Zeng S, Tang C, Wu X, Chen J, Ding J, Pennycook S J, Ariando A, Venkatesan T V and Chow G M 2018 Adv. Funct. Mater. 28 1801766
  * [34] Wu X, Lan D, Hwang I, Sun C, Zhou H, Yu X, Yang P, Yu X, Liu C, Chen P, Ding J, Chen J and Chow G M 2023 J. Alloy. Compd. 932 167582