Further author information: (Send correspondence to Y.T.)
Y.T.: E-mail<EMAIL_ADDRESS>Telephone: +81 (0)52 789 2846

# Wavefront sensor for millimeter/submillimeter-wave adaptive optics based on
aperture-plane interferometry

Yoichi Tamura  Nagoya University, Nagoya, 464-8602 Japan  Ryohei Kawabe
National Astronomical Observatory of Japan, Mitaka, Tokyo, 181-8588 Japan
Yuhei Fukasaku  University of Tsukuba, Tsukuba, Ibaraki, 305-8573 Japan
Kimihiro Kimura  Nagoya University, Nagoya, 464-8602 Japan  Tetsutaro Ueda
Nagoya University, Nagoya, 464-8602 Japan  Akio Taniguchi  Nagoya University,
Nagoya, 464-8602 Japan  Nozomi Okada  Dept. Physical Science, Osaka Prefecture
University, Sakai, Osaka, 599-8531 Japan  Hideo Ogawa  Dept. Physical Science,
Osaka Prefecture University, Sakai, Osaka, 599-8531 Japan  Ikumi Hashimoto
Dept. Aerospace Engineering, Osaka Prefecture University, Sakai, Osaka,
599-8531 Japan  Tetsuhiro Minamidani  National Astronomical Observatory of
Japan, Mitaka, Tokyo, 181-8588 Japan  Nobeyama Radio Observatory, Nagano,
384-1305 Japan  Noriyuki Kawaguchi  National Astronomical Observatory of
Japan, Mitaka, Tokyo, 181-8588 Japan  Nario Kuno  University of Tsukuba,
Tsukuba, Ibaraki, 305-8573 Japan  Yohei Togami  Nagoya University, Nagoya,
464-8602 Japan  Masato Hagimoto  Nagoya University, Nagoya, 464-8602 Japan
Satoya Nakano  Nagoya University, Nagoya, 464-8602 Japan  Keiichi Matsuda
Nagoya University, Nagoya, 464-8602 Japan  Sachiko Okumura  Japan Women’s
University, Bunkyo, Tokyo, 112-8681 Japan  Tomoko Nakamura  Japan Women’s
University, Bunkyo, Tokyo, 112-8681 Japan  Mikio Kurita  Kyoto University,
Kyoto, 606-8501 Japan  Tatsuya Takekoshi  Kitami Institute of Technology,
Kitami, Hokkaido, 090-8507 Japan  Institute of Astronomy, The University of
Tokyo, Mitaka, Tokyo, 181-0015 Japan  Tai Oshima  National Astronomical
Observatory of Japan, Mitaka, Tokyo, 181-8588 Japan  Toshikazu Onishi  Dept.
Physical Science, Osaka Prefecture University, Sakai, Osaka, 599-8531 Japan
Kotaro Kohno  Institute of Astronomy, The University of Tokyo, Mitaka, Tokyo,
181-0015 Japan

###### Abstract

We present a concept of a millimeter wavefront sensor that allows real-time
sensing of the surface of a ground-based millimeter/submillimeter telescope.
It is becoming important for ground-based millimeter/submillimeter astronomy
to make telescopes larger with keeping their surface accurate. To establish
‘millimetric adaptive optics (MAO)’ that instantaneously corrects the
wavefront degradation induced by deformation of telescope optics, our
wavefront sensor based on radio interferometry measures changes in excess path
lengths from characteristic positions on the primary mirror surface to the
focal plane. This plays a fundamental role in planed 50-m class submillimeter
telescopes such as LST and AtLAST.

###### keywords:

Submillimeter, single-dish telescope, adaptive optics, aperture-plane
interferometry

##  1 INTRODUCTION

The increase in size of telescopes, i.e., the acquisition of larger collecting
area and spatial resolution, is the very history of observational astronomy.
For millimeter and submillimeter-wave telescopes (mm/submm telescopes,
hereafter), increasing the antenna diameter with keeping the mirror surface
accurate is important in developing new astronomical fields. The advent of an
aperture synthesis interferometer has resulted in significant improvement in
collecting area and angular resolution. For example, the Atacama Large
Millimeter/submillimeter Array (ALMA) realized the performance by keeping the
size of the array element antennas moderate ($D=7$ and 12 m) and increasing
stiffness of them to keep their surface down to $20$ $\mu$m r.m.s or even
better.

However, aperture synthesis interferometry is not a panacea in realizing a
larger-scale mm/submm telescope. The focal plane instruments of radio
interferometers are limited to coherent receivers (e.g., heterodyne receivers)
that can detect the phase of celestial signals. On the other hand, large-
format arrays of direct photon detectors (e.g., cameras and integrated
superconducting spectrometers[1]), which have prospered in recent years,
cannot be accommodated as receivers for an aperture synthesis interferometer,
because phase detection is not possible. Thus, it becomes increasingly
important to have a large-aperture single-dish mm/submm telescope (e.g., Large
Submillimeter Telescope[2] and Atacama Large Aperture Submillimeter
Telescope[3]), especially in the field in which direct photon detector arrays
play a crucial role.

A limiting factor in size and operating frequency of mm/submm telescopes is
the deterioration of optical performance due to changes in the environment
surrounding the telescope, such as wind load and thermal deformation in
addition to gravitational deformation.111In the millimeter and submillimeter
regime, the inhomogeneities in column density of water vapor in the
troposphere also cause wavefront degradation. However, the spatial scale at
which the water vapor fluctuations dominate the wavefront degradation is
typically much longer than the telescope’s aperture, for example, $\gtrsim
10^{2}$ m at the ALMA site. For this reason, in most cases, the Earth’s
atmosphere is not a major source of wavefront degradation for single-dish
mm/submm telescopes. This is unlike the visible/near-infrared (NIR) regime,
where the major cause of wavefront degradation is due to temperature
inhomogeneities of the tropospheric atmosphere on the order of tens of cm.
Mm/submm telescopes have been manufactured on the premise of large ($>10$ m)
structures in order to ensure the large collecting area and high spatial
resolution. For this reason, building a huge dome to cover it would require
enormous costs. Therefore, it was an issue to construct an ‘exposed’ antenna
outside and ensure its mirror surface with an accuracy of $\sim 20$ $\mu$m.

With the existing technology, it is possible to measure the mirror surface
shape in advance by the photogrammetry or radio holography methods, and to
correct the mirror surface down to the level of $\sim 10$ $\mu$m (r.m.s.) by
adjusting the positions of the primary mirror panels. In recent years, the
weight reduction and enhanced stiffness of the primary mirror structure and
the active surface control of the primary mirror panels by the motorized
actuators have been used to keep the primary mirror an ideal shape with
respect to changes in telescope elevation angle.[4]

The problem is, however, the deformation of the primary mirror surface and the
secondary mirror support structure due to wind load and thermal deformation.
This is because the time-scale of the deformation is fairly short compared
with a typical time-scale of astronomical observations. It is difficult to
measure the mirror surface shape in real time during the observations. It is
likely that most of the deformation across the primary mirror are dominated by
lower-order deformation modes. Their typical spatial scale is found to be a
fraction of the aperture diameter (1–10 m)[5]. The time-scale of the
deformation is determined by the natural frequency of the primary mirror
structure to the wind load and is typically $\sim 10^{-1}$–1 s (e.g., $\simeq
1.1$ s for the Nobeyama 45 m mm-wave telescope[6]). In this case, the control
technology of the adaptive optics does not become a major obstacle. The
spatial and time-scales of adaptive optics control, which is already
established in NIR astronomy, are on the order of $\sim 10$ cm and $\sim 1$
kHz, respectively, which are several orders of magnitude higher in spatial and
temporal frequencies than those of large mm/submm telescopes. Therefore, in
the radio telescope, the wavefront compensation can be performed at the
primary or secondary mirror equipped with the active surface control, instead
of a high-speed deformable mirror employed in NIR telescopes.

Then, the major obstacle in the mm/submm regime lies in _how to measure the
wavefront in real time_. Wavefront sensing using a Shack-Hartmann sensor is
the mainstream in the optical and NIR. A similar approach, however, is not
applicable to MAO, since no large-format detector array that can be
manufactured cheaply is available in the mm/submm. Instead, radio astronomy
can exploit radio interferometry as a native wavefront sensing technology,
which measures the difference between the arrival times of the wavefronts
coming through two independent optical paths.

In this paper, we present a concept of a mm wavefront sensor that allows real-
time sensing of the primary mirror surface of a ground-based large-aperture
mm/submm telescope.

Figure 1: A schematic diagram of aperture-plane interferometry. (a) A simple
interferometer with a correlator which measures a difference in arrival time
between two signals coming from a common reference source. (b) The same as (a)
but one of the optical paths coming through the free space as shown with the
dashed line. (c) The same as (b) but the path through the free space goes
through the telescope optics. The lengths of the other path is adjusted by a
delay line shown as a loop. (d) The same as (c) with a switch followed by an
array of radiators which enables to measure the distances from arbitrary
positions to the focal point.

##  2 APERTURE-PLANE INTERFEROMETRY FOR WAVEFRONT SENSING

As seen in the previous section, what needs to be measured first is the
relative time variation of the deformation of the optical system, rather than
wavefronts of celestial signals coming through the atmosphere. This is
equivalent to the errors in excess path length (EPL) from the focus to
arbitrary points on the primary surface. Therefore, the wavefront sensor
employs multiple reference microwave/mm sources placed on the aperture of the
telescope as a phase standard. The goal of the surface accuracy we measure
with a prototype sensor (see later sections) is 40 $\mu$m r.m.s. This is close
to $\approx 100~{}\mu$m r.m.s., which is routinely achieved by holographic
measurements of the Nobeyama 45 m mm-wave telescope. This measurement accuracy
is equivalent to the phase accuracy of 1∘ r.m.s. when operating the wavefront
sensor at 20 GHz.

Figure 1 shows the schematic diagram representing the principle of aperture-
plane interferometry. A simple interferometer comprising a correlator (figure
1a, b) allows us to measure a difference in arrival time between two signals
generated by a common reference source. The interferometer in which one of the
paths goes through the telescope optics (figure 1c) measures the EPL from a
certain point on the primary surface to the focus of the telescope. This is
similar to what is used for phase calibration (the p-cal method, REF) of
phase-referencing receivers used in VLBI Exploration of Radio Astronomy
(VERA[7]). The method is a scalable technique and can be time-multiplexed by
an intervening switch which is followed by a series of radiators placed across
the surface of the primary mirror (figure 1d).

The required number of reference points (i.e., radiators) depends on to what
order of deformation modes needs to be characterized. As the wavefront sensor
detects the relative deformation with respect to the ideal surface realized
when the telescope structure is static, the surface needs to be adjusted by a
conventional method such as the radio holography in advance. The relative
deformation of the structure induced by wind and thermal loads emerges mostly
on large spatial scales (e.g., a fraction of an aperture size).[5] Therefore,
a couple of tens references across the aperture will be good enough.

We opt broadband noise as the common reference source, rather than a
continuous wave (CW) or narrowband signals. This is because in the case of CW
or narrowband signals, the stray lights through multiple paths can interfere
with each other, resulting in considerable systematic errors in measurements
of EPLs. On the other hand, broadband noises are more robust against the
multi-path interference, and thus they are often used for delay calibration in
aperture-synthesis interferometers such as ALMA. The EPL is imprinted as a
phase slope in a cross power spectrum (CPS) of the broadband noise. Let
$C(\nu;\,\mathrm{EPL})=A(\nu)\exp{[i\phi(\nu;\,\mathrm{EPL})]}$ be a measured
CPS. Then, a calibrated CPS is obtained as
$C(\nu;\,\mathrm{EPL})/C(\nu;0)=\exp{[i\\{\phi(\nu;\,\mathrm{EPL})-\phi(\nu;\,0)\\}]}$.
The phase of a calibrated CPS,
$\Delta\phi(\nu;\,\mathrm{EPL})\equiv\phi(\nu;\,\mathrm{EPL})-\phi(\nu;\,0)$,
induced by a certain EPL is expressed as a function of frequency $\nu$ as

$\displaystyle\Delta\phi(\nu)=2\pi\tau\nu=\frac{2\pi\,\mathrm{EPL}}{c}\nu,$
(1)

where $\tau$ is the time delay induced by EPL. This means that the slope of
the CPS phase immediately gives the EPL.

Figure 2: (a) The system block diagram of the prototype wavefront sensor for
millimetric adaptive optics. (b) The configuration of the wavefront sensor
system mounted on the Nobeyama 45 m telescope.

##  3 SYSTEM

Here we describe the system of the prototype wavefront sensor for
demonstration. Although this prototype only has several elements of radiators
operating at 20 GHz, which is relatively low frequency and easy to handle, the
aperture-plane interferometry is scalable up to tens to hundreds of elements
at higher frequencies.

###  3.1 Requirement and specifications

The top-level requirement on the system is to instantaneously measure the
deviation from the ideal mirror surface with an accuracy of 40 $\mu$m r.m.s.
with a time-resolution of 100 ms, well below the natural frequency of the
primary mirror structure ($\approx 1$ Hz). This requirement corresponds to the
phase accuracy of 1 deg r.m.s. for an operating frequency of 20 GHz. The error
budget is split into (1) the statistical error of thermal noises arising from
the reference source and the receiver noise (appendix A) and (2) the
systematic error due to short term ($\sim 10$ s) stability of EPL
measurements, whereas long term stability needs to be taken into account in
the actual implementation of the system. Table 1 summarizes the specifications
of the prototype wavefront sensor.

Table 1: Specifications of the prototype wavefront sensor. Item | Value
---|---
Frequency (GHz) | 16–24
Number of elements | 5

###  3.2 Subsystems

Figure 2 shows the schematic block diagram and configuration describing how
the wavefront sensor works. The system is threefold; (1) the transmitter
subsystem, which generates and transfers a reference signal, and injects it
into the telescope optics from multiple positions across the primary surface;
(2) the receiver subsystem, which collects the signal coming through the
telescope optics and sends it to the correlator; and (3) the correlator
subsystem, which differentiates the arrival time of the reference signals
coming directly from the reference generator and through the telescope optics.

Figure 3: (a) The setup of laboratory evaluation. (b) A feed horn mounted on a
micrometer stage. (c) The OCTAD-M correlator.

####  3.2.1 Transmitter subsystem

The transmitter subsystem comprises a reference noise generator, optical
transmitter and switch, and radiators on the primary mirror.

Reference noise generator. For the reference broadband noise, we use Johnson-
Nyquist noise of a microwave terminator amplified by a cascade of power
amplifiers which is followed by a bandpass filter (17.3–23.6 GHz). The excess
noise ratio (ENR) is $\approx 70$ dB. The reference noise is divided, one of
which goes directly into the correlator while the other is sent to the
telescope primary mirror through the optical fiber.

Optical fiber transmission and distribution. The reference signal is
transmitted on optical fiber, since the distances to the radiators are too
long to maintain the phase stability and gain of the reference signal. We use
a radio-frequency on fiber (RFoF) system (Optilab, RFLL-20-H), a pair of
optical modulator (E/O) and demodulator (O/E). We use a phase-stabilized
single-mode fiber (Sumitomo Electric) with a linear expansion coefficient of
$4.7\times 10^{-6}$ K-1. The optical signal leaving the E/O is coupled with
one of the five O/Es placed on the backup structure of the primary mirror by
switching over the fibers with a MEMS-based optical switch placed closed to
the central hub of the primary mirror.

Feed horn on the primary mirror. The reference signal is demodulated with the
O/E and is radiated by one of five identical feed horns placed on the surface
of the primary mirror. We choose a linear-polarized feed horn as a reference
signal radiator as the operating frequency is broad ($\Delta f/f=0.4$). The
gain and half-power beam-width at 20 GHz are 20 dBi and $10^{\circ}$,
respectively. The feed horns are placed at a radius of 16 m with position
angles of $-18^{\circ}$ (top), $90^{\circ}$ (left), $180^{\circ}$ (bottom),
and $270^{\circ}$ (right), and at $r=5$ m with a position angle of $0^{\circ}$
(center; see also figure 2b).

####  3.2.2 Receiver subsystem

We use the H22
receiver,222https://www.nro.nao.ac.jp/~nro45mrt/html/prop/status/Status_latest.html
a science-grade cryogenic coherent receiver operating at 20 GHz. The typical
system noise temperature is $T_{\mathrm{sys}}=100$ K. H22 detects two circular
polarization, while we only use one of them for phase measurements. The
receiver output is amplified and the bandwidth is constrained by the anti-
aliasing filter before the signal goes into the correlator.

####  3.2.3 Correlator subsystem

We develop a FX-type digital correlator, OCTAD-M (ELECS Inc.), which is based
on the architecture initially developed for the FPGA-based fast Fourier
transform (FFT) spectrometer OCTAD-S[8]. OCTAD-M is equipped with two 3-bit
analog-to-digital converters with 16.384 GSa/s. We opt the third-order mode
(16.384–24.576 GHz) of analog input signals, allowing us to directly sample
the H22 receiver output with no down conversion. Two input signals are FFTed
first and are multiplied to obtain the cross-power spectrum. The spectra are
accumulated for 5 or 10 msec. OCTAD-M has a digital delay capability with up
to $2^{15}$ times sampling clock time, allowing to digitally insert an
instrumental delay of 0–2 $\mu$sec, corresponding to 0–599.585 m for free-
space geometrical delay.

####  3.2.4 Computing/control subsystem

To synchronize the optical switch and the correlator sampling, we use common 1
PPS and 10 MHz reference signals fed by the observatory GPS server. System
control and data acquisition are performed by a single Linux server.

###  3.3 Ancillary devices

To help the wavefront sensing experiment, we use ancillary devices placed on
and around the Nobeyama 45 m telescope. We have six piezoelectric
accelerometers attached to the backup structure of the primary mirror[9],
which are placed close to the edge of the primary mirror with position angles
of $0^{\circ}$ (top), $45^{\circ}$, $90^{\circ}$ (left), $180^{\circ}$
(bottom), $270^{\circ}$ (right), and close to the central hub (see figure 2b).
We also use a weather monitor placed at the top of a 50-m tall meteorological
tower, which is located 75 m north of the
telescope.333https://www.nro.nao.ac.jp/~nro45mrt/html/obs/weather/

Figure 4: (a) The cross power spectrum (CPS) phases measured for different
excess path lengths (EPLs) induced by mechanical deviation of the micrometer
stage $z$. (b) The EPL measured by the phase slope as a function of $z$. (c)
The residual, EPL$-z$, as a function of $z$. The error bar on each plot
accounts for the $1\sigma$ error in linear regression of the CPS phase. The
average of the residuals is 4.8 $\mu$m with a standard deviation of 17.4
$\mu$m. The grey region represents the system requirement ($\pm 40$ $\mu$m).

##  4 DEMONSTRATION

###  4.1 Laboratory evaluation of linearity and accuracy in excess path length
measurements

In advance to the on-site demonstration, we evaluate the linearity, accuracy,
and stability of the system with respect to mechanical change in EPL. Although
neither the actual telescope optics nor the H22 receiver is available in the
lab, we use almost the same configuration of the other subsystems which are
used for the on-site measurements. The configuration of the lab experiment is
shown in figure 3. Here we use a feed horn followed by a power amplifier and a
bandpass filter instead of the H22 receiver. We placed the two feed horns of
the radiator and receiver $\approx 2$ m away from each other. The radiator
feed horn is placed on the micrometer stage that slides on a straight rail to
mechanically change the EPL between the radiator and receiver, mimicking the
deformation of the telescope optics. The mechanical accuracy of the micrometer
stage is 3 $\mu$m.

We measure the CPS phase at $z=0$, 1, 2, …, 15 mm, where $z$ is the EPL
induced by the mechanical deviation of the micrometer stage. The CPS is
calibrated for the complex bandpass (i.e., amplitude and phase) by dividing
the CPS by that of $z=0$. Figure 5a shows the calibrated CPS phases obtained
at $z=1$, 2, 3, …, 15 mm. The CPS is running-averaged with a rectangle window
of 300 MHz. We see a linear slope for each CPS phase although a small
fluctuation remains around the best-fitting line. Since the statistical error
of the CPS phase is almost negligible, the fluctuation could be dominated by
systematic errors such as standing waves, multi-path interference, or phase
drift of the system. As shown in figure 5b, the EPL derived from equation 2 is
linear with respect to the mechanical deviation $z$. Figure 5c shows the
residual of the measured EPL subtracted by $z$ as a function of $z$. The
average of the residuals is 4.8 $\mu$m with a standard deviation of 17.4
$\mu$m, which is well below the system requirement of $<40$ $\mu$m, confirming
the system linearity and accuracy.

Figure 5: (a) The cross power spectrum (CPS) phase obtained under a very windy
(5–10 m s-1) condition and (b) under a moderately windy ($\sim 5$ m s-1)
condition. The time origins are 3:21:00 on 2020 November 22 UTC and 2:09:00 on
2020 November 23 UTC, respectively. (Top) The CPS phase measured at the ‘top’
position of the 45 m primary mirror (blue curve). As the CPS is measured at 20
GHz ($\lambda\approx 15$ mm), a phase shift of 10 deg approximately
corresponds to a surface deviation of 400 $\mu$m. The grey curve represents
its second-order differential in arbitrary units. (Center) The acceleration
measured with two accelerometers attached to the top edge (blue) and the
center (purple) of the 45 m primary mirror. (Bottom) The wind speed measured
with a weather monitor placed atop the 50-m tall tower located 75 m north of
the 45 m telescope.

###  4.2 Demonstration at the Nobeyama 45 m telescope

We carried out a commissioning campaign of the MAO wavefront sensor with the
Nobeyama 45 m telescope in 2020 November. The goals are to confirm that the
wavefront sensor actually works as designed and to demonstrate that aperture
plane interferometry allows us to measure an EPL in real time. In this
campaign, only 2 radiators were installed at the ‘top’ and ‘center’ positions.

The top pannels of figures 5a and 5b show the temporal changes in CPS phase at
20 GHz, which were taken under very windy (5–10 m s-1) and moderate ($\sim 5$
m s-1) conditions, respectively. A large phase drift and small $\sim 3$ Hz
ripples in it are clearly seen in the very windy condition, while the drift is
rather suppressed under the moderate condition. This fact suggests that the
CPS phase measures the EPL change induced by the wind load, although further
analyses are necessary. This is also supported by the measurements of the
piezoelectric accelerometer at the ‘top’ position; a constant $\sim 3$ Hz
oscillation in addition to intermittent slow ($\sim 1.5$ Hz), high-amplitude
oscillations due to a series of gust are evident in figure 5a, whereas no
clear oscillation is seen in figure 5b.

##  5 CONCLUSIONS

We present a concept of a millimeter wavefront sensor which allows real-time
sensing of the primary mirror surface of a future 50-m class mm/submm
telescope, such as LST and AtLAST. To establish millimetric adaptive optics
(MAO) that instantaneously corrects the wavefront degradation induced by
deformation of telescope optics due to wind and thermal loads, our wavefront
sensor employs aperture-plane interferometry to measure real-time changes in
EPLs from characteristic points on the primary surface to the focal plain. The
proposed wavefront sensor operates at 16–24 GHz, which is cost-effective and
is accurate enough to measure the EPL down to the level of 40 $\mu$m r.m.s. as
demonstrated in the laboratory evaluation. Although the verification is in
progress, we have demonstrated that a 2-element pilot wavefront sensor with
the Nobeyama 45 m telescope worked properly and detected the EPL change
induced by wind load on the telescope structure.

##  Appendix A SENSITIVITY

In general, the variance of the cross power spectrum is given by

$\sigma^{2}=\frac{(S_{1}+N_{1})\cdot(S_{2}+N_{2})+S_{1}S_{2}}{2\Delta\nu\,\Delta
t},$

where $S_{i}$ and $N_{i}$ ($i=1,\,2$) are the correlated and non-correlated
components of the cross power spectrum of two incident signals, respectively.
Therefore, the S/N of the cross power spectrum is

$\mathrm{SNR}=\frac{\sqrt{S_{1}S_{2}}\sqrt{2\Delta\nu\,\Delta
t}}{\sqrt{(S_{1}+N_{1})\cdot(S_{2}+N_{2})+S_{1}S_{2}}}.$

Here, the correlated and non-correlated components input directly to the
correlator from the reference signal source (reference continuum source) are
$S_{1},\,N_{1}$, and components input to the correlator via the transmission
system and optical system are $S_{2},\,N_{2}$, then $S_{1}\gg N_{1}$, so the
S/N per spectral channel of the cross power spectrum can be approximated as

$\mathrm{SNR}\approx\frac{\sqrt{2\Delta\nu\,\Delta
t}}{\sqrt{N_{2}/S_{2}+2}}\sim\sqrt{2\Delta\nu\,\Delta t\,(S_{2}/N_{2})}\approx
8.16\left(\frac{\Delta\nu}{\rm 1~{}MHz}\right)^{0.5}\left(\frac{\Delta t}{\rm
10~{}ms}\right)^{0.5}\left(\frac{N_{2}/S_{2}}{300}\right)^{-0.5},$

where $N_{2}/S_{2}\sim O(10^{2})\gg 2$ is assumed for the second
approximation. If we assume that S/N of signal estimated using all
$n_{\mathrm{s}}$ spectral channels is scaled with $\sqrt{n_{\mathrm{s}}}$ for
S/N per channel, then

$\mathrm{SNR_{tot}}\approx 739\left(\frac{\Delta\nu}{\rm
8192~{}MHz}\right)^{0.5}\left(\frac{\Delta t}{\rm
10~{}ms}\right)^{0.5}\left(\frac{N_{2}/S_{2}}{300}\right)^{-0.5}.$

Assuming that thermal noise dominates $N_{2}$ and the cross power spectrum has
isotropic dispersion on the complex plane, the uncertainty of the phase
$\theta$ is $\sigma_{\theta}\sim 1/\mathrm{SNR_{tot}}$ (radian). Therefore,

$\sigma_{\theta}\approx 0.078\ \mathrm{(deg)}\times\left(\frac{\Delta\nu}{\rm
8192~{}MHz}\right)^{-0.5}\left(\frac{\Delta t}{\rm
10~{}ms}\right)^{-0.5}\left(\frac{N_{2}/S_{2}}{300}\right)^{0.5}$

If the value is fiducial, sufficient phase measurement accuracy can be
achieved. In other words, it is necessary to obtain a correlator input signal
that satisfies $N_{2}/S_{2}<300$.

###### Acknowledgements.

We acknowledge K. Handa, C. Miyazawa, T. Kanzawa, T. Wada, and M. Saito for
their support. This work was supported by JSPS KAKENHI (Grant No. 17H06206)
and NAOJ Research Coordination Committee, NINS.

## References

  * [1] Endo, A., Karatsu, K., Laguna, A. P., Mirzaei, B., Huiting, R., Thoen, D. J., Murugesan, V., Yates, S. J. C., Bueno, J., Marrewijk, N. v., Bosma, S., Yurduseven, O., Llombart, N., Suzuki, J., Naruse, M., de Visser, P. J., van der Werf, P. P., Klapwijk, T. M., and Baselmans, J. J. A., “Wideband on-chip terahertz spectrometer based on a superconducting filterbank,” Journal of Astronomical Telescopes, Instruments, and Systems 5, 035004 (July 2019).
  * [2] Kawabe, R., Kohno, K., Tamura, Y., Takekoshi, T., Oshima, T., and Ishii, S., “New 50-m-class single-dish telescope: Large Submillimeter Telescope (LST),” in [Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 9906, 990626 (2016).
  * [3] Klaassen, P. D., Mroczkowski, T. K., Cicone, C., Hatziminaoglou, E., Sartori, S., De Breuck, C., Bryan, S., Dicker, S. R., Duran, C., Groppi, C., Kaercher, H., Kawabe, R., Kohno, K., and Geach, J., “The Atacama Large Aperture Submillimeter Telescope (AtLAST),” in [Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 11445, 114452F (Dec. 2020).
  * [4] Hughes, D. H., Jáuregui Correa, J.-C., Schloerb, F. P., Erickson, N., Romero, J. G., Heyer, M., Reynoso, D. H., Narayanan, G., Perez-Grovas, A. S., Souccar, K., Wilson, G., and Yun, M., “The Large Millimeter Telescope,” in [Ground-based and Airborne Telescopes III ], Stepp, L. M., Gilmozzi, R., and Hall, H. J., eds., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 7733, 773312 (July 2010).
  * [5] Levy, R., [Structural engineering of microwave antennas for electrical, mechanical, and civil engineers ] (1996).
  * [6] Smith, D. R., Paglione, T. A., Lovell, A. J., Ukita, N., and Matsuo, H., “Measurements of dynamic pointing variations of a large radio telescope,” in [Radio Telescopes ], Butcher, H. R., ed., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 4015, 467–475 (July 2000).
  * [7] Kobayashi, H., Sasao, T., Kawaguchi, N., Manabe, S., Omodaka, T., Kameya, O., Shibata, K. M., Miyaji, T., Honma, M., Tamura, Y., Hirota, T., Kuji, S., Horiai, K., Sakai, S., Sato, K., Iwadate, K., Kanya, Y., Ujihara, H., Jike, T., Fujii, T., Motiduki, N., Oyama, T., Kurayama, H., Kamohara, R., Suda, H., and Kasuga, T., [VERA: A New VLBI Instrument Free from the Atmosphere ], vol. 306 of Astronomical Society of the Pacific Conference Series, 367 (2003).
  * [8] Iwai, K., Kubo, Y., Ishibashi, H., Naoi, T., Harada, K., Ema, K., Hayashi, Y., and Chikahiro, Y., “OCTAD-S: digital fast Fourier transform spectrometers by FPGA,” Earth, Planets, and Space 69, 95 (July 2017).
  * [9] Hashimoto, I., Chiba, M., Okada, N., Ogawa, H., Kawabe, R., Minamidani, M., Tamura, Y., and Kimura, K., “Wind- and Operation-Induced Vibration Measurements of the Main Reflector of the Nobeyama 45 m Radio Telescope,” Journal of Vibration Engineering & Technologies 8, 909 (2020).


compat=1.0.0 compat=1.0.0

# Probing a quantum spin liquid with equilibrium and non-equilibrium hole
dynamics

J. H. Nyhegn  Center for Complex Quantum Systems, Department of Physics and
Astronomy, Aarhus University, Ny Munkegade, DK-8000 Aarhus C, Denmark  K.
Knakkergaard Nielsen  Max-Planck Institute for Quantum Optics, Hans-
Kopfermann-Str. 1, D-85748 Garching, Germany  G. M. Bruun  Center for Complex
Quantum Systems, Department of Physics and Astronomy, Aarhus University, Ny
Munkegade, DK-8000 Aarhus C, Denmark

###### Abstract

The properties and experimental identification of quantum spin liquids (QSL)
remains an important topic with many fundamental questions. Here, we explore
the dynamics of a single charge dopant (hole) in a $t$-$J_{1}$-$J_{2}$ model
on a square lattice, which realises a gapless $\mathbb{Z}_{2}$ QSL at half-
filling. Using a field theory approach based on the parton construction, which
includes an infinite number of scatterings between the low-energy
quasiparticle excitations of the QSL via a self-consistent Born approximation,
we calculate both the equilibrium and non-equilibrium properties of the hole
for weak and strong interactions. Quasiparticle branches as well as string-
like excitations of the hole are identified, and we furthermore explore the
time-dependent spreading of a hole throughout the QSL after it has been
injected at a given lattice site. The final ballistic expension speed is shown
to exhibit a non-monotonic behavior as a quantum phase transition between an
antiferromagnetic (AFM) and the QSL phase is crossed, which is caused by a
qualitative change in the fundamental kinematics of the interactions between
the hole and the surrounding spins. Our results demonstrate how charge dopants
can be used as a quantum probe for QSLs and are directly relevant to optical
lattice experiments with single site resolution, as well as the rapidly
increasingly number of new quanum materials such as transition metal
dichacogenides.

##  I Introduction

The fundamental nature of QSLs continues to attract great interest since these
phases were first identified more than fifty years ago [1]. Quantum spin
liquid phases such as resonating valence bond states are typically described
as a fluid of strongly correlated spin singlets. The liquid or resonating
character of the ground state means that it consists of a macroscopically
large superposition of spin singlets, which can give rise to exotic properties
such as long-range quantum entanglement, quasiparticle excitations with
nontrivial statistics, and spin-charge separation [2, 3, 4, 5, 6]. Thanks to
sophisticated analytical and numerical methods, our understanding of QSLs
keeps improving [7, 8, 9, 10, 11, 12]. There is also great progress on the
experimental side, where a spin liquid phase has recently been realised using
the new powerful quantum simulation platform based on Rydberg atoms in optical
tweezer arrays [13]. There are furthermore promising prospects of realising
spin liquids in optical lattice experiments [14, 15, 16, 17, 18]. Finally, the
number of strongly correlated two-dimensional (2D) materials is growing
including atomically thin of transition metal dichacogenide bilayers, which
realise highly tunable moiré lattice systems where spin liquids are predicted
to emerge [19, 20, 21].

Figure 1: (a) As a holon (orange circle) moves through a QSL, it distorts the
spin singlets in its surroundings (blue ovals) leaving unpaired spins in its
wake (arrows). These processes give rise to quasiparticle states as well as
string like excitations of the holon. (b) The velocity of a hole spreading
ballistically as a quasiparticle after being injected at a given lattice site.
The expansion velocity increases with $J_{2}/J_{1}$ for $J_{2}/J_{1}<0.46$
(vertical line) where the ground state is an AFM. For $J_{2}/J_{1}>0.46$, the
ground state is a $\mathbb{Z}_{2}$ QSL and the expansion velocity is much
smaller. This gives rise to a non-monotonic behavior of the velocity with a
maximum likely to be close to the quantum phase transition point
$J_{2}/J_{1}=0.46$.

A feature of QSLs is that they elude local magnetic order down to zero
temperature [2]. This invalidates standard approaches such as spin-wave theory
in handling these strongly interacting phases. Except for the Kitaev spin
liquid offering an exact solution [22], this means that numerical tools such
as variational Monte Carlo (VMC) and Matrix Product States (MPSs) have played
an important role exploring QSLs [23, 10, 12, 24, 25]. Morever, a smoking gun
observation of a QSL remains a great challenge. Recent MPS studies suggest
that the dynamical response from inserting charge dopants can be used as
probes for these phases [26, 27]. The remarkable single site resolution
possible in optical lattices has moreover enabled the real space observation
of magnetic polarons formed by a dopant in an AF [28, 29, 17, 18].
Understanding the dynamics of dopants in spin liquids is however challenging
and has mostly been approached phenomenologically or at a mean-field level
[30, 31, 5, 32, 33]. In these studies, the focus has furthermore primarily
been on the interactions with the gauge fluctuations.

These experimental advances makes it timely to explore the dynamics of a
dopant in a QSL. In this paper, we focus on a hole in a so-called
$\mathbb{Z}_{2}$ QSL realized in a $t$-$J_{1}$-$J_{2}$ model at half-filling.
We develop a field theoretical approach for analyzing this, which is capable
of describing the strong coupling regime by including an infinite number of
scatterings between the hole and its surrounding spins via a SCBA. Using this,
we identify quasiparticles as well as string-like excitations of the hole. The
resulting spectral function is shown to be quite different from that of a hole
in an AFM, which indicates that it can be used to detect the presence of a
QSL. We furthermore explore the non-equilibrium spreading of a hole throughout
the lattice ensuing its injection at a specific lattice site. For long times,
the hole is shown to expand ballistically as a quasiparticle whose speed
depends non-monotonically on the ratio $J_{2}/J_{1}$, with a maximum likely
close to the quantum phase transition between the AFM and the QSL phases. This
non-trivial dependence provides another promising way to detect the QSL.

The paper is organized as follows. In Sec. II, we introduce the
$t$-$J_{1}$-$J_{2}$ model. We then discuss the parton construction and our
field theoretical approach. The spectral properties of a holon are discussed
in Sec. III. In Sec. IV, we explore the spectral properties of a physical hole
as well as its non-equilibrium expansion dynamics. By comparing with the same
properties of a hole in an AFM, we discuss in Sec. V two ways where the hole
can be used as a quantum probe of a QSL. We end with conclusions and outlook
in Sec. VI.

##  II Model and field theory

We study the behavior of a single dopant in a $\mathbb{Z}_{2}$ QSL in the
paradigmatic $t$–$J$ model with both nearest (NN) and next-nearest neighbor
(NNN) interactions

$\hat{H}=\hat{H}_{t}+\hat{H}_{J}.$ (1)

Here

$\displaystyle\hat{H}_{t}=$ $\displaystyle-t_{1}\sum_{\braket{{\bf i},{\bf
j}},\sigma}\tilde{c}^{\dagger}_{{\bf i},\sigma}\tilde{c}_{{\bf
j},\sigma}-t_{2}\sum_{\braket{\braket{{\bf i},{\bf
j}}},\sigma}\tilde{c}^{\dagger}_{{\bf i},\sigma}\tilde{c}_{{\bf
j},\sigma}+{\rm H.c.}$ (2)

is the hopping Hamiltonian with NN and NNN parameters $t_{1},t_{2}$, and

$\displaystyle\\!\\!\\!\hat{H}_{J}=$ $\displaystyle\ J_{1}\sum_{\braket{{\bf
i}{\bf j}}}\left[\hat{\mathbf{S}}_{{\bf i}}\cdot\hat{\mathbf{S}}_{{\bf
j}}-\frac{1}{4}\hat{n}_{{\bf i}}\hat{n}_{{\bf
j}}\right]+J_{2}\\!\\!\sum_{\braket{\braket{{\bf i}{\bf
j}}}}\\!\left[\hat{\mathbf{S}}_{{\bf i}}\cdot\hat{\mathbf{S}}_{{\bf
i}}-\frac{1}{4}\hat{n}_{{\bf i}}\hat{n}_{{\bf j}}\right]\\!$ (3)

describes AFM Heisenberg type NN and NNN spin couplings $J_{1},J_{2}>0$. The
restricted fermionic operators $\tilde{c}^{\dagger}_{{\bf
i},\sigma}=\hat{c}^{\dagger}_{{\bf i},\sigma}(1-\hat{n}_{{\bf i}})$ ensure no
double occupancy in terms of the onsite density operator $\hat{n}_{{\bf
i}}=\sum_{\sigma}\hat{c}^{\dagger}_{{\bf i},\sigma}\hat{c}_{{\bf i},\sigma}$.
Here, $\hat{c}_{{\bf i}\sigma}$ removes a fermion with spin $\sigma$ at site
${\bf i}$ and $\hat{{\bf S}}_{\bf
i}=\frac{1}{2}\sum_{\sigma,\sigma^{\prime}}\hat{c}^{\dagger}_{{\bf
i},\sigma}\boldsymbol{\sigma}_{\sigma\sigma^{\prime}}\hat{c}_{{\bf
i},\sigma^{\prime}}$. The spin coupling in Eq. (3) naturally arises as an
effective low-energy description of super-exchange processes in the Fermi-
Hubbard model for large onsite repulsion $U$ close to half-filling, where
$J_{1}=4t_{1}^{2}/U$ and $J_{2}=4t_{2}^{2}/U$. As this is relevant for
strongly correlated 2D materials and experiments with atoms in optical
lattices, we will in the following uphold these relations such that
$J_{2}/J_{1}=(t_{2}/t_{1})^{2}$.

At half-filling, the Heisenberg model supports four distinct phases [34, 35,
36, 37]. With no NNN interaction, the model exhibits antiferromagnetic (AFM)
order. Turning on a non-zero $J_{2}$, we initiate a competition that
frustrates the magnetic order, and by increasing $J_{2}/J_{1}$ the system
first enters a gapless $\mathbb{Z}_{2}$ QSL at $J_{2}/J_{1}\sim 0.46$, then a
valence bond solid (VBS) phase at $J_{2}/J_{1}\sim 0.55$, and at last the
magnetically ordered columnar phase at $J_{2}/J_{1}\sim 0.6$. To explore the
QSL phase, one, thus, needs to be able to achieve NNN hoppings
$|t_{2}|\simeq|t_{1}|/\sqrt{2}$. While NNN hoppings in standard optical
lattice experiments are inherently much smaller than the NN hoppings, optical
superlattices have been shown [38] to allow the desired control over
$t_{2}/t_{1}$. Focusing on the QSL phase, we now turn to the parton
construction.

###  II.1 Parton construction

Following Refs. [39, 40, 5], we decompose the fermions into partons as

$\displaystyle\hat{c}_{{\bf i},\sigma}=\hat{b}^{\dagger}_{{\bf
i}}\hat{f}_{{\bf i},\sigma},$ (4)

where $\hat{b}_{{\bf i}}$ removes a bosonic holon at site ${\bf i}$ describing
the charge degree of freedom, and $\hat{f}_{{\bf i},\sigma}$ removes a
fermionic spinon describing the spin degree of freedom. Using this parton
decomposition in Eq. (3), the mean-field theory most accurately describing the
ground state for $J_{2}\simeq 0.5J_{1}$ at half-filling with no holons is
found to be a gapless $\mathbb{Z}_{2}$ spin liquid [35, 10, 41, 42, 43]. From
VMC results, the optimal mean-field parton Hamiltonian was found to be

$\displaystyle\hat{H}_{J}=\sum_{{\bf k},\sigma}\epsilon_{{\bf
k}}\hat{f}^{\dagger}_{{\bf k},\sigma}\hat{f}_{{\bf k},\sigma}+\sum_{{\bf
k}}\left(\Delta_{{\bf k}}\hat{f}^{\dagger}_{{\bf
k},\uparrow}\hat{f}^{\dagger}_{-{\bf k},\downarrow}+\text{h.c.}\right),$ (5)

with the BCS-type pairing functions [8, 43, 10]

$\displaystyle\epsilon_{{\bf k}}$ $\displaystyle=\epsilon(\cos k_{x}+\cos
k_{y})$ $\displaystyle\Delta_{{\bf k}}$ $\displaystyle=\Delta_{1}(\cos
k_{x}-\cos k_{y})+\Delta_{2}\sin{(2k_{x})}\sin{(2k_{y})},\\!\\!$ (6)

where we have taken the lattice constant to unity. Within the theory of the
projective-symmetry group, this spin liquid is denoted Z2Azz13 [8]. Instead of
determining the parameters in Eq. (6) via their self-consistent mean-field
relations to the Heisenberg parameters $J_{1}$ and $J_{2}$ [44], we use unless
otherwise stated the values $\Delta_{1}=1.1\epsilon$ and
$\Delta_{2}=1.8\epsilon$, found as the optimal parameters in the VMC study
[43] for $J_{2}/J_{1}=1/2$. Note that Eq. (5) is defined in an enlarged
Hilbert space where unphysical doubly occupied sites are allowed keeping only
the average number of particles per site fixed at unity. These unphysical
states can be projected out by using the Gutzwiller operator
$\hat{\mathcal{P}}=\prod_{{\bf i}}(\hat{n}_{{\bf i},\downarrow}-\hat{n}_{{\bf
i},\uparrow})^{2}$ [44, 24, 10, 37]. It has however been shown that the the
mean-field description given by Eq. (5) without a subsequent Gutzwiller
projection accurately captures the low-energy excitations of the spin liquid
state of the Heisenberg Hamiltonian $\hat{H}_{J}$ taking $\epsilon\simeq
J_{1}/3$ [45]. Not having to perform a Gutzwiller projection dramatically
simplifies the description, and we therefore use this mean-field approach in
the following to describe the QSL. Diagonalizing Eq. (5) then yields

$\displaystyle\hat{H}_{J}=$ $\displaystyle\sum_{{\bf
k},\sigma}\omega^{s}_{{\bf k}}\hat{\gamma}^{\dagger}_{{\bf
k},\sigma}\hat{\gamma}_{{\bf k},\sigma},$ (7)

with the spinon spectrum

$\displaystyle\omega^{s}_{{\bf k}}=\sqrt{\epsilon_{{\bf k}}^{2}+\Delta_{{\bf
k}}^{2}},$ (8)

shown in Fig. 2(a). The BCS-type ground state is defined by
$\hat{\gamma}_{{\bf k},\sigma}\ket{\text{GS}}=0$. Note that the Cooper pairs
in this mean-field theory correspond to the spin singlets.

Figure 2: (a) Spinon spectrum given by Eq. (8) with $\Delta_{1}=1.1\epsilon$
and $\Delta_{2}=1.8\epsilon$. (b) Bare holon spectrum given by Eq. (11).

The hopping term $\hat{H}_{t}$ comes into play when we add a single dopant to
the system. Using the parton construction yields

$\displaystyle\hat{H}_{t}=-t_{1}\sum_{\left\langle{\bf i}{\bf
j}\right\rangle}\hat{b}_{{\bf i}}\hat{f}^{\dagger}_{{\bf
i},\sigma}\hat{f}_{{\bf j},\sigma}\hat{b}^{\dagger}_{{\bf
j}}-t_{2}\sum_{\left\langle\left\langle{\bf i}{\bf
j}\right\rangle\right\rangle}\hat{b}_{{\bf i}}\hat{f}^{\dagger}_{{\bf
i},\sigma}\hat{f}_{{\bf j},\sigma}\hat{b}^{\dagger}_{{\bf j}},$ (9)

where we have removed the projection operators since we only consider a single
dopant and allow multiple spinons on a given site [5, 39, 46]. Expressing this
in terms of the spinon operators $\hat{\gamma}^{\dagger}_{{\bf k},\sigma}$
diagonalising $\hat{H}_{J}$ giving Eq. (7), we can write Eq. (9) as

$\displaystyle\hat{H}_{t}=$ $\displaystyle\sum_{{\bf k}}\omega^{b}_{{\bf
k}}\hat{b}^{\dagger}_{{\bf k}}\hat{b}_{{\bf k}}+\sum_{{\bf k},{\bf
k}^{\prime},{\bf q},\sigma}h_{{\bf k},{\bf k}^{\prime},{\bf
q}}\hat{b}^{\dagger}_{{\bf k}-{\bf q}}\hat{b}_{{\bf
k}}\hat{\gamma}^{\dagger}_{{\bf k}^{\prime},\sigma}\hat{\gamma}_{{\bf
k}^{\prime}-{\bf q},\sigma}$ $\displaystyle+$ $\displaystyle\sum_{{\bf k},{\bf
k}^{\prime},{\bf q}}g_{{\bf k},{\bf k}^{\prime},{\bf
q}}\left(\hat{b}^{\dagger}_{{\bf k}-{\bf q}}\hat{b}_{{\bf
k}}\hat{\gamma}^{\dagger}_{{\bf
k}^{\prime},\uparrow}\hat{\gamma}^{\dagger}_{-{\bf k}^{\prime}+{\bf
q},\downarrow}+\text{h.c.}\right).$ (10)

Here

$\omega^{b}_{\bf k}=2\sum_{{\bf k}^{\prime}}\Gamma_{{\bf k}+{\bf
k}^{\prime}}v_{{\bf k}^{\prime}}^{2}$ (11)

is the dispersion of a bare holon shown in Fig. 2(b), and

$\displaystyle g_{{\bf k},{\bf k}^{\prime},{\bf q}}$
$\displaystyle=-\left[\Gamma_{{\bf k}-{\bf k}^{\prime}}u_{{\bf
k}^{\prime}}v_{{\bf k}^{\prime}-{\bf q}}+\Gamma_{{\bf k}+{\bf k}^{\prime}-{\bf
q}}v_{{\bf k}^{\prime}}u_{{\bf k}^{\prime}-{\bf q}}\right]$ $\displaystyle
h_{{\bf k},{\bf k}^{\prime},{\bf q}}$ $\displaystyle=\Gamma_{{\bf k}-{\bf
k}^{\prime}}u_{{\bf k}^{\prime}}u_{{\bf k}^{\prime}-{\bf q}}-\Gamma_{{\bf
k}+{\bf k}^{\prime}-{\bf q}}v_{{\bf k}^{\prime}}v_{{\bf k}^{\prime}-{\bf q}}$
(12)

are the vertex functions with

$\displaystyle\Gamma_{\bf k}=$
$\displaystyle\frac{2}{N}(t_{1}\left[\cos{k_{x}}+\cos{k_{y}}\right]$
$\displaystyle+$ $\displaystyle
t_{2}\left[\cos{(k_{x}+k_{y})}+\cos{(k_{y}-k_{y})}\right])$ (13)

coming from NN and NNN structure factors. The number of lattice sites is $N$,
and the coherence factors are $u_{{\bf k}}=\sqrt{(1+\epsilon_{{\bf
k}}/\omega^{s}_{{\bf k}})/2}$ and $v_{{\bf k}}=\text{sgn}\Delta_{{\bf
k}}\sqrt{(1-\epsilon_{{\bf k}}/\omega^{s}_{{\bf k}})/2}$. There are two kinds
of interaction vertices between the holon and the spinons in Eq. (10):
$h_{{\bf k},{\bf k}^{\prime},{\bf q}}$ describes the scattering of a holon on
a spinon whereas $g_{{\bf k},{\bf k}^{\prime},{\bf q}}$ describes the
emission/absorption of two spinons by a holon.

With Eqs. (7) and (10), we have arrived at our final Hamiltonian describing a
dopant in a $\mathbb{Z}_{2}$ QSL. Staying true to the Fermi-Hubbard relation
$J_{i}=4t_{i}^{2}/U$, we use $t_{2}/t_{1}=-\sqrt{J_{2}/J_{1}}=-1/\sqrt{2}$
throughout corresponding to the case of a single hole doping
($t_{2}/t_{1}=1/\sqrt{2}$ for a doublon [47]). We also keep the ratio
$\Delta_{2}/\Delta_{1}=1.8/1.1$ fixed and $\Delta_{1}=1.1\epsilon$ unless
otherwise stated, which leaves $\epsilon/t_{1}$ as a free parameter in our
theory. Strong interaction with the holon creating many spinon excitations
corresponds to $\epsilon/t_{1}\ll 1$ whereas $\epsilon/t_{1}\gg 1$ corresponds
to weak interaction.

###  II.2 Quantum field theory

We will now describe our quantum field theory approach to describe the motion
of the hole in the QSL. In the parton construction, a first natural step is to
analyse the motion of the holon quasiparticle [30, 33] as described with the
retarded holon Green’s function

$G({\bf k},\tau)=-i\theta(\tau)\langle\\{\hat{b}_{{\bf
k}}(\tau),\hat{b}^{\dagger}_{{\bf k}}(0)\\}\rangle,$ (14)

where $\hat{A}(\tau)=\exp(i\hat{H}\tau)\hat{A}\exp(-i\hat{H}\tau)$ is the
operator at time $\tau$ in the Heisenberg picture and
$\\{\hat{A},\hat{B}\\}=\hat{A}\hat{B}+\hat{B}\hat{A}$. In frequency/momentum
space, a formal solution for the Green’s function is

$\displaystyle G({\bf k},\omega)=\frac{1}{\omega-\omega^{b}_{{\bf
k}}-\Sigma({\bf k},\omega)+i0^{+}},$ (15)

where $\Sigma({\bf k},\omega)$ is the holon self-energy. In order to calculate
the self-energy, one in general needs to resort to approximations. Here, we
will apply the self-consistent Born approximation (SCBA) [46], which has
proven to be very accurate for describing a mobile hole in an AFM even for
strong interactions when compared to Monte-Carlo calculations and exact
diagonalisation [48, 49]. Remarkably, this holds even for describing non-
equilibrium hole dynamics [50]. In the present case with the two kinds of
spinon-hole interaction vertices given by Eq. (12), the SCBA amounts to
including the two types of diagrams for the holon self-energy shown in Fig. 3.
These two diagrams give the give the exact solution of the self-energy to
third order in the hopping amplitudes, and by using the self-consistent holon
Green’s function we include an infinite class of diagrams in order to describe
the limit $\epsilon/t_{1}\ll 1$ of strong interactions. The self-consistent
diagrams in Fig. 3 give the self-energy

$\displaystyle\Sigma({\bf k},\omega)=\sum_{{\bf k}^{\prime},{\bf
q}_{1}}g_{{\bf k},{\bf k}^{\prime},{\bf q}_{1}}^{2}G({\bf
K}_{1},\omega-\omega^{s}_{{\bf k}^{\prime}}-\omega^{s}_{-{\bf k}^{\prime}+{\bf
q}_{1}})$ $\displaystyle-2\sum_{{\bf q}_{2}}g_{{\bf k},{\bf k}^{\prime},{\bf
q}_{1}}G({\bf K}_{1},\omega-\omega^{s}_{{\bf k}^{\prime}}-\omega^{s}_{-{\bf
k}^{\prime}+{\bf q}_{1}})h_{{\bf k}-{\bf q}_{1},{\bf k}^{\prime}+{\bf
q}_{2},{\bf q}_{2}}$ $\displaystyle\times G({\bf
K}_{2},\omega-\omega^{s}_{{\bf k}^{\prime}+{\bf q}_{2}}-\omega^{s}_{-{\bf
k}^{\prime}+{\bf q}_{1}})g_{{\bf k},{\bf k}^{\prime}+{\bf q}_{2},{\bf
q}_{1}+{\bf q}_{2}},\\!\\!$ (16)

with ${\bf K}_{n}={\bf k}-{\bf q}_{1}-...-{\bf q}_{n}$. Thus, our field
theoretical approach amounts to calculating the holon Green’s function by
solving Eqs. (15) and (16) self-consistently, which includes an infinite
number of spinon excitations created by the holon.

Figure 3: Feynman diagrams included in the extended SCBA for the holon self-
energy. Double wavy line are the full holon propagator and the single lines
are spinon propagators. The blue and red dots are the $g$\- and $h$-vertices
in Eq. (12) respectively.

##  III Holon properties

We now describe the properties of a holon dynamics captured by its spectral
function

$\displaystyle A({\bf k},\omega)=-2\text{Im}\left[G({\bf k},\omega)\right],$
(17)

which gives the energy spectrum of a holon injected with momentum ${\bf k}$
into the spin liquid. As detailed in App. A, it turns out that the second
diagram in Fig. 3 has only small effects on the self-energy and since it is
computationally quite demanding to evaluate, it is omitted in the following.

Initiating the investigation in the weak coupling regime, we plot in Fig. 4(a)
the holon spectral function for $\epsilon/t_{1}=1$ and two different momenta
obtained by numerically solving Eqs. (15) and (16) self-consistently on a
$32\times 32$ lattice. We see that the spectral function exhibits a clear
quasiparticle peak for both momenta. Figure 4(b) shows that the quasiparticle
is well-defined for all crystal momenta along the specified path in the
Brillouin zone. The energy of this quasiparticle is due to interactions with
the spinons shifted by a constant relative to the bare holon dispersion given
by the green dashed line. We note that for all results, the spectral function
obeys the sum rule $\int d\omega A_{h}({\bf p},\omega)=2\pi$ with less than 1%
deviation confirming the accuracy of our numerics.

Figure 4: (a) Spectral function of the holon for two different crystal momenta
and the interaction strength $\epsilon/t_{1}=1$. (b) Spectral function for
crystal momenta along the path in the Brillouin zone shown in Fig. 2(b). The
green dashed line is the bare holon dispersion given by Eq. (11).

Figure 5(a) shows that for intermediate interactions $\epsilon/t_{1}=0.5$, the
quasiparticle peak in holon spectral function has lost significant spectral
weight to a many-body continuum. Indeed, the quasiparticle is completely
washed out for momenta around ${\bf k}=(0,0)$, which can be attributed to the
emissions of spinon pairs as described by the first Feynman diagram in Fig. 3.
The on-shell energies of these emissions form a continuum between
$\text{min/max}\\{\omega^{b}_{{\bf k}-{\bf q}}+\omega^{s}_{{\bf
k}^{\prime}}+\omega^{s}_{-{\bf k}^{\prime}+{\bf q}}\\}$, shown as orange
dashed lines in Fig. 5(a). We clearly see that the quasiparticle is heavily
damped when its energy is inside this continuum.

Figure 5: Spectral function of the holon along the same path in the Brillouin
zone as in Fig. 4, but now for intermediate interaction strength
$\epsilon/t_{1}=0.5$ (a) and strong interaction strength $\epsilon/t_{1}=0.1$
(b). The dashed green line in (a) is again the bare holon dispersion given by
Eq. (11) whereas the orange lines indicate the boundaries of the spinon-holon
many-body continuum. The red dashed line in (b) is a fit to the holon
dispersion, see App. E.

Increasing the interacting strength further to $\epsilon/t_{1}=0.1$, Fig. 5(b)
shows that the quasiparticle peak is damped even more due to decay into spinon
modes as expected. We find that the width of the quasiparticle band is
approximately given by $3.9\epsilon/t_{1}$, i.e. increasing the interaction
strength decreases the bandwidth. This reflects that increasing the dressing
of the holon by spinon modes results in a heavier quasiparticle. Such
significant dressing cannot be captured by mean-field approaches, which in
stark contrast predict a holon quasiparticle with a dispersion scaling as
$\sim t_{1}$ [30, 51]. Our results are on the other hand consistent with MPS
calculations for a chiral spin liquid predicting similar features for the
physical hole [26]. We will analyse the properties of a physical hole in Sec.
IV.

###  III.1 String excitations

Intriguingly, one can also see several additional bands at higher energy
energies in Fig. 5(b). To investigate the origin of these bands, we plot Fig.
6(a) the spectral function for momentum ${\bf p}=(\pi,\pi)$ as a function of
$\epsilon/t_{1}$, clearly demonstrating that more excited states emerge with
increasing interaction strength. Figure 6(b) shows that the energy of the
ground state and these excited states scale as $(\epsilon/t_{1})^{2/3}$ for
low $\epsilon/t_{1}$. For a hole in an AFM, such a scaling is interpreted as a
smoking gun for the presence of so-called string excitations [52]. The basic
mechanism is that the hole leaves a trace of frustrated spins in its wake as
it moves through the AFM, which means that it experiences a linear potential.
This gives rise to string states with the characteristic energy scaling in the
continuum limit. Based on this, we conjecture that the bands in Fig. 5(b) also
arise from string excitations of the holon. While there is no local magnetic
order in a spin liquid, the holon still distorts its surroundings by creating
pairs of spinon excitations as it moves though the lattice. This is explicitly
evident in Eq. (10) and Fig. 3, and is illustrated in Fig. 1(a). The result is
that the holon experiences a linear potential in real space leading to the
excited string states visible in Fig. 5(b).

Figure 6: (a) The holon spectral function for crystal momentum $(\pi,\pi)$ as
a function of the interaction strength $\epsilon/t_{1}$. (b) Energy of the
three lowest bands in (a) replotted as a function of $(\epsilon/t_{1})^{2/3}$.

To support this hypothesis further, we plot in Fig. 6 the energies of the
excited states as a function of the spinon excitation energy averaged over the
Brillouin zone $\omega^{s}_{avg}=\sum_{{\bf k}}\omega^{s}_{{\bf k}}/N$. The
different colors correspond to calculations with different ratios of
$\Delta_{1}/\epsilon$. Figure 6 shows that these different calculations
collapse on a single line corresponding to the excitation energy scaling as
$(\omega^{s}_{avg}/t_{1})^{2/3}$. This confirms that the low-energy states can
be understood as string excitations coming from the holon breaking spin
singlets into individual spinons so that it experiences a linear potential
with a strength $\propto\omega^{s}_{avg}$. In App. B, we further investigate
the nature of these string excitations and find that they are suppressed if
one increases the average length of the spin singlets.

Figure 7: Excitation energy $E_{1}-E_{0}$ of the lowest string band for
momentum $(\pi,\pi)$ as a function of $(\omega^{s}_{avg}/t_{1})^{2/3}$.
$\omega^{s}_{avg}$ is the average energy of a spinon, and $E_{0}$ is the
energy of excitation at $\epsilon/t_{1}=0$.

##  IV Properties of physical holes

We are now ready to explore the properties of the physical holes, which are
the ones observable in experiments. To do this, we consider the hole Green’s
function

$\displaystyle G_{h\sigma}({\bf
k},\tau)=-i\theta(\tau)\langle\\{\hat{c}^{\dagger}_{{\bf
k},\sigma}(\tau),\hat{c}_{{\bf k},\sigma}(0)\\}\rangle$
$\displaystyle=-i\theta(\tau)\langle\hat{c}^{\dagger}_{{\bf
k},\sigma}(\tau)\hat{c}_{{\bf k},\sigma}(0)\rangle=\theta(\tau)G_{h}^{>}({\bf
k},\tau).$ (18)

The second identify follows from the fact that only the single hole injected
is present in the lattice, and it is crucial for enabling us to describe the
non-equilibrium dynamics without having to use more cumbersome Keldysh
techniques, as will be exploited in Sec. IV.2.

Using the parton construction, the hole Green’s function becomes a two-body
Green’s function

$\displaystyle G_{h}({\bf k},\tau)=-\frac{1}{N}\\!\sum_{{\bf q}_{1},{\bf
q}_{2}}\\!\left\langle\hat{f}^{\dagger}_{{\bf k}+{\bf
q}_{2},\sigma}(\tau)\hat{b}_{{\bf q}_{2}}(\tau)\hat{b}^{\dagger}_{{\bf
q}_{1}}\hat{f}_{{\bf k}+{\bf q}_{1},\sigma}\right\rangle,\\!$ (19)

describing a holon-spinon ”particle-hole” bubble. Here and in the following,
we suppress the spin index on the Green’s function since it does not depend on
the direction of the spin due to the SU$(2)$ spin symmetry of the Fermi-
Hubbard model. Writing it in terms of the rotated spinons gives

$\displaystyle G_{h}({\bf k},\tau)=\ $ $\displaystyle\frac{-1}{N}\sum_{{\bf
q}_{1},{\bf q}_{2}}v_{{\bf k}+{\bf q}_{1}}v_{{\bf k}+{\bf q}_{2}}$
$\displaystyle\left\langle\hat{\gamma}_{-{\bf k}-{\bf
q}_{2},\downarrow}(\tau)\hat{b}_{{\bf q}_{2}}(\tau)\hat{b}^{\dagger}_{{\bf
q}_{1}}\hat{\gamma}^{\dagger}_{-{\bf p}-{\bf q}_{1},\downarrow}\right\rangle.$
(20)

Since $G_{h}({\bf k},\tau)$ is now transformed into a two-body Green’s
function due to the parton construction, it is in general very challenging to
calculate. Here, we resort to approximating it with the lowest-order diagram

(21)

which describes the uncorrelated motion of a spinon and a holon. We find that
this approximation obeys the frequency sum rule for the resulting hole
spectral function within $1\%$ and we expect it to describe the most important
physics for the present purposes. It does not capture exotic features such as
the possibility of spinon-holon bound states, which however likely will have a
small spectral weight due to their small Franck-Condon overlap with a free
hole. Such bound states will be examined in a future work.

###  IV.1 Spectral properties

To explore the spectral properties of the physical hole, measurable with ARPES
[26, 27], we plot in Fig. 8(a) its spectral function

$\displaystyle A_{h}({\bf k},\omega)=-2\text{Im}\left[G_{h}({\bf
k},\omega)\right]$ (22)

calculated numerically from Eqs. (21)-(22) for the case of weak interactions
$\epsilon/t_{1}=1$. We clearly see a quasiparticle like band, which however is
significantly broadened compared to the corresponding holon quasiparticle band
in Fig. 4(b). It follows from Eq. (21) that the hole spectral function
$A_{h}({\bf k},\omega)$ is non-zero for energies satisfying
$\omega\approx\omega^{\rm QP}_{{\bf q}}+\omega^{s}_{{\bf k}+{\bf q}}$. The
dispersion of the hole band is in fact to a very good approximation given by

$\displaystyle\omega^{h}_{{\bf k}}=\omega^{\rm QP}_{{\bf
k}+\boldsymbol{\pi}}+\omega^{s}_{\boldsymbol{\pi}}$ (23)

with $\boldsymbol{\pi}=(\pi,\pi)$ as shown by the green dashed line in Fig.
8(a). This can be understood from the fact that the spinon spectrum shown in
Fig. 2(a) has two Van Hove singularities at the center and corners of the
Brillouin zone. While the coherence factor $v^{2}_{{\bf k}+{\bf q}}$ in Eq.
(21) vanishes at the origin suppressing its contribution to the hole spectral
function, the Van Hove singularity at $\boldsymbol{\pi}$ gives rise to a large
contribution to the hole spectral function, which in turn gives a band with a
dispersion given by Eq. (23). Similar features were reported in Ref. [27],
where a hole in a Kitaev spin liquid was investigated using MPS methods.

Figure 8:  (a) Spectral function of a physical hole along the same path in the
Brillouin zone as in Fig. 4 for interaction strength $\epsilon/t_{1}=1$ and
$\Delta_{1}=1.1\epsilon$. The green dashed line is Eq. (23). In panel (b) the
pairing field is taken to be ten times smaller with the top green dashed line
showing Eq. (23) and the bottom green dashed line $\omega=\omega^{\rm
QP}_{{\bf k}+(0,\pi)}+\omega^{s}_{(0,\pi)}$.

In order to explore further how the spinons affect the properties of the hole,
we plot in Fig. 8(b) its spectral function where we have used the same
parameters as in Fig. 8(a) except for a smaller pairing potential
$\Delta_{{\bf k}}\rightarrow 0.1\Delta_{{\bf k}}$. Since the length of the
singlets is determined by the BCS coherence length $\propto 1/\Delta$, this
makes them longer, see App. B. Figure 8(b) shows that spectral function of the
physical hole in turn becomes much broader no with clearly visible band except
for small momenta, reflecting the increased dressing of the underlying holon
by spinons. The many-body continuum of the hole spectral function is
approximately bounded from below by the energy $\omega^{\rm QP}_{{\bf
k}+(0,\pi)}+\omega^{s}_{(0,\pi)}$, which comes form a new Van Hove singularity
in the spinon spectrum at ${\bf k}=(0,\pi)$. These results show that the hole
spectrum is broadened with less visible energy bands for spin liquids with
longer singlets and a softer spinon spectrum.

In Sec. III.1, we saw that the holon spectrum exhibits string excitations and
we now investigate whether they are present also for the physical hole. Figure
9 shows the hole spectral function for strong interactions
$\epsilon/t_{1}=0.1$ again taking $\Delta_{1}/\epsilon=1.8$ where the spectral
function of the holon in Fig. 5(b) exhibits clear string excitation bands.
While the spectrum as expected is quite broad, there is still a band visible
with an energy given by Eq. (23) shown by the green dashed line. Importantly,
we also see distinct bands at higher energies and as analysed further in App.
D, they are a direct consequence of the string excitations of the holon.
Hence, we conclude that the hole exhibits string excitations when it interacts
strongly with a surrounding spin liquids. We note that similar peaks at
excited energies have previously been interpreted as bound states between a
spinon and a holon [30, 31, 51].

Figure 9: (a) Spectral function of a physical hole for momentum ${\bf p}={\bf
0}$ and strong coupling $\epsilon/t_{1}=0.1$. The arrow indicates a string
excitation. In (b), the spectral function is plotted along the same path in
the Brillouin zone as in Fig. 4. The green dashed line shows Eq. (23).

###  IV.2 Real space non-equilibrium dynamics

In a recent ground breaking optical lattice experiment, the non-equilibrium
motion of a hole after its release from a given lattice site in an AFM was
explored in real space with single site resolution [29]. In light of this
experiment and the prospects of realising spin liquids in optical lattices
[14, 15, 16, 17, 18], we therefore now turn to investigating the real-space
dynamics following a sudden insertion of a hole.

Concretely, we calculate the hole Green’s function $G_{h}^{>}({\bf
r},\tau)=-i\langle\hat{c}^{\dagger}_{{\bf
r},\sigma}(\tau)\hat{c}_{{\mathbf{0}},\sigma}\rangle$ in real space and time,
which can be obtained from $G_{h}^{>}({\bf k},\omega)$ by Fourier
transforming. As discussed below Eq. (18), this gives us access to the non-
equilibrium dynamics described by the overlap of taking out a bare hole at
site ${\bf r}$ and time $\tau>0$ after it was inserted at the origin at time
$\tau=0$. Note that in order to describe all dynamics, one needs the full
time-dependent many-body wave function which is a much harder problem [50].
However, the bare hole dynamics described by $G_{h}^{>}({\bf r},\tau)$ is
sufficient to explore the long time dynamics as has been successfully
exploited both for magnetic polarons in AFMs [53, 50] and for polarons in
atomic gases [54, 55].

In Fig. 10, we show an example of such calculations. We plot

Figure 10: The expansion of the hole throughout the lattice at times
$\tau=5/t_{1}$ (a) and $\tau=10/t_{1}$ (b) after it was injected at the
origin, as parameterized by $|G_{h}^{>}|$. The interaction strength is
$\epsilon/t_{1}=0.25$.

the the overlap $|G_{h}^{>}|$ in real space for two different times after the
hole has been injected at the origin at time $\tau=0$ for the interaction
strength $\epsilon/t_{1}=0.25$. This clearly shows how the hole spreads out
through the lattice in a dynamical fashion. Note that the accumulated
probability obtained by integrating the signal over the entire lattice is much
larger in Fig. 10(a) than in Fig. 10(b). This reflects that the hole becomes
increasingly dressed with spinons as it moves through the lattice, so that the
total probability of observing a _bare_ hole decreases. The slight symmetry
breaking barely visible in the lower panel is an artifact of the gauge choice
in the mean-field description of the QSL.

In Fig. 11(a), we plot the rms distance of the hole from its initial position
at the origin defined as

$\displaystyle d_{\text{rms}}(\tau)=\sqrt{\frac{\sum_{{\bf
r}}r^{2}|G^{>}_{h}({\bf r},\tau)|^{2})}{\sum_{{\bf r}}|G^{>}_{h}({\bf
r},\tau)|^{2}}}.$ (24)

The denominator takes care of the decreasing probability to observe a bare
hole with time, and it is easy to see that $d_{\text{rms}}$ is gauge
invariant. Figure 11(a) clearly show how the expansion of the hole slows down
with increasing interaction strength $\epsilon/t_{1}$. Physically, the reason
is that the hole becomes increasingly dressed by spinons making it heavier.
Indeed, similar results were found regarding the expansion of a hole in an AFM
both theoretically [56, 53, 50] and experimentally [29].

To quantify this, we plot in Fig. 11(a) the expansion speed
$v_{\mathrm{rms}}=\partial_{\tau}d_{\text{rms}}$ as a function of time for
different interaction strengths. This shows that after some initial non-
trivial dynamics, the hole expands ballistically for long times with a speed
that decreases with increasing interaction strength. Figure 11(b) plots the
final velocity of this ballistic expansion as a function of the interaction
strength, which clearly shows this decrease. We also plot in Fig. 11(b) the
bandwidth of the lowest band of the hole spectral function obtained from plots
like those in Sec. IV.1. The agreement between these two quantities in Fig.
11(b) demonstrates that for long times, the hole expands as a quasiparticle
with a dispersion determined by the lowest band in its spectral function. As
this band becomes flatter with increasing interaction strength reflecting to
dressing of the bare hole with spinons, it slows down.

Figure 11: (a) The rms distance from its initial position of an expanding hole
defined by Eq. (24) as a function of time for different interaction strengths.
(b) The associated long-time expansion speed as a function of interaction
strength (blue line). The red line gives the bandwidth of the lowest hole
quasiparticle band.  Figure 12: The rms distance as a function of time for the
interaction strength $\epsilon/t_{1}=0.25$ and different pairing strengths
$\Delta_{1}/\epsilon$.

##  V The hole as a quantum probe

A central problem regarding QSLs is how to detect them. Armed with our field
theoretical formalism capable of describing both equilibrium and non-
equilibrium dynamics, we will in this section suggest two ways to use a hole
as a quantum probe for spin liquids.

In order to do this, we compare the spectral as well as expansion properties
of a hole in a QSL with those in an AFM, which is the ground state of the
Heisenberg model $\hat{H}_{J}$ for small $J_{2}$ [34, 35, 41]. To calculate
the properties of a hole in an AFM, we use linear spin wave theory combined
with the SCBA, which has been shown to be remarkably accurate both for the
equilibrium and non-equilibrium properties [46, 57, 48, 58, 59, 60, 61, 49,
62, 50, 53, 63]. We refer the reader to Refs. [62, 50, 53, 63] and App. C for
the technical details of this method.

Let us first compare the spectral properties of a hole in an AFM and in a QSL.
In Fig. 13, we show the spectral function of the hole in an AFM for
$J_{2}/J_{1}=0.15$, $t_{2}/t_{1}=-\sqrt{J_{2}/J_{1}}$, and $J_{1}/t_{1}=0.3$.
Using $J_{1}\simeq 3\epsilon$ for the QSL [45], the latter ratio corresponds
to $\epsilon/t_{1}\sim 0.1$. It is therefore natural to compare Fig. 9(b),
which shows the spectral function of a hole in a QSL with $\epsilon/t_{1}\sim
0.1$ and $J_{2}/J_{1}=1/2$. We see that these two spectral functions are
qualitatively different. Whereas the spectral function in the AFM phase has
several sharp quasiparticle bands corresponding to a ground state
quasiparticle and string excitations, which take up significant spectral
weight, the spectral function in the QSL phase has significantly broader bands
with a many-body continuum taking up much more weight. From this we conclude
that the spectral function of a hole provides a promising way to detect a spin
liquid. Similar results were reported in Ref. [26], where the spectral
function for a hole in a chiral spin liquid was calculated using MPS methods.

Figure 13:  Spectral function of a hole along the same path in the Brillouin
zone as in Fig. 4 for $J_{1}/t_{1}=0.3$, $J_{2}/J_{1}=0.15$, and
$t_{2}/t_{1}=-\sqrt{J_{2}/J_{1}}$, where the system is in an AFM phase.

Next, we compare the expansion dynamics of a hole in an AFM with that in a
spin liquid, which can be measured in cold atom experiments. Our main result
is shown in Fig. 1, which plots the long time ballistic expansion speed
$v_{\text{rms}}$ as a function of the ratio $J_{2}/J_{1}$, where we have
assumed that the spins form an AFM for $J_{2}/J_{1}<0.46$ indicated by the
vertical black line, and a spin liquid for $J_{2}/J_{1}>0.46$ [43, 35, 41, 36,
34]. We see that the expansion speed of the hole increases with $J_{2}/J_{1}$
in the AFM phase. This may seem surprising at first since increasing
$J_{2}/J_{1}$ leads to a softer spin-wave spectrum due to increased magnetic
frustration, from which it would be tempting to conclude that the hole becomes
more dressed and therefore slower in analogy with what happens in the QSL. The
reason this does not occur is that a non-zero value of
$t_{2}=-\sqrt{J_{2}/J_{1}}t_{1}$ means that the hole can hop to the next-
nearest neighbors without distorting the AFM background. This gives rise to a
bare hole dispersion, which starts to dominate the dynamics as the AFM order
decreases with $J_{2}/J_{1}$ eventually leading to an increase in ballistic
expansion of a bare hole. At a fundamental level, the hole interacts with the
AFM background by emitting or absorbing a _single_ spin wave [46], and these
scattering processes become off-resonant when the characteristic energy of the
bare hole dispersion is much larger than the spin wave energies, from which it
follows that hole becomes less dressed.

The situation is qualitatively different in the QSL. Here, a basic interaction
vertex corresponds to the emission or absorption of _two_ spinons as can be
seen in Eq. (10) and Fig. (3). This makes resonant holon-spinon scattering
possible even when the characteristic energy of the spinon spectrum is much
smaller than that of the bare holon. It follows that the holon and therefore
the hole is always heavily dressed by spinons for strong interactions, even
when the characteristic energy of the bare holon dispersion is large, which
slows it down. While we unfortunately cannot calculate the hole velocity for
more values of $J_{2}/J_{1}$ inside the QSL, since we are only aware of
numerical results providing a mapping to the mean-field description for
$J_{2}/J_{1}=1/2$ [43], our results do allow us to predict that the ballistic
expansion velocity of a hole depends _non-monotonically_ on the ratio
$J_{2}/J_{1}$: In the AFM phase it increases whereas it becomes smaller again
in the spin liquid phases as can be seen in Fig. 1. The maximum velocity of
the hole is likely to occur close to the quantum phase transition between the
AFM and the QSL, and since this result is due to the different kinematics in
the hole-medium scattering in the two phases, we believe it is robust and not
an artefact of any approximations. Hence, we conclude that observing the
expansion dynamics of a hole is a promising way to detect the transition
between and AFM and a QSL.

##  VI Conclusions and outlook

Inspired by the remarkable advances with experimental platforms based on
optical lattices and 2D quantum materials, we explored the equilibrium and
non-equilibrium properties of a hole in a QSL. Using a field theory approach,
we identified the presence of quasiparticle string-like states of the hole,
originating from breaking up singlets into their constituent spinons. The hole
was furthermore shown to expand ballistically as a quasiparticle throughout
the QSL for long times after its injection at a given lattice site, with a
velocity exhibiting a characteristic non-monotonic behavior as a quantum phase
transition between the AFM and the QSL phases is crossed. As such, our results
suggests promising new ways to use the hole as a quantum probe for QSLs.

In addition to being relevant for new experiments exploring atoms in optical
lattices and 2D quantum materials, our paper opens up several new research
directions. This includes examining the possible presence of holon-spinon
bound states, which in our field theory approach should show up as poles in
the holon-spinon scattering matrix. It would be interesting to connect such
calculations with earlier results stating that bound spinon-holon states are
caused by an attractive interaction mediated by gauge fluctuations [5].
Another intriguing problem concerns the role of gauge fluctuations on the hole
properties. While such gauge fluctuations likely are not important for the
present $\mathbb{Z}_{2}$ QSL where they are gapped [8] and the low-energy
spectrum of the mean-field description is similar to that after projection
[45], it would be interesting to compare with calculations performed using VMC
[64] and MPS methods [26, 27]. Gauge fluctuations are moreover crucial for
describing other QSLs such as those realised in triangular lattices [8].
Finally, it would be very interesting to explore the binding of two holes in a
QSL, which could lead to new and exotic superconductors/superfluids for a non-
zero hole concentration.

###### Acknowledgements.

This work has been supported by the Danish National Research Foundation
through the Center of Excellence “CCQ” (Grant Agreement No. DNRF156), as well
as by the Carlsberg Foundation through a Carlsberg Internationalisation
Fellowship, Grant No. CF21_0410. J.H.N. would like to thank Leon Balents and
Federico Becca for useful comments and discussions.

##  Appendix A Extended self-consistent Born approximation

In this appendix, we investigate what happens if we include the type 2
diagrams from Fig. 3 when we calculate the Green’s function of the holon, Eq.
(15). Looking at type 1 and 2, we see that type 2 is a higher order diagram,
with the order referring to the number of vertices. When including these
diagrams, we, therefore, ensure that we keep the same order of diagrams. We do
this by first including both types of diagrams in the iterative calculations,
and then after 2/3 of the iterations we only include type 1. By including the
type 2 diagrams, we cannot go much beyond a lattice size of 8x8, but in Fig.
14 we see that including both types compared to only type 1 does not change
the overall behavior of the spectral function. These two results represent the
general tendencies.

Figure 14: Spectral function of the holon calculated using both types of
diagrams from Fig. 3 (full line) and only type 1 (dashed line). These
calculations are performed with $\epsilon/t_{1}=1/20$ for ${\bf k}=(\pi,\pi)$
in panel (a), ${\bf k}=(0,0)$ in panel (b), and a lattice size of 8x8. Except
for small discrepancies at the top of the spectrum, we find that the overall
behavior does not change by including the type 2 diagrams. This demonstrates
the general tendency.

For the calculation in the main text we hence only include the type 1 diagrams
in order to perform the calculations for larger system sizes.

##  Appendix B String excitations

In this Appendix, we further investigate the nature of the string excitations.
We do this by tuning the pairing potentials $\Delta_{1}/\epsilon$ while
keeping $\Delta_{1}/\Delta_{2}$ constant. In doing this the average excitation
energy of the spinons, $\omega^{s}_{avg}$ decreases, and we retrieve a softer
spinon spectrum. This results in a background with longer cooper pairs and
hence longer singlets, which is seen in Fig. 15 by plotting the average length
of the cooper pairs,

$\displaystyle
d=\sqrt{\frac{\sum_{\mathbf{d}}|\Delta_{\mathbf{d}}|^{2}|\mathbf{d}|^{2}}{\sum_{\mathbf{d}}|\Delta_{\mathbf{d}}|^{2}}},\
\Delta_{\mathbf{d}}=\frac{1}{N}\sum_{{\bf r}}\left\langle\hat{f}_{{\bf
r},\uparrow}\hat{f}_{{\bf r}+\mathbf{d},\downarrow}\right\rangle,$ (25)

as a function of $\Delta_{1}/\epsilon$.

Figure 15: Cooper pair length, Eq. (25), For changing pairing potentials. As
$\Delta_{{\bf k}}/t$ is increased, the average length of the cooper pairs is
lowered and flattens out around $\Delta_{1}/\epsilon=1$. In the inset, we find
the strength of string potential plotted as a function of
$\Delta_{1}/\epsilon=1$. This depends linearly on $\Delta_{1}/\epsilon=1$,
indicating that the energy penalty of breaking up the cooper pairs increases
with $\Delta_{1}/\epsilon=1$.

In an environment with more low-lying excitations, we also find a much heavier
quasiparticle, See Fig. 16. With the average cost of exciting the spinons
decreasing, we would then expect the linear potential, given by the slope in
Fig. 6, to decrease. This slope is plotted in the inset in Fig. 16, and we see
just that.

Figure 16: Holon spectral function with $\epsilon/t_{1}=0.5$,
$t_{2}/t_{1}=-\sqrt{0.5}$, and $\Delta_{{\bf k}}\rightarrow 0.1\Delta_{{\bf
k}}$ for a system size of 32x32. Comparing with Fig. 5 we see that the QP
bandwidth has decreased with a factor of $3$ and that fewer crystal momenta
state has a well-defined peak.

Another feature found is that, by decreasing the pairing potential, we need a
larger interaction strength before the string excitations start to appear. In
Fig. 17 we show results for $\Delta_{{\bf k}}\rightarrow 0.1\Delta_{{\bf k}}$
similar to those presented in Fig. 6(b).

Figure 17: Illustrating the importante of the pairing potential, we plot in
panel (a) the poles of the Green’s function as a function of the interaction
strength with $\Delta_{{\bf k}}\rightarrow 0.1\Delta_{{\bf k}}$. The excited
states do not show before $\epsilon/t_{1}\sim 1/16$, which is at a three times
larger interaction strength compared to the results presented in Fig. 6(b).
Panel (b) summarizes the behavior by showing the interaction strength needed
for the string excitations to show when the pairing potentials are changed.
The red line illustrates when then string excitations show, and is the
boundary between the blue region where no string excitations are present and
the beige-colored region with string excitations present.

Here, we see that the first excited states at an interaction strength below
$\epsilon/t_{1}\sim 1/16$, and the ones corresponding to the third pole in
Fig. 6(a) do not show for the $\epsilon/t_{1}$ values shown. This behavior is
summarized in Fig. 17(b), indicating for which parameters string excitations
occur. This, therefore, means that the tightly bound cooper pairs promote the
formation of these states. We believe the reason for this is the much heavier
holon quasiparticle which hinders the formation of the strings.

##  Appendix C Linear spin wave theory

In this appendix, we will use the Linear Spin Wave Theory (LSWT) on the
$t$-$J_{1}$-$J_{2}$ given by Eq. (2) and (3). Considering a single dopant, we
will proceed as done [46, 48, 58, 57, 62, 50, 53, 63], which utilizes a
generalized Holstein-Primarkorff. This relys on the system supporting AFM
order, such that we can define two sublattices where the spins predominantly
point up on sublattice A and down on sublattice B. Using this, we express the
operators on sublattice A as
$\hat{S}^{z}_{l,{\mathbf{i}}}=(1-\hat{h}^{\dagger}_{l,{\mathbf{i}}}\hat{h}_{l,{\mathbf{i}}})/2-\hat{s}^{\dagger}_{l,{\mathbf{i}}}\hat{s}_{l,{\mathbf{i}}}$
and
$\hat{S}_{l,{\mathbf{i}}}^{-}=\hat{S}^{x}_{l,{\mathbf{i}}}-i\hat{S}^{y}_{l,{\mathbf{i}}}=\hat{s}^{\dagger}_{l,{\mathbf{i}}}(1-\hat{s}^{\dagger}_{l,{\mathbf{i}}}\hat{s}_{l,{\mathbf{i}}}-\hat{h}^{\dagger}_{l,{\mathbf{i}}}\hat{h}_{l,{\mathbf{i}}})^{1/2}$,
and the creation operators are $\tilde{c}_{l,{\bf
i},{\downarrow}}=\hat{h}^{\dagger}_{l,{\mathbf{i}}}\hat{s}_{l,{\mathbf{i}}}$
and $\tilde{c}_{l,{\bf
i},{\uparrow}}=\hat{h}^{\dagger}_{l,{\mathbf{i}}}(1-\hat{s}^{\dagger}_{l,{\mathbf{i}}}\hat{s}_{l,{\mathbf{i}}}-\hat{h}^{\dagger}_{l,{\mathbf{i}}}\hat{h}_{l,{\mathbf{i}}})^{1/2}$.
Opposite to the slave-boson description used in the main text,
$\hat{h}^{\dagger}_{l,{\mathbf{i}}}$ is a fermionic creation operator of a
hole, and $\hat{s}^{\dagger}_{l,{\mathbf{i}}}$ is a bosonic creation operator
describing derivation from the AFM ordered state. Similar relations are used
for sublattice B.

Using these expressions and keeping only linear terms, we can diagonalize
$J$-part of the Hamiltonian to find

$\displaystyle\hat{H}_{J}=\sum_{{\bf k}}\omega^{m}_{{\bf
k}}\hat{b}^{\dagger}_{{\bf k}}\hat{b}_{{\bf k}}+E_{0},$ (26)

with

$\displaystyle\omega^{m}_{{\bf k}}$ $\displaystyle=\
\frac{1}{2}\sqrt{\tilde{J}^{2}-(\alpha_{1}J_{1}z\gamma_{{\bf k}})^{2}},$
$\displaystyle\tilde{J}_{{\bf k}}$ $\displaystyle\equiv
z\left(J_{1}+J_{2}\left[\alpha_{2}\gamma^{\prime}_{{\bf k}}-1\right]\right),$
(27)

for $z=4$ ands $\gamma_{{\bf k}}=\sum_{\boldsymbol{\delta}}e^{i{\bf
k}\cdot\boldsymbol{\delta}}/z$ and $\gamma^{\prime}_{{\bf
k}}=\sum_{\boldsymbol{\delta}^{\prime}}e^{i{\bf
k}\cdot\boldsymbol{\delta}^{\prime}}/z$ being the structure factor for NN and
NNN terms respectively. The spin wave operators, $\hat{b}$, is related the
spin fluctuations, $\hat{s}$, by a Bogoliubov transformation

$\displaystyle\begin{bmatrix}\hat{s}_{{\bf k}}\\\ \hat{s}^{\dagger}_{-{\bf
k}}\end{bmatrix}=\begin{bmatrix}u_{{\bf k}}&-v_{{\bf k}}\\\ -v_{{\bf
k}}&u_{{\bf k}}\end{bmatrix}\begin{bmatrix}\hat{b}_{{\bf k}}\\\
\hat{b}^{\dagger}_{-{\bf k}}\end{bmatrix},$ (28)

where $u_{{\bf k}}=\sqrt{\frac{1}{2}\left(\tilde{J}/(2\omega^{m}_{{\bf
k}})+1\right)}$ and $v_{{\bf k}}={\rm sgn}\left[\gamma_{\bf
k}\right]\sqrt{\frac{1}{2}\left(\tilde{J}/(2\omega^{m}_{{\bf k}})-1\right)}$
are the coherence factors. Looking at Fig. 18, we see that this approach
qualitative captures the behavior when $J_{2}$ is increased. This predicts a
phase transition to occur at $J_{2}/J_{1}=0.38$, which is not far from the
predictions made with more advanced methods. This yields it to be around
$J_{2}/J_{1}\sim 0.46$ [34, 35, 36, 37].

Figure 18: The average magnetization, $2|\left\langle
S_{z}\right\rangle|=1-2\sum_{{\bf k}}v^{2}_{{\bf k}}/N$, calculated within
LSWT. This predicts a phase transition to occur at $J_{2}/J_{1}=0.38$.

With the described transformations, the kinetic term takes the form

$\displaystyle\hat{H}_{t}=$ $\displaystyle\sum_{\mathbf{k}}\omega^{p}_{{\bf
k}}\hat{h}_{\mathbf{k}}^{\dagger}\hat{h}_{\mathbf{k}}+\sum_{{\bf k},{\bf
p}}g({\bf p},{\bf k})\hat{h}_{{\bf p}+{\bf k}}^{\dagger}\hat{h}_{{\bf
p}}\left(\hat{b}^{\dagger}_{-{\bf k}}+\hat{b}_{{\bf k}}\right)$ (29)

with

$\displaystyle g({\bf p},{\bf k})=\frac{zt_{1}}{\sqrt{N}}\left[u_{{\bf
k}}\gamma_{{\bf p}+{\bf k}}-v_{{\bf k}}\gamma_{\bf p}\right],$ (30)

and

$\displaystyle\omega^{p}_{{\bf
k}}=2t_{2}\left[\cos{(k_{x}+k_{y})}+\cos{(k_{x}-k_{y})}\right].$ (31)

For a finite $J_{2}$ we hence get a bare dispersion corresponding to the hole
jumping within the same sublattice. Proceeding as done in the literature [46,
48, 58, 57, 62, 50, 53, 63], we can use the above expressions to calculate the
Green’s function,

$\displaystyle G({\bf k},\tau)=-i\left\langle T_{\tau}\left[\hat{h}_{{\bf
k}}(\tau)\hat{h}^{\dagger}_{{\bf k}}(0)\right]\right\rangle,$ (32)

of a single hole with the SCBA approximation.

Investigation the nature of the rms speed-up as we approach the phase
transition, we find that by neglecting the bare dispersion of the hole we
retrieve an rms-velocity which decreases with interaction strength, see Fig.
19. This we ascribe to the fact that, with an increasing $J_{2}/J_{1}$ the
spinon spectrum gets increasingly softer around ${\bf p}=(\pi,0)$ and ${\bf
p}=(0,\pi)$.

Figure 19: Same plot as in Fig. 1, but where we have omitted the bare
dispersion of the hole in the AFM phase. As we approach the phase transition
from the AFM side, the hole experiences a slow-down instead of a speed-up as
in Fig. 1.

##  Appendix D String excitations for hole

In this appendix, we show that indeed the higher excitation peaks in the
hole’s spectral function have the characteristic two-third-scaling associated
with string excitations. This is seen in Fig. 20, where we plot the energy of
the peak as a function of $(\epsilon/t_{1})^{2/3}$. Here, we see the linear
behavior and in addition, we see that the energy is given by the energy of the
holon string peak plus $\omega_{(\pi,\pi)}$ corresponding to the Van Hove
singularity discussed in the main part. This shows that the excitation peaks
in the hole spectral function come from the string excitations of the holon.

Figure 20: Excitation energy of the excitation peak from the hole spectral
function for ${\bf p}=(0,0)$ in orange. Equivalent to Eq. (23), the blue data
points show the energy of the holon excitation peak for ${\bf p}=(\pi,\pi)$
added to $\omega_{(\pi,\pi)}$ coming from the Van Hove singularity discussed
in the main part. This is plotted as a function of $(\epsilon/t_{1})^{2/3}$.

##  Appendix E Quasiparticle band fit

Interestingly, we find that the lowest quasiparticle band can be fitted to the
form

$\displaystyle\omega^{\rm QP}_{{\bf k}}=$
$\displaystyle\tilde{t}_{1}(\cos{k_{x}}+\cos{k_{y}})+\tilde{t}_{2}[\cos{(k_{x}+k_{y})}+\cos{(k_{x}-k_{y})}]$
$\displaystyle+\tilde{t}_{3}(\cos{2k_{x}}+\cos{2k_{y}}),$ (33)

shown as the red dashed line in Fig. 5(b), where
$\tilde{t}_{2}/\tilde{t}_{1}\approx 0.15$ and
$\tilde{t}_{3}/\tilde{t}_{1}\approx-0.05$. We find that these ratios are
approximately constant for strong interaction strengths $\epsilon/t_{1}<0.1$.

## References

  * Anderson [1973] P. Anderson, Resonating valence bonds: A new kind of insulator?, Materials Research Bulletin 8, 153 (1973).
  * Savary and Balents [2017] L. Savary and L. Balents, Quantum spin liquids: A review, Reports on Progress in Physics 80, 016502 (2017).
  * Broholm  _et al._ [2020] C. Broholm, R. J. Cava, S. A. Kivelson, D. G. Nocera, M. R. Norman, and T. Senthil, Quantum spin liquids, Science 367, eaay0668 (2020).
  * Kalmeyer and Laughlin [1987] V. Kalmeyer and R. B. Laughlin, Equivalence of the resonating-valence-bond and fractional quantum Hall states, Physical Review Letters 59, 2095 (1987).
  * Lee  _et al._ [2006] P. A. Lee, N. Nagaosa, and X.-G. Wen, Doping a Mott insulator: Physics of high-temperature superconductivity, Reviews of Modern Physics 78, 17 (2006).
  * Knolle and Moessner [2019] J. Knolle and R. Moessner, A Field Guide to Spin Liquids, Annual Review of Condensed Matter Physics 10, 451 (2019).
  * Baskaran  _et al._ [1987] G. Baskaran, Z. Zou, and P. Anderson, The resonating valence bond state and high-Tc superconductivity — A mean field theory, Solid State Communications 63, 973 (1987).
  * Wen [2002] X.-G. Wen, Quantum orders and symmetric spin liquids, Physical Review B 65, 165113 (2002).
  * Poilblanc  _et al._ [2012] D. Poilblanc, N. Schuch, D. Pérez-García, and J. I. Cirac, Topological and entanglement properties of resonating valence bond wave functions, Physical Review B 86, 014404 (2012).
  * Hu  _et al._ [2013] W.-J. Hu, F. Becca, A. Parola, and S. Sorella, Direct evidence for a gapless Z 2 spin liquid by frustrating Néel antiferromagnetism, Physical Review B 88, 060402 (2013).
  * Meng  _et al._ [2010] Z. Y. Meng, T. C. Lang, S. Wessel, F. F. Assaad, and A. Muramatsu, Quantum spin liquid emerging in two-dimensional correlated Dirac fermions, Nature 464, 847 (2010).
  * Yan  _et al._ [2011] S. Yan, D. A. Huse, and S. R. White, Spin-Liquid Ground State of the _S_ = 1/2 Kagome Heisenberg Antiferromagnet, Science 332, 1173 (2011).
  * 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, A. Omran, S. Sachdev, A. Vishwanath, M. Greiner, V. Vuletić, and M. D. Lukin, Probing topological spin liquids on a programmable quantum simulator, Science 374, 1242 (2021).
  * Yamamoto  _et al._ [2020] R. Yamamoto, H. Ozawa, D. C. Nak, I. Nakamura, and T. Fukuhara, Single-site-resolved imaging of ultracold atoms in a triangular optical lattice, New Journal of Physics 22, 123028 (2020).
  * Yang  _et al._ [2021] J. Yang, L. Liu, J. Mongkolkiattichai, and P. Schauss, Site-resolved imaging of ultracold fermions in a triangular-lattice quantum gas microscope, PRX Quantum 2, 020344 (2021).
  * Sun  _et al._ [2023] B.-Y. Sun, N. Goldman, M. Aidelsburger, and M. Bukov, Engineering and Probing Non-Abelian Chiral Spin Liquids Using Periodically Driven Ultracold Atoms, PRX Quantum 4, 020329 (2023).
  * Lebrat  _et al._ [2024] M. Lebrat, M. Xu, L. H. Kendrick, A. Kale, Y. Gang, P. Seetharaman, I. Morera, E. Khatami, E. Demler, and M. Greiner, Observation of nagaoka polarons in a fermi–hubbard quantum simulator, Nature 629, 317 (2024).
  * Prichard  _et al._ [2024] M. L. Prichard, B. M. Spar, I. Morera, E. Demler, Z. Z. Yan, and W. S. Bakr, Directly imaging spin polarons in a kinetically frustrated hubbard system, Nature 629, 323 (2024).
  * Pan  _et al._ [2020] H. Pan, F. Wu, and S. Das Sarma, Quantum phase diagram of a moiré-hubbard model, Phys. Rev. B 102, 201104 (2020).
  * Zang  _et al._ [2021] J. Zang, J. Wang, J. Cano, and A. J. Millis, Hartree-fock study of the moiré hubbard model for twisted bilayer transition metal dichalcogenides, Phys. Rev. B 104, 075150 (2021).
  * Morales-Durán  _et al._ [2022] N. Morales-Durán, N. C. Hu, P. Potasz, and A. H. MacDonald, Nonlocal interactions in moiré hubbard systems, Phys. Rev. Lett. 128, 217202 (2022).
  * Kitaev [2006] A. Kitaev, Anyons in an exactly solved model and beyond, Annals of Physics 321, 2 (2006).
  * Jiang  _et al._ [2012] H.-C. Jiang, H. Yao, and L. Balents, Spin liquid ground state of the spin- 1 2 square J 1 - J 2 Heisenberg model, Physical Review B 86, 024424 (2012).
  * Iqbal  _et al._ [2016] Y. Iqbal, W.-J. Hu, R. Thomale, D. Poilblanc, and F. Becca, Spin liquid nature in the Heisenberg J 1 - J 2 triangular antiferromagnet, Physical Review B 93, 144411 (2016).
  * Zhu and White [2015] Z. Zhu and S. R. White, Spin liquid phase of the S = 1 2 J 1 - J 2 Heisenberg model on the triangular lattice, Physical Review B 92, 041105 (2015).
  * Kadow  _et al._ [2022] W. Kadow, L. Vanderstraeten, and M. Knap, Hole spectral function of a chiral spin liquid in the triangular lattice hubbard model, Phys. Rev. B 106, 094417 (2022).
  * kad [2024] Single-hole spectra of kitaev spin liquids: from dynamical nagaoka ferromagnetism to spin-hole fractionalization, npj Quantum Materials 9, 32 (2024).
  * Koepsell  _et al._ [2019] J. Koepsell, J. Vijayan, P. Sompet, F. Grusdt, T. A. Hilker, E. Demler, G. Salomon, I. Bloch, and C. Gross, Imaging magnetic polarons in the doped Fermi–Hubbard model, Nature 572, 358 (2019).
  * Ji  _et al._ [2021] G. Ji, M. Xu, L. H. Kendrick, C. S. Chiu, J. C. Brüggenjürgen, D. Greif, A. Bohrdt, F. Grusdt, E. Demler, M. Lebrat, and M. Greiner, Coupling a Mobile Hole to an Antiferromagnetic Spin Background: Transient Dynamics of a Magnetic Polaron, Physical Review X 11, 021022 (2021).
  * Laughlin [1995] R. B. Laughlin, Numerical evidence for electron decay in the two-dimensional t-J model, Journal of Low Temperature Physics 99, 443 (1995).
  * Béran  _et al._ [1996] P. Béran, D. Poilblanc, and R. Laughlin, Evidence for composite nature of quasiparticles in the 2D t-J model, Nuclear Physics B 473, 707 (1996).
  * Tang  _et al._ [2013] E. Tang, M. P. A. Fisher, and P. A. Lee, Low-energy behavior of spin-liquid electron spectral functions, Physical Review B 87, 045119 (2013).
  * Rantner and Wen [2001] W. Rantner and X.-G. Wen, Electron Spectral Function and Algebraic Spin Liquid for the Normal State of Underdoped High T c Superconductors, Physical Review Letters 86, 3871 (2001).
  * Liu  _et al._ [2022] W.-Y. Liu, S.-S. Gong, Y.-B. Li, D. Poilblanc, W.-Q. Chen, and Z.-C. Gu, Gapless quantum spin liquid and global phase diagram of the spin-1/2 J 1 - J 2 square antiferromagnetic Heisenberg model, Science Bulletin 67, 1034 (2022).
  * Wang and Sandvik [2018] L. Wang and A. W. Sandvik, Critical Level Crossings and Gapless Spin Liquid in the Square-Lattice Spin- 1 / 2 J 1 - J 2 Heisenberg Antiferromagnet, Physical Review Letters 121, 107202 (2018).
  * Nomura and Imada [2021] Y. Nomura and M. Imada, Dirac-Type Nodal Spin Liquid Revealed by Refined Quantum Many-Body Solver Using Neural-Network Wave Function, Correlation Ratio, and Level Spectroscopy, Physical Review X 11, 031034 (2021).
  * Ferrari and Becca [2020a] F. Ferrari and F. Becca, Gapless spin liquid and valence-bond solid in the J 1 - J 2 Heisenberg model on the square lattice: Insights from singlet and triplet excitations, Physical Review B 102, 014417 (2020a).
  * Chalopin  _et al._ [2024] T. Chalopin, P. Bojović, D. Bourgund, S. Wang, T. Franz, I. Bloch, and T. Hilker, Optical superlattice for engineering Hubbard couplings in quantum simulation (2024).
  * Song  _et al._ [2021] X.-Y. Song, A. Vishwanath, and Y.-H. Zhang, Doping the chiral spin liquid: Topological superconductor or chiral metal, Physical Review B 103, 165138 (2021).
  * Brunkert and Punk [2020] J. Brunkert and M. Punk, Slave-boson description of pseudogap metals in t - J models, Physical Review Research 2, 043019 (2020).
  * Ferrari and Becca [2018] F. Ferrari and F. Becca, Spectral signatures of fractionalization in the frustrated Heisenberg model on the square lattice, Physical Review B 98, 100405 (2018).
  * Ferrari and Becca [2020b] F. Ferrari and F. Becca, Gapless spin liquid and valence-bond solid in the J 1 - J 2 Heisenberg model on the square lattice: Insights from singlet and triplet excitations, Physical Review B 102, 014417 (2020b).
  * Yu  _et al._ [2018] S.-L. Yu, W. Wang, Z.-Y. Dong, Z.-J. Yao, and J.-X. Li, Deconfinement of spinons in frustrated spin systems: Spectral perspective, Physical Review B 98, 134410 (2018).
  * Wen [2010] X.-G. Wen,  _Quantum Field Theory of Many-Body Systems: From the Origin of Sound to an Origin of Light and Electrons_ , repr ed., Oxford Graduate Texts (Oxford University Press, Oxford, 2010).
  * Ferrari and Becca [2019] F. Ferrari and F. Becca, Dynamical Structure Factor of the J 1 - J 2 Heisenberg Model on the Triangular Lattice: Magnons, Spinons, and Gauge Fields, Physical Review X 9, 031026 (2019).
  * Kane  _et al._ [1989] C. L. Kane, P. A. Lee, and N. Read, Motion of a single hole in a quantum antiferromagnet, Physical Review B 39, 6880 (1989).
  * Jiang  _et al._ [2021] S. Jiang, D. J. Scalapino, and S. R. White, Ground-state phase diagram of the _t-t_ $\prime$ _-__J_ model, Proceedings of the National Academy of Sciences 118, e2109978118 (2021).
  * Martinez and Horsch [1991] G. Martinez and P. Horsch, Spin polarons in the _t_ \- _J_ model, Physical Review B 44, 317 (1991).
  * Diamantis and Manousakis [2021] N. G. Diamantis and E. Manousakis, Dynamics of string-like states of a hole in a quantum antiferromagnet: A diagrammatic Monte Carlo simulation, New Journal of Physics 23, 123005 (2021).
  * Nielsen  _et al._ [2022] K. K. Nielsen, T. Pohl, and G. M. Bruun, Nonequilibrium Hole Dynamics in Antiferromagnets: Damped Strings and Polarons, Physical Review Letters 129, 246601 (2022).
  * Laughlin [1997] R. B. Laughlin, Evidence for Quasiparticle Decay in Photoemission from Underdoped Cuprates, Physical Review Letters 79, 1726 (1997).
  * Brinkman and Rice [1970] W. F. Brinkman and T. M. Rice, Single-Particle Excitations in Magnetic Insulators, Physical Review B 2, 1324 (1970).
  * Nyhegn  _et al._ [2022] J. H. Nyhegn, K. K. Nielsen, and G. M. Bruun, Equilibrium and nonequilibrium dynamics of a hole in a bilayer antiferromagnet, Physical Review B 106, 155160 (2022).
  * Skou  _et al._ [2021] M. G. Skou, T. G. Skov, N. B. Jørgensen, K. K. Nielsen, A. Camacho-Guardian, T. Pohl, G. M. Bruun, and J. J. Arlt, Non-equilibrium quantum dynamics and formation of the Bose polaron, Nature Physics 17, 731 (2021).
  * Nielsen  _et al._ [2019] K. K. Nielsen, L. A. P. Ardila, G. M. Bruun, and T. Pohl, Critical slowdown of non-equilibrium polaron dynamics, New Journal of Physics 21, 043014 (2019).
  * Bohrdt  _et al._ [2020] A. Bohrdt, F. Grusdt, and M. Knap, Dynamical formation of a magnetic polaron in a two-dimensional quantum antiferromagnet, New Journal of Physics 22, 123023 (2020).
  * Schmitt-Rink  _et al._ [1988] S. Schmitt-Rink, C. M. Varma, and A. E. Ruckenstein, Spectral Function of Holes in a Quantum Antiferromagnet, Physical Review Letters 60, 2793 (1988).
  * Liu and Manousakis [1991] Z. Liu and E. Manousakis, Spectral function of a hole in the _t_ \- _J_ model, Physical Review B 44, 2414 (1991).
  * Marsiglio  _et al._ [1991] F. Marsiglio, A. E. Ruckenstein, S. Schmitt-Rink, and C. M. Varma, Spectral function of a single hole in a two-dimensional quantum antiferromagnet, Physical Review B 43, 10882 (1991).
  * Liu and Manousakis [1992] Z. Liu and E. Manousakis, Dynamical properties of a hole in a Heisenberg antiferromagnet, Physical Review B 45, 2425 (1992).
  * Chernyshev and Leung [1999] A. L. Chernyshev and P. W. Leung, Holes in the t - J z model: A diagrammatic study, Physical Review B 60, 1592 (1999).
  * Nielsen  _et al._ [2021] K. K. Nielsen, M. A. Bastarrachea-Magnani, T. Pohl, and G. M. Bruun, Spatial structure of magnetic polarons in strongly interacting antiferromagnets, Physical Review B 104, 155136 (2021).
  * Nyhegn  _et al._ [2023] J. H. Nyhegn, G. M. Bruun, and K. K. Nielsen, Wave function and spatial structure of polarons in an antiferromagnetic bilayer 10.48550/ARXIV.2305.05980 (2023).
  * Charlebois and Imada [2020] M. Charlebois and M. Imada, Single-Particle Spectral Function Formulated and Calculated by Variational Monte Carlo Method with Application to d -Wave Superconducting State, Physical Review X 10, 041023 (2020).