11institutetext: DOT-Lab Research Center, University of Quebec (TELUQ),
Montreal, Canada
11email<EMAIL_ADDRESS>
# Unicode at Gigabytes per Second

Daniel Lemire

###### Abstract

We often represent text using Unicode formats (UTF-8 and UTF-16). The UTF-8
format is increasingly popular, especially on the web (XML, HTML, JSON, Rust,
Go, Swift, Ruby). The UTF-16 format is most common in Java, .NET, and inside
operating systems such as Windows.

Software systems frequently have to convert text from one Unicode format to
the other. While recent disks have bandwidths of
$5\text{\,}\mathrm{GiB}\text{/}\mathrm{s}$ or more, conventional approaches
transcode non-ASCII text at a fraction of a gigabyte per second.

We show that we can validate and transcode Unicode text at gigabytes per
second on current systems (x64 and ARM) without sacrificing safety. Our open-
source library can be ten times faster than the popular ICU library on non-
ASCII strings and even faster on ASCII strings.

###### Keywords:

Unicode Vectorization Internationalization.

##  1 Introduction

††footnotetext: This manuscript is based on a forthcoming long-form article
written with Wojciech Muła, _Transcoding Billions of Unicode Characters per
Second with SIMD Instructions_.

From the early days of computing, programmers have had to represent characters
in software. They needed to agree to standards so that different software
written by different vendors would be interoperable. One of the earliest such
standards is ASCII—first specified in the early 1960s. The ASCII standard is
still popular today: it uses one byte per character—with the most significant
bit set to zero. Unfortunately, ASCII could only ever represent up to 128
characters—far less than needed.

Thus many diverging standards emerged for representing characters in software.
The existence of multiple incompatible formats made the production of
interoperable localized software difficult. Conversion between some formats
could sometimes be lossy or ambiguous.

Unicode arose in the late 1980s as an attempt to provide a single agreed-upon
standard. Initially, it was believed that using 16 bits per character would be
sufficient. However, engineers realized over time that a wider range of
characters should be supported—if the standard was to be universal. Thus the
Unicode standard was extended to potentially include up to $1\,114\,112$
characters. Characters are sometimes called code points and represented as
integer values between 0 and $1\,114\,112$. In practice, only a small fraction
of all possible code points have been assigned, but more are assigned over
time with each Unicode revision. The Unicode standard is an extension of the
ASCII standard in the sense that the first 128 Unicode code-point values match
the ASCII characters.

There are several ways to represent Unicode characters in bytes. Due to the
original expectation that Unicode would fit in 16-bit space, a format based on
16-bit words (UTF-16) format was published in 1996 and formalized in 2000 [3].
It may use either 16-bit or 32-bit per character. The UTF-16 format was
adopted by programming languages such as Java, and became a default under
Windows.

Unfortunately, UTF-16 is not backward compatible with ASCII at a byte level.
Thus an ASCII compatible format was proposed and formalized in 2003: UTF-8
[8]. Over time, it became widely used. Text interchange formats such as JSON,
HTML or XML are expected to be in the UTF-8 format. Programming languages such
as Go, Rust and Swift use UTF-8 by default. When used as part of data
interchange documents, UTF-8 is commonly more concise due to its ability to
use one byte per character to represent ASCII text.

Though UTF-8 dominates in many applications, it does not make UTF-16 obsolete.
The UTF-16 format has advantages. Indeed, most text represented in the UTF-16
format has exactly 2-byte per character, except for the occasional special
character (e.g., an emoji). Having a flat 2-byte per character makes some
operations faster. Both formats require validation: not all arrays of bytes
are valid. However, the UTF-8 format is more expensive to validate. In some
cases, UTF-16 may be even more concise as well (e.g., when representing
Chinese text).

Thus, for the foreseeable future, we need to validate both UTF-16 and UTF-8
strings, and to convert (transcode) text between the two formats.

We should expect these operations to be fast, and they are. However, speed and
efficiency are relative. Cloud vendors offer high bandwidth between node
instances (e.g., $3.3\text{\,}\mathrm{GiB}\text{/}\mathrm{s}$) and the
bandwidth of disks is rising fast (e.g.,
$5\text{\,}\mathrm{GiB}\text{/}\mathrm{s}$ with PCIe 4.0) [1]. We believe that
we should be able to match such high speeds inside the processor when
transcoding text.

For fast processing, we should seek to make the best possible use out of our
processors. Commodity processors support single-instruction-multiple-data
(SIMD) instructions. These instructions operate on several words at once
unlike regular instructions. Starting with the Haswell microarchitecture
(2013), Intel and AMD processors support the AVX2 instruction set and 256-bit
vector registers. Most mobile phones tablets have 64-bit ARM processors
(aarch64) with NEON instructions (128-bit registers). Hence, on recent x64
processors, we can compare two strings of 32 characters in a single
instruction. Algorithms designed for SIMD instructions typically require fewer
instructions per byte. Software that retires fewer instructions tends to use
less power and to be faster.

##  2 Related Work

There are many ways to validate and transcode Unicode text. One may use series
of branches or finite-state approaches [4]. We are most interested in the
fastest techniques. Keiser and Lemire [6] describe a fast SIMD-based UTF-8
validation directly on byte streams. We make use of their approach (see § 4).

Cameron [2] proposed that we transform text using _bit streams_. Given byte
arrays, the bit-stream approach creates eight bit arrays, each one
corresponding to a bit position within a byte. There is one bit stream for the
least significant bit, one for the second least significant bit, and so forth.
We must first transform the input data into such bit streams and then convert
the data back from bit streams to an array of bytes. Cameron applied this
strategy to UTF-8 to UTF-16 transcoding.

Independently, Inoue et al. [5] proposed a UTF-8 to UTF-16 SIMD-accelerated
transcoder, but it does not validate its inputs. The authors did not make
their implementation (for PowerPC processors) available.

We are not aware of any further work in the scientific literature regarding
the application of SIMD instruction to the validation or transcoding of
Unicode text. Except for fast ASCII paths, we do not know of any widespread
use of SIMD instructions for validating or transcoding Unicode text.

##  3 The Unicode Formats

ASCII characters require one byte with UTF-8 and two bytes with UTF-16. UTF-16
can represent all characters—except for the supplemental characters such as
emojis—using two bytes. The UTF-8 format uses two bytes for Latin, Hebrew and
Arabic alphabets. Asiatic characters (including Chinese and Japanese) require
three UTF-8 bytes. Both UTF-8 and UTF-16 require 4 bytes for the supplemental
characters. We often represent Unicode characters using its integer value in
hexadecimal as, for example, U+7F (for 127).

UTF-8 encodes values in sequences of one to four bytes. We refer to the first
byte of a sequence as a leading byte; the most significant bits of the leading
byte indicate the length of the sequence:

  * •

If the most significant bit is zero, we have a sequence of one byte (ASCII).

  * •

If the three most significant bits are 110, we have a two-byte sequence.

  * •

If the four most significant bits are 1110, we have a three-byte sequence.

  * •

Finally, if the five most significant bits are 11110, we have a four-byte
sequence.

All bytes following the leading byte in a sequence are continuation bytes, and
they must have their two most significant bits as 10. Except for the required
most significant bit sequences, other bits (from 7 bits to 21 bits) provide
the code-point value. The most significant bits of the code-point value are in
the leading byte, followed by lesser significant bits in the second byte and
so forth, with the least significant bits in the last byte of the sequence.
Valid UTF-8 sequences must follow the following exhaustive rules:

  1. 1.

The five most significant bits of any byte cannot be all ones.

  2. 2.

The leading byte must be followed by the right number of continuation bytes.

  3. 3.

A continuation byte must be preceded by a leading byte.

  4. 4.

The decoded character must be larger than U+7F for two-byte sequences, larger
than U+7FF for three-byte sequences, and larger than U+FFFF for four-byte
sequences.

  5. 5.

The decoded code-point value must be less than $1\,114\,112$.

  6. 6.

The code-point value must not be in the range U+D800…DFFF.

In the UTF-16 format, characters in U+0000…D7FF and U+E000…FFFF are stored as
16-bit values—using two bytes. The characters in the range U+010000…10FFFF
require two 16-bit words (a surrogate pair). The first word in the pair is in
0xD800…DBFF whereas the second word is in 0xDC00…DFFF. The code-point value is
made of the 10 least significant bits of the two words—using the second word
as least significant—adding 0x10000 to the result. During validation, only the
possible surrogate pairs require attention.

##  4 Algorithms

All commodity software with SIMD instructions (e.g., x64, ARM, POWER) have
fast instructions to permute bytes within a SIMD register according to a
sequence of indexes. Our transcoding techniques depend critically on this
feature: we code in a table the necessary parameters—including the indexes
(sometimes called shuffle masks)—necessary to process a variety of incoming
characters.

Our accelerated UTF-8 to UTF-16 transcoding algorithm processes up to 12 input
UTF-8 bytes at a time. From the input bytes, we can quickly determine the
leading bytes and thus the end of each character. We use a 12-bit word as a
key in a 1024-entry table. Each entry in the table contains the number of
UTF-8 bytes that will be consumed and an index into another table where we
find shuffle masks. The value of the index into the other table also
determines one of three possible code paths. The first 64 index values
indicate that we have 6 characters spanning between one and two bytes. Index
values in $[64,145)$ indicate that we have 4 characters spanning between one
and three bytes. The remaining indexes represent the general case: 3
characters spanning between one and four bytes. The shuffle mask can then be
applied to the 12 input bytes to form a vector register that can be
transformed efficiently. We use this 12-byte routine inside 64-byte blocks.
After loading a 64-byte block, we apply the Keiser-Lemire validation routine
[6]. Afterward, we identify the leading bytes, and then process the block in
multiple iterations, using up to 12 bytes each time. In the special case where
all 64 bytes are ASCII, we use a fast path. For even greater efficiency, we
have three other fast paths within the 12-byte routine: we check whether the
next 16 bytes are ASCII bytes, whether they are all two-byte characters, or
all three-byte characters.

Our UTF-16 to UTF-8 algorithm iteratively reads a block of input bytes in a
SIMD register. If all 16-bit words in the loaded SIMD register are in the
range U+0000…007F, we use a fast routine to convert the 16 input bytes into
eight equivalent ASCII bytes. If all 16-bit words are in the range
U+0000…07FF, then we use a fast routine to produce sequences of one-byte or
two-byte UTF-8 characters. Given an 8-bit bitset which indicates which 16-bit
words are ASCII, we load a byte value from a table indicating how many bytes
will be written, and a 16-byte shuffle mask. If all 16-bit words are in the
ranges U+0000…D7FF, U+E000…FFFF, we use another similar specialized routine to
produce sequences of one-byte, two-byte and three-byte UTF-8 characters.
Otherwise, when we detect that the input register contains at least one part
of a surrogate pair, we fall back to a conventional code path.

##  5 Experiments

We make available our software as a portable open-source C++
library.††https://github.com/simdutf/simdutf As a benchmarking system, we use
a recent AMD processor (AMD EPYC 7262, Zen 2 microarchitecture,
$3.39\text{\,}\mathrm{GHz}\text{/}$) and GCC 10. We compare against a popular
library: International Components for Unicode (ICU) [icu] (version 67.1). We
also use the u8u16 library [2]. Unlike ICU and our own work, the u8u16 library
only provides UTF-8 to UTF-16 transcoding. For our experiments, we use lipsum
text in various languages.††https://github.com/rusticstuff/simdutf8 All of our
transcoding tests include validation. To measure the speed, we record the time
by repeating the task $2000$ times. We compare the average time with the
minimal time and find that we have an accuracy of at least 1%. We divide the
input volume by the time required for the transcoding. Fig. 1 shows our
results. Our UTF-8 to UTF-16 transcoding speed exceeds
$4\text{\,}\mathrm{GiB}\text{/}\mathrm{s}$ for Chinese and Japanese texts,
which is about four times faster than ICU. In our tests, the u8u16 library
only surpasses ICU significantly for Arabic. Our UTF-16 to UTF-8 transcoding
speed is nearly $6\text{\,}\mathrm{GiB}\text{/}\mathrm{s}$ in all tests which
is nearly ten times faster than ICU.

For ASCII transcoding (not shown in the figures), we achieve
$36\text{\,}\mathrm{GiB}\text{/}\mathrm{s}$ for UTF-16 to UTF-8 transcoding,
and $20\text{\,}\mathrm{GiB}\text{/}\mathrm{s}$ for UTF-8 to UTF-16
transcoding. Effectively, we are so fast that we are nearly limited by memory
bandwidth. Comparatively, ICU delivers
$2\text{\,}\mathrm{GiB}\text{/}\mathrm{s}$ and
$1\text{\,}\mathrm{GiB}\text{/}\mathrm{s}$ in our tests.

(a) UTF-8 to UTF-16 transcoding (b) UTF-16 to UTF-8 transcoding

Figure 1:  Transcoding speeds for various test files.

##  6 Conclusion

Our SIMD-based transcoders can surpass popular transcoders (e.g., ICU) by a
wide margin (e.g., $4\times$). Our UTF-16 to UTF-8 transcoder achieves speed
of about $6\text{\,}\mathrm{GiB}\text{/}\mathrm{s}$ for many Asiatic languages
using a recent x64 processor. In some cases, we achieve
$4\text{\,}\mathrm{GiB}\text{/}\mathrm{s}$ for UTF-8 to UTF-16 transcoding
with full validation. For ASCII inputs, we achieve tens of gigabytes per
second.

## References

  * [1] Barr, J.: The Floodgates Are Open – Increased Network Bandwidth for EC2 Instances. https://aws.amazon.com/blogs/aws/the-floodgates-are-open-increased-network-bandwidth-for-ec2-instances/ [last checked July 2021] (2018)
  * [2] Cameron, R.D.: A case study in SIMD text processing with parallel bit streams: UTF-8 to UTF-16 transcoding. In: Proceedings of the 13th ACM SIGPLAN Symposium on Principles and practice of parallel programming. pp. 91–98. ACM (2008)
  * [3] Hoffman, P., Yergeau, F.: UTF-16, an encoding of ISO 10646. https://tools.ietf.org/html/rfc2781 [last checked July 2021] (2000), Internet Engineering Task Force, Request for Comments: 3629
  * [4] Höhrmann, B.: Flexible and Economical UTF-8 Decoder. http://bjoern.hoehrmann.de/utf-8/decoder/dfa/ [last checked July 2021] (2010)
  * [5] Inoue, H., Komatsu, H., Nakatani, T.: Accelerating UTF-8 Decoding Using SIMD Instructions (in Japanese). Information Processing Society of Japan Transactions on Programming 1(2), 1–8 (2008)
  * [6] Keiser, J., Lemire, D.: Validating utf-8 in less than one instruction per byte. Software: Practice and Experience 51(5) (2021)
  * [7] International Components for Unicode (UCI). http://site.icu-project.org [last checked July 2021] (2010)
  * [8] Yergeau, F.: UTF-8, a transformation format of ISO 10646. https://tools.ietf.org/html/rfc3629 [last checked July 2021] (2003), Internet Engineering Task Force, Request for Comments: 3629


* [23] Yaozhong Hu, David Nualart, and Jian Song.  Feynman-Kac formula for heat equation driven by fractional white noise.  The Annals of Probability, 39(1):291–326, 2011.
  * [24] David A Huse and Christopher L Henley.  Pinning and roughening of domain walls in Ising systems due to random impurities.  Physical review letters, 54(25):2708, 1985.
  * [25] I. A. Ibragimov and Yu. V. Linnik.  Independent and stationary sequences of random variables.  Wolters-Noordhoff Publishing, Groningen, 1971.  With a supplementary chapter by I. A. Ibragimov and V. V. Petrov, Translation from the Russian edited by J. F. C. Kingman.
  * [26] Svante Janson et al.  Gaussian Hilbert spaces.  Number 129. Cambridge university press, 1997.
  * [27] Maria Jolis.  On a multiple Stratonovich-type integral for some Gaussian processes.  Journal of Theoretical Probability, 19(1):121–133, 2006.
  * [28] Mathew Joseph, Davar Khoshnevisan, and Carl Mueller.  Strong invariance and noise-comparison principles for some parabolic stochastic pdes.  The Annals of Probability, 45(1):377–403, 2017.
  * [29] Ida Kruk, Francesco Russo, and Ciprian A Tudor.  Wiener integrals, Malliavin calculus and covariance measure structure.  Journal of Functional Analysis, 249(1):92–142, 2007.
  * [30] Hubert Lacoin.  Influence of spatial correlation for directed polymers.  The Annals of Probability, 39(1):139–175, 2011.
  * [31] Jani Lukkarinen and Matteo Marcozzi.  Wick polynomials and time-evolution of cumulants.  Journal of Mathematical Physics, 57(8):083301, 2016.
  * [32] Ernesto Medina, Terence Hwa, Mehran Kardar, and Yi-Cheng Zhang.  Burgers equation with correlated noise: Renormalization-group analysis and applications to directed polymers and interface growth.  Physical Review A, 39(6):3053, 1989.
  * [33] Jean Mémin, Yulia Mishura, and Esko Valkeila.  Inequalities for the moments of Wiener integrals with respect to a fractional Brownian motion.  Statist. Probab. Lett., 51(2):197–206, 2001.
  * [34] Elchanan Mossel, Ryan O’Donnell, and Krzysztof Oleszkiewicz.  Noise stability of functions with low influences: invariance and optimality.  Ann. of Math. (2), 171(1):295–341, 2010.
  * [35] David Nualart.  The Malliavin calculus and related topics, volume 1995.  Springer, 2006.
  * [36] Vladas Pipiras and Murad S Taqqu.  Integration questions related to fractional Brownian motion.  Probability theory and related fields, 118(2):251–291, 2000.
  * [37] Guanglin Rang.  From directed polymers in spatial-correlated environment to stochastic heat equations driven by fractional noise in 1+ 1 dimensions.  Stochastic Processes and their Applications, 130(6):3408–3444, 2020\.
  * [38] Jay Rosen.  Random walks and intersection local time.  The Annals of Probability, pages 959–977, 1990.
  * [39] Carles Rovira and Samy Tindel.  On the Brownian-directed polymer in a Gaussian random environment.  Journal of Functional Analysis, 222(1):178–201, 2005.
  * [40] EL Rvaceva.  On domains of attraction of multi-dimensional distributions.  Select. Transl. Math. Statist. and Probability, 2:183–205, 1962\.
  * [41] Hao Shen, Jian Song, Rongfeng Sun, and Lihu Xu.  Scaling limit of a directed polymer among a Poisson field of independent walks.  Journal of Functional Analysis, 281(5):109066, 2021.
  * [42] Jian Song.  On a class of stochastic partial differential equations.  Stochastic Processes and their Applications, 127(1):37–79, 2017\.
  * [43] D Surgailis.  On Poisson multiple stochastic integrals and associated equilibrium Markov processes.  In Theory and application of random fields, pages 233–248. Springer, 1983.
  * [44] Ran Wei.  On the long-range directed polymer model.  Journal of Statistical Physics, 165(2):320–350, 2016.
  * [45] Gian-Carlo Wick.  The evaluation of the collision matrix.  Physical review, 80(2):268, 1950.


# Tuning the Many-body Interactions in a Helical Luttinger Liquid

Junxiang Jia  Elizabeth Marcellina  Division of Physics and Applied Physics,
School of Physical and Mathematical Sciences, Nanyang Technological
University, Singapore 637371, Singapore  Anirban Das  Department of Physics,
Indian Institute Of Technology Madras, Chennai, Tamil Nadu 600036, India
Center for Atomistic Modelling and Materials Design, Indian Institute Of
Technology Madras, Chennai, Tamil Nadu 600036, India  Michael S. Lodge
Division of Physics and Applied Physics, School of Physical and Mathematical
Sciences, Nanyang Technological University, Singapore 637371, Singapore
BaoKai Wang  Department of Physics, Northeastern University, Boston,
Massachusetts 02115, USA  Duc Quan Ho  Riddhi Biswas  Tuan Anh Pham  Wei Tao
Division of Physics and Applied Physics, School of Physical and Mathematical
Sciences, Nanyang Technological University, Singapore 637371, Singapore
Cheng-Yi Huang  Department of Physics, Northeastern University, Boston,
Massachusetts 02115, USA  Hsin Lin  Institute of Physics, Academia Sinica,
Taipei 115201, Taiwan  Arun Bansil  Department of Physics, Northeastern
University, Boston, Massachusetts 02115, USA  Shantanu Mukherjee  Department
of Physics, Indian Institute Of Technology Madras, Chennai, Tamil Nadu 600036,
India  Center for Atomistic Modelling and Materials Design, Indian Institute
Of Technology Madras, Chennai, Tamil Nadu 600036, India  Quantum Centre for
Diamond and Emergent Materials, Indian Institute of Technology Madras,
Chennai, Tamil Nadu 600036, India  Bent Weber<EMAIL_ADDRESS>Division of
Physics and Applied Physics, School of Physical and Mathematical Sciences,
Nanyang Technological University, Singapore 637371, Singapore  ARC Centre of
Excellence for Future Low-Energy Electronics Technologies (FLEET), School of
Physics & Astronomy, Monash University, Clayton 3800, Australia

Luttinger Liquid, quantum spin Hall, topological insulator, helical edge
state, Scanning tunnelling microscopy

††preprint: APS/123-QED

In one-dimensional (1D) systems, electronic interactions lead to a breakdown
of Fermi liquid theory and the formation of a Tomonaga-Luttinger Liquid (TLL)
haldane1981luttinger ; giamarchi2003quantum . The strength of its many-body
correlations can be quantified by a single dimensionless parameter, the
Luttinger parameter $K$ voit1995one , characterising the competition between
the electrons’ kinetic and electrostatic energies. Recently, signatures of a
TLL have been reported for the topological edge states of quantum spin Hall
(QSH) insulators li2015observation ; stuhler2020tomonaga , strictly 1D
electronic structures with linear (Dirac) dispersion and spin-momentum locking
kane2005quantum . Here we show that the many-body interactions in such helical
Luttinger Liquid can be effectively controlled by the edge state’s dielectric
environment. This is reflected in a tunability of the Luttinger parameter $K$,
distinct on different edges of the crystal, and extracted to high accuracy
from the statistics of tunnelling spectra at tens of tunneling points. The
interplay of topology and many-body correlations in 1D helical systems has
been suggested as a potential avenue towards realising non-Abelian
parafermions Orth2015 ; hsu2021helical ; Lodge2021 .

Figure 1: Edge states in monolayer 1T’-WTe2. a, STM topography of a monolayer
1T’-WTe2 crystal (scale bar: 5 nm). Insert: Crystal structure of the 1T’-WTe2
crystal structure with W atoms (blue) and Te atoms (orange). Blue arrows
represent different edge terminations. b,c, Corresponding differential
conductance (d$I$/d$V$) maps of the same island as in (a), measured in
constant height mode. d-f, Spatial profile of the measured local density of
state (LDOS) at 4.5 K, across the monolayer edge along the black arrow in (a),
and compared to a corresponding STM height profile (d). We extract a gap of
$\sim$70 meV in the 2D bulk. The data in (f) has been obtained by integration
of the LDOS in (e) between $E=-20$ meV to $E=-40$ meV and shows an overall
enhancement near the edge with superimposed charge density modulations also
observed in (b). The solid black line shows a fit to an exponential decay from
which we extract a decay length of $(2.1\pm 0.2)$ nm. Insert: Comparison of
point spectra at bulk and edge, highlighting the bulk gap (black) and LDOS
enhancement at the edge (blue) corresponding to the edge state. g, Evolution
of the edge state LDOS from the edge into the bulk (see corresponding markers
in (a,d)). The solid black lines are fits to TLL theory (see text). All
spectra have been offset for clarity, with the zero position indicated
(horizontal dashed lines).

A one-dimensional (1D) electronic system is described as a Tomonaga-Luttinger
liquid (TLL) rather than a Landau-Fermi liquid, since the strictly 1D nature
of the Fermi surface with only two Fermi points ($\pm k_{\mathrm{F}}$) leads
to a breakdown of screening haldane1981luttinger . As a consequence,
electronic correlations are strong, even for weak electron-electron
interactions, leading to collective excitations voit1995one . The elementary
excitations of a TLL can be described as bosons, rather than fermions, and are
given by charge and spin density waves with spin-charge separation
haldane1981luttinger ; giamarchi2003quantum . For weak interactions, the
correlation strength can be approximated by a single dimensionless parameter,
the Luttinger parameter bockrath1999balents ,

$K\approx\left[1+\frac{W(q=0)}{\pi\hbar v_{\mathrm{F}}}\right]^{-1/2}.$ (1)

Here, $\pi\hbar v_{\mathrm{F}}$ is the inverse density of states for 1D Dirac
electrons, with a constant Fermi velocity $v_{\mathrm{F}}$, and $W(q=0)$ is
the energy of electrostatic interactions in the long wavelength limit. The
former depends on the detail of the electronic band dispersion directly, while
the latter can be expected to be susceptible to the dielectric environment of
the 1D charge distribution, screening the many-body interactions. The
Luttinger parameter can thus be expected to be highly tunable, ranging within
its theoretical bounds, $0<K<1$, for repulsive (e.g. Coulombic) interactions,
where $K=1$ [$W(q=0)=0$] represents the limit of non-interacting electrons. A
threshold $K\approx 1/2$ is naturally given where the Coulomb and kinetic
energy terms become equal, thus defining $K<1/2$ as the limit of strong
interactions.

Luttinger liquids have been detected in several realisations of 1D electronic
structure, including linear chain organic conductors Schwartz1998 , metallic
single-wall carbon nanotubes (SWNTs) bockrath1999balents ; lee2004real ,
quantum Hall edge channels chang1996 , cold atoms paredes2004tonks , atomic
wires blumenstein2011atomically , and even crystalline mirror-twin boundaries
in semiconducting transition metal dichalcogenides jolie2019tomonaga .

A key experimental signature of a TLL is a pseudogap suppression in the local
density of states (LDOS), measured via tunnelling into the 1D electronic
system, e.g. from the tip of a scanning tunnelling microscope (STM). Such
pseudogap suppression is often referred to as zero-bias anomaly (ZBA) as it is
strictly centred at the Fermi energy $E_{\mathrm{F}}$ stuhler2020tomonaga ;
blumenstein2011atomically ; jolie2019tomonaga ; lee2004real . In a TLL, the
tunnelling conductivity around such ZBA scales universally, that is, as a
power law in both bias voltage and temperature, with the same scaling exponent
$\alpha$, related to the Luttinger parameter as
$\alpha=C\left(K+K^{-1}-2\right)$ Braunecker2012 . Different from the spin-
degenerate (spinful) parabolic dispersions of conventional 1D metallic systems
($C=1/4$), quantum spin Hall (QSH) insulators host linearly dispersing 1D edge
states, in which the spin polarity is locked to the crystal momentum
(helicity) kane2005quantum ; bernevig2006quantum . For a helical Luttinger
Liquid, strong spin orbit coupling causes the spin and chirality indices to
coincide, and the bosonized Hamiltonian can be expressed as an effectively
“spinless” Luttinger liquid with $C=1/2$. TLL formation in such 1D helical
systems is not only of fundamental interest as a testbed to study the
interplay of topology and strong correlations, but may further constitute a
potential host platform to realise non-Abelian parafermions qian2014quantum ;
zhang2014time ; Orth2015 ; hsu2021helical ; Lodge2021 .

Tuning of the many-body interaction strength has so far been demonstrated only
for spinful Luttinger liquids, for instance, by applied gate voltages
Kanda2004 ; Prokudina2014 , or magnetic fields Klanjsek2008 . In this work, we
demonstrate control of the Luttinger parameter $K$ in a helical TLL at the
edges of the quantum spin Hall insulator 1T’-WTe2, in which $K$ is found
distinct for different crystal edges, and can be controlled by dielectric
screening via the substrate. This becomes evident from a dependence of the
Luttinger parameter on substrate dielectric constant and crystal edge, ranging
between $K=(0.21\pm 0.01)$ and $K=(0.33\pm 0.01)$ in the strong coupling
limit. Our conclusions are based on temperature-dependent scanning tunnelling
spectroscopy, demonstrating universal scaling, combined with a detailed
statistical analysis of tens of tunneling points along WTe2 edge states,
showing that values of $K$ fall into normal distributions of distinct
statistical mean for different substrates and crystal edges.

The experimental extraction of $K$ is further supported by both a numerical
model to estimate the strength of screened Coulomb interactions for the edge
state’s charge distribution and a material-specific real-space tight-binding
model, which takes charge distribution and WTe2 edge dispersion self-
consistently into account. Both models are able to reproduce the trend of the
Luttinger parameter qualitatively in different screening environments.
Interestingly, despite being in the strong coupling limit ($K<0.5$), where
deviations from Eq. (1) may be expected due to strong interactions and a
renormalised Fermi velocity, we show that the perturbative approximation can
allow for qualitative predictions of its fundamental dependencies in terms of
the Coulomb interaction strength and edge state dispersion. In the limit of
very strong interactions ($K<1/4$) reached in some of our samples, two-
particle back-scattering of helical electrons may become significant
li2015observation ; Ziani2015 , potentially leading to additional effects such
as charge density wave formation, or even fractional Wigner crystallisation
Ziani2015 . Eq. (1) can thus serve as an important indicator for the strength
and tunability of Coulomb correlations in a Tomonaga-Luttinger liquid.

Figure 2: Demonstration of universal scaling for tunnelling into a Tomonaga-
Luttinger Liquid (TLL). a,b, Temperature dependent point spectra of the edge
state local density of states (LDOS) up to $T=25$ K, measured at Y-edge on
WTe2/HOPG and WTe2/BLG, respectively. The tunnelling spectra are offset by 0.6
on the $y$-axis with the zero position highlighted by horizontal dashed lines.
Solid black lines are fits to TLL theory, taking into account thermal
broadening (Section S3 of the Supplementary Information). The corresponding
bare power-law fits are shown as dotted curves. The shaded bands indicate the
conduction band tail. c,d, The data in (a) and (b) can be collapsed onto a
single universal curve with common scaling exponent $\alpha=0.45$ (WTe2/HOPG)
and $\alpha=0.98$ (WTe2/BLG). Shaded bands indicate statistical confidence
defined by the mean absolute error across the individual fits at different
measurement temperature. The shaded band in (d) indicates the conduction band
tail. e,f Temperature dependence of the measured zero-bias
($E=E_{\mathrm{F}}$) conductance in (a)-(d) (blue circles), compared with the
TLL model (solid black lines) and Coulomb gap (CG) model (solid red lines).
Fitting parameters of the CG model, such as the disorder parameter $\Omega$,
have been determined from a best fit to each data at $T=4.5$ K (see
Supplementary Fig. 4). The error bars are obtained from the Gaussian width of
the TLL fitting residual. Shaded bands indicate statistical confidence as
defined for (c,d). Figure 3: Dielectric screening of the TLL many-body
interactions. a,b, Histograms of the Luttinger parameter $K$, recorded across
69 locations along the X- and Y-edges of on WTe2/HOPG and WTe2/BLG islands.
Solid lines are fits to Gaussian normal distributions characterised to extract
statistical mean and standard deviation. Arrows indicate data points from
literature tang2017quantum ; chen2018large ; zhao2020strain ;
xu2018observation . c, Luttinger parameter $K$, extracted across different
samples, and literature, plotted as a function of the substrate dielectric
constant $\epsilon_{r}$. Red (blue) circles are the mean and standard error of
the mean extracted for the distributions in (a and b). Shaded bands and dashed
lines are numerical calculations of $K$ based on the electrostatic energy
$W(q=0)$ stored in a screened edge charge distribution of varying density (see
text), as well as the calculated edge dependent Fermi velocity $v_{F}$
extracted from (d) and (e). The bands have been modeled to coincide in width
with the extracted standard deviation of the distributions in a and b. The
solid red line represents tight-binding calculations of the Luttinger
parameter. Insert: Representative spectra measured at the Y- and X-edge,
respectively, on HOPG, with extracted Luttinger parameters $K\approx 0.35$
(Y-edge) and $K\approx 0.24$ (X-edge) from TLL fits shown as black lines. For
the dielectric constants, we assume $\epsilon_{r}=6$ (BLG/SiC)
bessler2019dielectric  and $\epsilon_{r}\rightarrow\infty$ (HOPG),
approximating the perfect metal limit). d,e, Tight-binding band structure of
the $\Gamma$-Y and $\Gamma$-X directions, respectively, assuming the atomic
terminations shown by coloured arrows in (insert to d). Grey shading indicates
the 2D bulk band edges with a band gap comparable to the data in Fig. 1e,f.
Calculations are based on the edge terminations highlighted in the insert to
d. Insert to d, 1T’-WTe2 crystal structure, showing the position of W atoms
(blue) and Te atoms (orange). Insert to e, Corresponding WTe2 Brillouin zone.

Helical TLLs have been reported for the 1D topological edge states of
InAs/GaSb semiconductor heterostructures li2015observation  and, more
recently, for the edge states of atomically-flat bismuthene
stuhler2020tomonaga . An alternative realisation of an atomically-thin QSH
insulator Lodge2021 ; hsu2021helical  is 1T’-WTe2 qian2014quantum ;
fei2017edge ; tang2017quantum ; jia2017direct ; Wu2018 , with topological band
gaps of a few tens of meV tang2017quantum ; jia2017direct , tunable by
electric fields, local strain. A typical 1T’-WTe2 monolayer crystal is shown
in Fig. 1a, synthesised by van-der-Waals molecular-beam epitaxy (MBE)
tang2017quantum ; jia2017direct  on highly-oriented pyrolytic graphite (HOPG)
(see Methods). We achieve typical island circumferences $L=100-400$ nm with
disordered edges. Locally, however, we find straight and well-defined
sections, with typical lengths between $10-20$ nm that are terminated either
perpendicular (Y-edge) or parallel (X-edge) to the atomic rows.

Corresponding differential conductance (d$I$/d$V$) maps of the same area are
shown in Fig. 1b,c, reflecting the expected QSH electronic signature of an
enhanced local density of states (LDOS) at the crystal edge around an
insulating 2D bulk tang2017quantum ; collins2018electric . A larger area of
the same crystal can be viewed in Supplementary Fig. 1, showing that the
electronic edge state is largely continuous along the island’s circumference,
even in the presence of local edge roughness and disorder, and regardless of
edge direction or termination, and at adjoining corners. The relevant length-
scale determining the electrostatic energy of the screened Coulomb
interactions within the edge state’s 1D charge distribution is therefore set
by the islands’ circumference. Superimposed spatial modulations in the charge
density can be shown to arise from edge roughness even in otherwise perfect
monolayer crystals (Supplementary Fig. 7). However, as these modulations –
more clearly pronounced on the X-edges – appear quasi-periodic, additional
mechanisms may be at play, including spatial inhomogeneity of the Fermi
velocity and/or interaction-mediated scattering of the edge state, leading to
charge density wave order.

We note that we do observe semi-metallic behavior and a variety of in-gap
states in some monolayer crystals, especially in the vicinity of lattice
defects and adsorbates (Supplementary Fig. 1). Some crystals also show a small
but finite 2D LDOS within the bulk and slight variations in gap size. For this
study, however, we focus on larger clean regions, such as shown in Fig. 1a-c,
in which a well-defined suppression of the bulk LDOS is observed.

A spatial profile of the local density of states (LDOS), measured along the
black arrow in Fig.1a, is shown in Fig. 1e, alongside a corresponding STM
height profile (Fig. 1d). The data confirms a well-defined energy gap over
$\sim 70$ meV in the 2D bulk, similar in magnitude to earlier reports
tang2017quantum ; jia2017direct , and in agreement with the tight-binding
model adopted in this work lau2019influence  (compare Fig. 3d,e). The gap is
seen to narrow and eventually disappear at the edge, consistent with the
presence of a metallic edge state. We identify the position of the edge,
electronically, from a profile of the energy-integrated LDOS (-20 to -40 meV)
in Fig. 1f, which shows a finite width $\Delta x\approx 2$ nm with
exponentially decay into the 2D bulk. From an exponential fit, we extract a
decay length of $(2.1\pm 0.2)$ nm, significantly larger compared to the
crystal lattice constant ($a=0.348$ nm). Spatial modulations, superimposed
onto the overall edge state LDOS enhancement further indicate that the
observed LDOS enhancement is not due to localised atomic orbitals, but rather
are an electronic effect due to the metallic edge state.

The difference between bulk and edge state LDOS is further highlighted by the
direct comparison of individual STS point spectra shown in the insert to Fig.
1f. For a linear (Dirac) dispersion in 1D, one would expect a constant
(energy-independent) edge state LDOS throughout the gap. Instead, we observe a
clear V-shaped suppression with minimum at the Fermi energy $E=E_{\mathrm{F}}$
(ZBA). This ZBA – a first indication for the presence of 1D correlations
tang2017quantum ; collins2018electric ; stuhler2020tomonaga  – is strictly
confined to the crystal edges as confirmed in Fig. 1g, where we plot
individual d$I$/d$V$ point spectra away from the Y-edge into the monolayer
bulk.

Black solid lines are fits to the theoretically predicted tunnelling density
of states into a TLL at 4.5 K, following the universal scaling relationship
bockrath1999balents ,

$\rho_{\text{TLL}}(E,T)\propto
T^{\alpha}\cosh\left(\frac{E}{2k_{\mathrm{B}}T}\right)\left|\Gamma\left(\frac{1+\alpha}{2}+i\frac{E}{2\pi
k_{\mathrm{B}}T}\right)\right|^{2}$ (2)

Here, $\alpha$ is the scaling exponent, $E-E_{\mathrm{F}}\equiv eV$ is the STM
bias voltage, $T$ is temperature, and $\Gamma\left(x\right)$ is the gamma
function. For our fits, we further take broadening of the LDOS into account,
consistent with the procedure in Ref. stuhler2020tomonaga  (as shown in
Supplementary Fig. 2).

Particular to 1D helical systems, the linear (Dirac) dispersion supports the
bosonisation picture of TLL theory over a much larger energy range.
Contrasting the spinful parabolic dispersions of conventional metals, which
deviate from a linear approximation already at moderately low energy, the
1T’-WTe2 edge state dispersion remains linear over a wide energy range within
the QSH band gap lau2019influence  (compare Fig. 3d-e), consistent with the
range of our TLL fits (typically -20 meV to +15 meV). From 33 edge spectra in
(2 of which are shown in Fig. 1g), we extract $\alpha=0.59\pm 0.02$ and
$K=0.35\pm 0.02$, assuming helicity of the QSH edge state ($C=1/2$), in good
agreement with values recently reported for the helical edges of bismuthene
($K\approx 0.42$) stuhler2020tomonaga , which features similarly tightly
confined 1D edge states, in contrast to the larger edge state decay lengths of
semiconductor-based QSH systems li2015observation .

We note that in principle, a narrow suppression in the LDOS of 1D electronic
systems haldane1981luttinger ; bockrath1999balents ; tang2017quantum ;
collins2018electric ; stuhler2020tomonaga  may have its origin in a number of
different mechanisms, most notably the formation of a Coulomb gap (CG) due to
long-range interactions bartosch2002zero . However, as we show in
Supplementary Fig. 4 and 5, a Coulomb gap (CG) model cannot explain the
detailed functional form of the ZBA observed. To further rule out the presence
of a CG, we turn to temperature-dependent scanning probe spectroscopy, to
confirm universal scaling of the TLL tunnelling conductance. In Fig. 2a-d, we
compare d$I$/d$V$ spectra measured on WTe2/HOPG with those taken on WTe2/BLG.
We limit the temperature range explored to $T<25$ K, as higher temperatures
are found to give rise to a significant bulk contribution to the overall LDOS
(Supplementary Fig. 5e). Black solid lines in Fig. 2a,b are fits to the
universal scaling relationship (Eq. 2), yielding comparable power-law
exponents $\alpha\approx 0.45$ (HOPG) and $\alpha\approx 0.98$ (BLG) at all
temperatures. The corresponding bare power law is indicated by dotted lines.
Universal scaling is confirmed by normalising the measured d$I$/d$V$ with
$T^{\alpha}$ and plotting the result as a function of the dimensionless energy
parameter $eV/k_{\mathrm{B}}T$. Following this procedure, all spectra collapse
onto a single universal curve (Fig. 2c,d) with power-law exponents
$\alpha=0.45\pm 0.03$ ($K=0.40\pm 0.01$) for WTe2/HOPG and $\alpha=0.98\pm
0.03$ ($K=0.27\pm 0.01$) for WTe2/BLG. The uncertainty in the fitted values of
$\alpha$, indicated by the shaded confidence bands in the figure, arises from
the mean absolute error of TLL fits across the different measurement
temperatures. As a further confirmation of universal scaling, we plot the
zero-bias ($E=E_{\mathrm{F}}$) tunnelling conductance as a function of
temperature in Fig. 2e-f, comparing the expected temperature dependencies of
TLL and CG theories, respectively. The temperature dependence of the TLL is in
excellent agreement with the universal scaling exponents extracted from Fig.
2a-d, while the CG model shows substantial deviations from the measured data
on both substrates, especially at temperatures $<10$ K.

Universal scaling as demonstrated is usually considered sufficient proof for
the presence of a TLL stuhler2020tomonaga . Further to this, however, below we
show that the observed differences in the Luttinger parameter across
substrates can be explained by substrate-induced screening the many-body
interactions. To this end, we have performed a detailed statistical analysis
of $K$ values extracted from a total of 69 locations, measured along the edges
states of samples on both substrates, plotted in the histograms of Fig. 3a,b.
We find that $K$ follows normal distributions with distinct statistical mean
but comparable standard deviation for both edge terminations and substrates.
This is further confirmed by fits to Gaussian normal distributions (solid
lines), from which we extract $K=0.27\pm 0.01$ (Y-edge) and $K=0.21\pm 0.01$
(X-edge) for WTe2/BLG and $K=0.33\pm 0.01$ (Y-edge) and $K=0.24\pm 0.01$
(X-edge) for WTe2/HOPG. The uncertainty is given by the standard error of the
mean $\sigma/\sqrt{n}$ which quantifies the accuracy of the extracted mean,
reflecting that $K$ was obtained from a population of $n$ sample points with
standard deviation $\sigma$. Further detail on the statistical analysis,
including normality tests, and statistical significance of the differences in
mean is listed in Supplementary Table 1. Finally, we include $K$ values
extracted from edge state tunnelling spectra in literature tang2017quantum ;
zhao2020strain ; xu2018observation ; chen2018large , and find that their
values fall within the same distributions as shown in Fig. 3a,b.

Mean and standard error of the mean for both edges are plotted as a function
of the substrate dielectric constant $\epsilon_{r}$ in Fig. 3c. Across all
data, the Luttinger parameter $K$ increases at low $\epsilon_{r}$ and
saturates for $\epsilon_{r}\rightarrow\infty$, where the latter corresponds to
HOPG Gramse2012 . We note that although there have also been reports of a
finite dielectric constant for HOPG ($\epsilon_{r}\approx 53$)
hong2015exfoliated , $K$ starts to saturate at ($\epsilon_{r}\approx 50$),
implying that $\epsilon_{r}\rightarrow\infty$ provides a reasonable
approximation in this context.

To explain this trend, we devise a qualitative model for $K$ teo2009critical ,
based on the perturbative approximation (Eq. (1)) and the assumption of a
substrate-dependent mirror charge, screening many-body Coulomb interactions
within the helical edge. In our model, edge state and mirror charge are
approximated by space-charge distributions $\rho({x,y,z})$ and
$\rho^{\prime}(x,y,z)=-\frac{\epsilon_{r}-1}{\epsilon_{r}+1}\rho(x,y,-z)$,
respectively, allowing to calculate the electrostatic energy $W(q=0)$ of the
total charge numerically (see Section S6 of the Supplementary Information). In
order to estimate $K$ for the respective Y- (X-) edges, we extract the Fermi
velocities $v_{\mathrm{F}}=1.16\times 10^{5}$ m/s ($2.69\times 10^{5}$ m/s)
from the tight-binding dispersions in Fig. 3b-e. The red and blue shaded bands
in Fig. 3c have been calculated for an average edge length $L=200$ nm, and
space-charge densities $\rho_{0,y}=0.0039\pm 0.0007$ e/nm3 and
$\rho_{0,x}=0.0082\pm 0.0009$ e/nm3) to reflect the standard deviations
consistent with Fig. 3a,b. The observed trend is further confirmed by our
material-specific tight-binding model, taking edge state charge density and
band dispersion self-consistently into account (Supplementary Fig. 7). The
solid red line in Fig. 3a represents the calculated the Luttinger parameter
$K$ along the Y-edge.

The ability of our model to reproduce the trends in $K$ may seem surprising at
first, given that Eq. (1) should strictly only apply in the perturbative limit
of weak interactions ($K>0.5$). However, we note that several factors can
extend the range of validity of Eq. (1), especially in helical systems. For
instance, Coulomb interactions do not allow for backscattering with large
momentum transfer $q=2k_{\mathrm{F}}$ in a helical TLL. This supports a
perturbative formalism which primarily includes the long wavelength forward
scattering modes ($q\approx 0$). More so, our observation of a finite width
TLL in which the LDOS falls exponentially from the edge into the QSH bulk
(compare Fig. 1f) leads to a Coulomb interaction extending up to the fifth
nearest neighbour ion of the lattice (Supplementary Fig. 8). This justifies
the approximation $g_{2}=g_{4}$ implied in our perturbative analysis where
$g_{2}$ and $g_{4}$ correspond to interaction terms in a Hamiltonian
contribution,
$g_{4}(\psi_{R\uparrow}^{\dagger}\psi_{R\uparrow})^{2}+g_{2}(\psi_{R\uparrow}^{\dagger}\psi_{R\uparrow}\psi_{L\downarrow}^{\dagger}\psi_{L\downarrow})$,
where $\psi_{R\uparrow}$ annihilates a right mover electron with spin up, and
$\psi_{L\downarrow}$ annihilates a left mover electron with spin down. We
conclude that while the precise quantitative match of $K$ may be fortuitous,
the perturbative treatment is able to correctly predict the trend in $K$, for
different edge dispersions and substrate mediated screening. This work
therefore presents a precision experimental determination of the effect of
kinetic and potential energy terms on the Luttinger parameter, and may form a
benchmark for future theoretical studies.

To summarise, we have demonstrated tunability of the many-body interactions in
a Tomonaga-Luttinger Liquid (TLL) at the helical edges of the quantum spin
Hall insulator 1T’-WTe2. Synthesising WTe2 on metallic van-der-Waals
substrates, we have been able to demonstrate effective screening of the many-
body Coulomb interactions, reflected in a significant enhancement of the
Luttinger parameter. The cleanest signatures of the TLL are observed on the
well-screened Y-edges of WTe2/HOPG, where $K=(0.33\pm 0.01)$. Here, the higher
degree of tunability observed would suggest that systems of lower Fermi
velocity and low carrier density may allow to explore a wider range of $K$,
something which may be confirmed in other QSH material systems. In turn, the
strongest many-body interactions are found on the unscreened X-edges for
WTe2/BLG where we measure $K=(0.21\pm 0.01)$ – the lowest reported for any
helical TLL system to date. Indeed, this value represents the limit of very
strong interactions ($K<0.25$), where additional mechanisms, such as charge
density wave formation due to two-particle scattering or even Wigner
crystallisation Ziani2015  may be expected. Indeed, the quasi-periodic charge
density modulations observed on X-edges are not inconsistent with such
phenomena and will be a focus for future investigation into the role of edge
disorder and scattering. We thus expect that our results will stimulate
further experimental and theoretical investigations of the interplay of
topology and strong electronic interactions in 1D, also in the superconducting
state, in which non-Abelian parafermions are predicted hsu2021helical .

##  I Methods

MBE growth. Monolayer 1T’-WTe2 was synthesised by molecular-beam epitaxy (MBE)
on highly oriented pyrolytic graphite (HOPG) and bilayer graphene (BLG) grown
on 6H-SiC(0001) substrates in an Omicron Lab10 ultrahigh vacuum (UHV) MBE
chamber with a base pressure of $2\times 10^{-10}$ mbar. Prior to MBE, the
HOPG substrates were mechanically cleaved in ambient conditions and then
degassed at 300 overnight. The 6H-SiC(0001) substrates were degassed at 600
overnight in UHV and then flash annealed at 2000 for a few cycles to
graphitise producing BLG. The WTe2 islands were grown by evaporating high-
purity W(99.998$\%$) and Te(99.999$\%$) with a flux ratio of 1:280 and
substrate temperature of 165 for 2 hours to cover approximately 50% of the
monolayer area. In total, we measured three 1T’-WTe2/HOPG samples and another
three of 1T’-WTe2/BLG.

Scanning tunnelling spectroscopy. Low-temperature scanning tunnelling
microscopy and spectroscopy (STM/STS) measurements were performed in an
Omicron low-temperature STM ($\simeq 4.5$ K) under UHV conditions ($\simeq
5\times 10^{-11}$ mbar). For temperature-dependent STM/STS, a heater element
fitted to the sample stage was used for counter-heating from 4.5 K. For all
spectroscopy measurements, we used a platinum/iridium tip calibrated against
the Au(111) Shockley surface state. The spectroscopy measurements were carried
out using standard lock-in techniques with a modulation amplitude of
$V_{\mathrm{ac}}=2$ mV and a modulation frequency of 1.331 kHz. The value of
$V_{\mathrm{ac}}$ was optimised to maximize the signal-to-noise ratio while
eliminating instrumental broadening (Supplementary Fig. 3). Differential
conductance maps were taken in constant height mode.

## Data availability

The scanning tunneling spectroscopy data generated in this study have been
deposited in the NTU research data repository DR-NTU (Data) database at
https://doi.org/10.21979/N9/J8BWVQ.

##  II Code availability

The computer code used for data analysis and modelling is available from BW or
SM upon reasonable request.

## References

  * (1) Haldane, F.  Luttinger liquid theory of one-dimensional quantum fluids. I. Properties of the Luttinger model and their extension to the general 1D interacting spinless Fermi gas.  _J. Phys. C: Solid State Phys._ 14, 2585 (1981).
  * (2) Voit, J.  One-dimensional Fermi liquids.  _Rep. Prog. Phys._ 58, 977 (1995).
  * (3) Giamarchi, T.  _Quantum physics in one dimension_ , vol. 121 (Clarendon press, 2003).
  * (4) Bockrath, M. _et al._ Luttinger-liquid behaviour in carbon nanotubes.  _Nature_ 397, 598–601 (1999).
  * (5) Schwartz, A. _et al._ On-chain electrodynamics of metallic $(\mathrm{TMTSF}{)}_{2}$x salts: Observation of tomonaga-luttinger liquid response.  _Phys. Rev. B_ 58, 1261–1271 (1998).
  * (6) Lee, J. _et al._ Real space imaging of one-dimensional standing waves: Direct evidence for a Luttinger liquid.  _Phys. Rev. Lett._ 93, 166403 (2004).
  * (7) Chang, A. M., Pfeiffer, L. & West, K.  Observation of chiral Luttinger behavior in electron tunneling into fractional Quantum Hall Edges.  _Phys. Rev. Lett._ 77, 2538 (1996).
  * (8) Paredes, B. _et al._ Tonks–Girardeau gas of ultracold atoms in an optical lattice.  _Nature_ 429, 277–281 (2004).
  * (9) Blumenstein, C. _et al._ Atomically controlled quantum chains hosting a Tomonaga–Luttinger liquid.  _Nat. Phys._ 7, 776–780 (2011).
  * (10) Jolie, W. _et al._ Tomonaga-Luttinger liquid in a box: electrons confined within MoS2 mirror-twin boundaries.  _Phys. Rev. X_ 9, 011055 (2019).
  * (11) Stühler, R. _et al._ Tomonaga–Luttinger liquid in the edge channels of a quantum spin Hall insulator.  _Nat. Phys._ 16, 47–51 (2020).
  * (12) Braunecker, B., Bena, C. & Simon, P.  Spectral properties of luttinger liquids: A comparative analysis of regular, helical, and spiral luttinger liquids.  _Phys. Rev. B_ 85, 035136 (2012).
  * (13) Kane, C. L. & Mele, E. J.  Quantum spin Hall effect in graphene.  _Phys. Rev. Lett._ 95, 226801 (2005).
  * (14) Bernevig, B. A., Hughes, T. L. & Zhang, S.-C.  Quantum spin Hall effect and topological phase transition in HgTe quantum wells.  _Science_ 314, 1757–1761 (2006).
  * (15) Qian, X., Liu, J., Fu, L. & Li, J.  Quantum spin Hall effect in two-dimensional transition metal dichalcogenides.  _Science_ 346, 1344–1347 (2014).
  * (16) Zhang, F. & Kane, C. L.  Time-reversal-invariant Z4 fractional Josephson effect.  _Phys. Rev. Lett._ 113, 036401 (2014).
  * (17) Orth, C. P., Tiwari, R. P., Meng, T. & Schmidt, T. L.  Non-Abelian parafermions in time-reversal-invariant interacting helical systems.  _Phys. Rev. B_ 91, 081406(R) (2015).
  * (18) Hsu, C.-H., Stano, P., Klinovaja, J. & Loss, D.  Helical Liquids in Semiconductors.  _Semicond. Sci. Technol._ 36, 12003 (2021).
  * (19) Lodge, M. S., Yang, S. A., Mukherjee, S. & Weber, B.  Atomically Thin Quantum Spin Hall Insulators.  _Adv. Mater._ 2008029 (2021).
  * (20) Kanda, A., Tsukagoshi, K., Aoyagi, Y. & Ootuka, Y.  Gate-Voltage Dependence of Zero-Bias Anomalies in Multiwall Carbon Nanotubes.  _Phys. Rev. Lett._ 92, 036801 (2004).
  * (21) Prokudina, M. G. _et al._ Tunable Nonequilibrium Luttinger Liquid Based on Counterpropagating Edge Channels.  _Phys. Rev. Lett._ 112, 216402 (2014).
  * (22) Klanšek, M. _et al._ Controlling Luttinger Liquid Physics in Spin Ladders under a Magnetic Field.  _Phys. Rev. Lett._ 101, 137207 (2008).
  * (23) Li, T. _et al._ Observation of a helical Luttinger liquid in InAs/GaSb quantum spin Hall edges.  _Phys. Rev. Lett._ 115, 136804 (2015).
  * (24) Ziani, N. T., Crépin, F. & Trauzettel, B.  Fractional Wigner Crystal in the Helical Luttinger Liquid.  _Phys. Rev. Lett._ 115, 206402 (2015).
  * (25) Fei, Z. _et al._ Edge conduction in monolayer WTe2.  _Nat. Phys._ 13, 677–682 (2017).
  * (26) Tang, S. _et al._ Quantum spin Hall state in monolayer 1T’-WTe2.  _Nat. Phys._ 13, 683–687 (2017).
  * (27) Jia, Z.-Y. _et al._ Direct visualization of a two-dimensional topological insulator in the single-layer 1T’-WTe2.  _Phys. Rev. B_ 96, 041108 (2017).
  * (28) Wu, S. _et al._ Observation of the quantum spin Hall effect up to 100 kelvin in a monolayer crystal.  _Science_ 359, 76 (2018).
  * (29) Collins, J. L. _et al._ Electric-field-tuned topological phase transition in ultrathin Na3Bi.  _Nature_ 564, 390–394 (2018).
  * (30) Lau, A., Ray, R., Varjas, D. & Akhmerov, A. R.  Influence of lattice termination on the edge states of the quantum spin Hall insulator monolayer 1T‘-WTe2.  _Phys. Rev. Mater._ 3, 054206 (2019).
  * (31) Bartosch, L. & Kopietz, P.  Zero bias anomaly in the density of states of low-dimensional metals.  _Eur. Phys. J. B_ 28, 29–36 (2002).
  * (32) Zhao, C. _et al._ Strain Tunable Semimetal–Topological-Insulator Transition in Monolayer 1T’-WTe2.  _Phys. Rev. Lett._ 125, 046801 (2020).
  * (33) Xu, H. _et al._ Observation of gap opening in 1T’ Phase MoS2 nanocrystals.  _Nano Lett._ 18, 5085–5090 (2018).
  * (34) Chen, P. _et al._ Large quantum-spin-Hall gap in single-layer 1T’ WSe2.  _Nat. Commun._ 9, 1–7 (2018).
  * (35) Gramse, G., Gomila, G. & Fumagalli, L.  Quantifying the dielectric constant of thick insulators by electrostatic force microscopy: effects of the microscopic parts of the probe.  _Nanotechnology_ 23, 205703 (2012).
  * (36) Hong, X. & Chung, D.  Exfoliated graphite with relative dielectric constant reaching 360, obtained by exfoliation of acid-intercalated graphite flakes without subsequent removal of the residual acidity.  _Carbon_ 91, 1–10 (2015).
  * (37) Teo, J. C. & Kane, C.  Critical behavior of a point contact in a quantum spin Hall insulator.  _Phys. Rev. B_ 79, 235321 (2009).
  * (38) Bessler, R., Duerig, U. & Koren, E.  The dielectric constant of a bilayer graphene interface.  _Nanoscale Advances_ 1, 1702–1706 (2019).

##  III Acknowledgement

This research is supported by the National Research Foundation (NRF)
Singapore, under the Competitive Research Programme “Towards On-Chip
Topological Quantum Devices” (NRF-CRP21-2018-0001) (BW) with further support
from the Singapore Ministry of Education (MOE) Academic Research Fund Tier 3
grant (MOE2018-T3-1-002) “Geometrical Quantum Materials” (BW). The work at
Northeastern University was supported by the US Department of Energy (DOE),
Office of Science, Basic Energy Sciences Grant No. DE-SC0022216 and benefited
from Northeastern University’s Advanced Scientific Computation Center and the
Discovery Cluster and the National Energy Research Scientific Computing Center
through DOE Grant No. DE-AC02-05CH11231 (AB). HL acknowledges the support by
the Ministry of Science and Technology (MOST) in Taiwan under grant number
MOST 109-2112-M-001-014-MY3. SM acknowledges financial support from SERB grant
number SP/2021/1008/PH/SERB/008839. BW acknowledges a Singapore National
Research Foundation (NRF) Fellowship (NRF-NRFF2017-11).

##  IV Author Contributions

JJ, DQH, WT, and TAP grew the WTe2 monolayers. JJ, DQH, and MSL performed the
scanning tunnelling spectroscopy experiments. AD, JJ, RB, CYH, BKW, HL, AB and
SM performed the theoretical calculations. JJ, EM, and BW analyzed the data.
HL and SM supervised the theory work. BW conceived and coordinated the
project. JJ, EM, SM and BW wrote the manuscript with input from all authors.

##  V Competing Interests

The authors declare no competing interests.


# Cross-correlation image analysis for real-time particle tracking

Leonardo R. Werneck  Department of Physics, University of Idaho, Moscow, ID
83843, USA  Cody Jessup  Department of Physics, Montana State University,
Bozeman, MT 59717, USA  Austin Brandenberger  Department of Physics, Montana
State University, Bozeman, MT 59717, USA  Tyler Knowles  Department of
Mathematics, West Virginia University, Morgantown, WV 26506, USA  Charles W.
Lewandowski  Department of Physics, Montana State University, Bozeman, MT
59717, USA  Space Dynamics Laboratory, Albuquerque, NM 87106, USA  Megan Nolan
Department of Physics, Montana State University, Bozeman, MT 59717, USA  Ken
Sible  Department of Physics and Astronomy, West Virginia University,
Morgantown, WV 26506, USA  Center for Gravitational Waves and Cosmology,
Chestnut Ridge Research Building, Morgantown, WV 26506, USA  Department of
Computer Science and Engineering, University of Notre Dame, South Bend, IN
46556, USA  Zachariah B. Etienne  Department of Physics, University of Idaho,
Moscow, ID 83843, USA  Department of Physics and Astronomy, West Virginia
University, Morgantown, WV 26506, USA  Center for Gravitational Waves and
Cosmology, Chestnut Ridge Research Building, Morgantown, WV 26506, USA  Brian
D’Urso  Department of Physics, Montana State University, Bozeman, MT 59717,
USA

###### Abstract

Accurately measuring translations between images is essential in many fields,
including biology, medicine, geography, and physics. Existing methods,
including the popular FFT-based cross-correlation, are not suitable for real-
time analysis, which is especially vital in feedback control systems. To fill
this gap, we introduce a new algorithm which approaches shot-noise limited
displacement detection and a GPU-based implementation for real-time image
analysis.

††journal: opticajournal

##  1 Introduction

Tracking the motion of particles in a sequence of images over time is a common
practice in various fields, ranging from biology [1] and geography [2] to
physics [3]. For applications requiring real-time feedback control, the
center-of-mass algorithm is often employed due to its simplicity [4]. However,
this method is primarily designed for measuring the position of a bright spot
against a dark background, making it susceptible to biases [5].

Accurately tracking a single particle against a relatively featureless
background can be achieved by measuring the apparent displacement of the
image’s content as the particle moves [6, 3]. Various codes for generalized
image analysis are available [7, 8, 9, 10, 11], but are typically not suitable
for real-time tracking.

In this letter, we introduce a real-time, cross-correlation-based image
analysis method that approaches the shot-noise limit of particle tracking. We
also provide an open-source implementation of this algorithm and demonstrate
its capabilities by tracking a microsphere levitated in a magneto-
gravitational trap [12]. Our implementation is suitable for commercial off-
the-shelf graphics processing units (GPUs), and can analyze hundreds of images
per second.

##  2 Basic Approach

A common method to find the location $\vec{R}$ of a bright spot on a dark
background is the center of mass (CM) calculation,
$\vec{R}=\sum_{\vec{r}}\vec{r}\,I(\vec{r})/\sum_{\vec{r}}I(\vec{r})$, with
$\vec{r}$ denoting the location of a pixel in an image $I(\vec{r})$. In the
ideal case of a nearly featureless background, this method can approach the
shot noise limit [13]. However, it is highly affected by background light and
image boundaries [5]. A more statistically robust approach, maximum likelihood
estimation (MLE), adjusts parameters in a fit distribution to maximize the
likelihood that the data comes from the distribution with Poisson statistics.
If Gaussian statistics are used as an approximation, we can instead minimize
$\chi^{2}$ relative to the displacement $\vec{r}_{0}$ between two images, as
expressed in the equation

$\chi^{2}=\sum\limits_{\vec{r}}\Biggl{[}\frac{I(\vec{r}-\vec{r}_{0})-E(\vec{r})}{\sigma(\vec{r},\vec{r}_{0})}\Biggr{]}^{2},$
(1)

where $I(\vec{r}-\vec{r}_{0})$ is the image $I(\vec{r})$ translated by
$\vec{r}_{0}$, $E(\vec{r})$ is the reference image (RI; see Sec. 3), and
$\sigma(\vec{r},\vec{r}_{0})$ models the noise in the images. Minimizing
$\chi^{2}$ directly is computationally intensive, so more efficient
strategies, such as cross-correlation (CC), are generally preferred.

The uniformly-weighted CC (CC-U) method, for example, assumes that
$\sigma(\vec{r},\vec{r}_{0})$ has a constant value denoted by $\sigma$. Under
this assumption, Eq. (1) simplifies to:

$\chi^{2}=\frac{1}{\sigma^{2}}\sum\limits_{\vec{r}}\Bigl{[}I(\vec{r}-\vec{r}_{0})^{2}-2I(\vec{r}-\vec{r}_{0})E(\vec{r})+E(\vec{r})^{2}\Bigr{]}\;.$
(2)

Since the sums over $I(\vec{r}-\vec{r}_{0})^{2}$ and $E(\vec{r})^{2}$ span the
entire images, they are independent of $\vec{r}_{0}$. The displacement
$\vec{r}_{0}$ is then obtained by finding the position of the maximum of the
CC between the shifted image $I(\vec{r}-\vec{r}_{0})$ and the RI $E(\vec{r})$,
i.e.,

$I{\star}E=\mathcal{F}^{-1}\Bigl{[}\mathcal{F}\bigl{[}I\bigr{]}\otimes\overline{\mathcal{F}\bigl{[}E\bigr{]}}\Bigr{]}\;,$
(3)

where $\mathcal{F}$ is the Fourier transform, $\mathcal{F}^{-1}$ its inverse,
$\otimes$ denotes element-wise multiplication, and an overline denotes complex
conjugation.

For CMOS cameras with adequate light levels, photon shot noise is often the
dominant noise source. The noise in each pixel follows a Poisson distribution,
which can be approximated by a Gaussian distribution with standard deviation
equal to the square root of the photon count in the pixel. In this case, a
shot-noise-weighted version of the CC method (CC-SN) is more appropriate. This
method assumes the RI provides an estimate for the _average_ pixel values of
the images in a data set, i.e.,

$\sigma(\vec{r},\vec{r}_{0})=\sqrt{E(\vec{r})}\;.$ (4)

For nearly zero pixel values, the square root poorly approximates the
uncertainty of the Poisson distribution. Negative pixel values may also arise
when subtracting the average mean background from images. To compensate, we
add a constant offset $\varepsilon$ to the new images (see Sec. 3 for how this
affects the RI), with larger offset values de-emphasizing darker pixels.

Adding the offset to the new image and using Eq. (4) in Eq. (1) we find

$\chi^{2}=\sum\limits_{\vec{r}}\Bigl{\\{}E(\vec{r})-2I^{\prime}(\vec{r}-\vec{r}_{0})+\bigl{[}I^{\prime}(\vec{r}-\vec{r}_{0})\bigr{]}^{2}/E(\vec{r})\Bigr{\\}}\;,$
(5)

where $I^{\prime}\equiv I+\varepsilon$. As with Eq. (2), the first two sums in
Eq. (5) are independent of $\vec{r}_{0}$, and the minimization only depends on
the CC between $I^{\prime 2}$ and the reciprocal of $E$. Therefore, the
displacement $\vec{r}_{0}$ that minimizes Eq. (5) is determined by the
position of the maximum of

$I^{\prime
2}{\star}E^{-1}=\mathcal{F}^{-1}\Bigl{[}\mathcal{F}\bigl{[}I^{\prime
2}\bigr{]}\otimes\overline{\mathcal{F}\bigl{[}1/E\bigr{]}}\Bigr{]}\;.$ (6)

##  3 Numerical Methods

The location of the maximum of the CC gives an integer pixel estimate of the
displacement between two images. Following [6, 7], we define a 1.5-pixel
square region centered on the CC’s maximum and sample it with $1.5u$ points in
both the horizontal and vertical directions, where $u$ is the upsampling
factor. Data from the element-wise product on the right-hand sides of Eqs. (3)
or (6) (depending on the method) is then interpolated to the upsampled region,
followed by an inverse discrete Fourier transform to generate an upsampled CC.
The location of the maximum in the upsampled CC yields a sub-pixel estimate of
the displacement.

The CC-U and CC-SN methods can be used to analyze images both in real-time or
with an iterative method that uses past and future images, with different
approaches for updating the RI in each case. The following discussion focuses
on the CC-U method, but it can be easily adapted for the CC-SN method by
replacing $I_{i}\to 1/I_{i}^{\prime}$ and $E_{n}\to 1/E_{n}$.

###  3.1 Offline Analysis

In the offline case [3], we iteratively update the RI as the average in
Fourier space of all images within a data set
$\bigl{\\{}I_{0},I_{1},\ldots,I_{N-1}\bigr{\\}}$. For each complete pass
through the images, the Fourier transform of the RI at the $n$-th iteration,
$\mathcal{F}[E_{n}]$, is computed using

$\mathcal{F}[E_{n}]=\frac{1}{N}\sum\limits_{i=0}^{N-1}\mathcal{F}\bigl{[}I_{i}(\vec{r}-\vec{r}_{0})\bigr{]}\;,$
(7)

where $\mathcal{F}[E_{0}]=\mathcal{F}[I_{0}]$. For a given iteration,
computing the displacement of each image in the data set with respect to the
constant RI is an embarrassingly parallel operation. The data set can thus be
split across multiple processing units for a significant speed-up.

###  3.2 Real-time Analysis

Real-time analysis is inherently a sequential process that requires
continuously updating the RI as each new image is processed. The computational
expense is large, as both the displacement and the RI must be updated at least
at the camera’s frame rate. As each new image is processed, the RI is updated
using a low-pass filter,

$\mathcal{F}[E_{i}]=a_{0}\mathcal{F}\bigl{[}I_{i}(\vec{r}-\vec{r}_{0})\bigr{]}+b_{1}\mathcal{F}\bigl{[}E_{i-1}\bigr{]}\;,$
(8)

where $\mathcal{F}[E_{0}]=\mathcal{F}[I_{0}]$, $a_{0}=1-e^{-T/\tau}$,
$b_{1}=e^{-T/\tau}$, $T$ is the inverse of the camera’s frame rate, and $\tau$
is the filter time constant (typically $\tau\gg T$).

This method may experience drift over time as a result of the gradual
accumulation of numerical errors, which can displace the particle in the RI
from its initial position. Similarly, drift in the particle’s actual position
can occur over time due to changes in its equilibrium position. Although these
drifts usually have no significant consequences, they can be inconvenient, and
it is often desirable to eliminate them from the detected motion, retaining
only the faster particle oscillation.

A natural solution is to subtract off the (possibly drifting) equilibrium
position, as determined by a low-pass filter. Given the $n$-th raw
displacement $\vec{r}_{n}$, we calculate the filtered displacement value
$\vec{r}_{n}^{\prime}=\vec{r}_{n}-\vec{R}_{n}$, where $\vec{R}_{n}$ is the
subtracted offset calculated by passing the raw data through a low pass
filter,

$\vec{R}_{n}=a_{0}\vec{r}_{n}+b_{1}\vec{R}_{n-1}\;,$ (9)

and $\vec{R}_{0}=a_{0}\vec{r}_{0}$.

If there is concern that Eq. (9) may be too aggressive at correcting drift and
could inadvertently remove an important signal, the shift in the offset can be
limited to a fraction $\alpha$ of the resolution of the CC analysis (the
inverse of the upsampling factor $u^{-1}$). In this case, the drift correction
effectively biases the round-off error in the CC towards a zero average
displacement, and the change in the correction offset can be written as

$\begin{split}d\vec{R}_{n}&=\bigl{(}a_{0}\vec{r}_{n}+b_{1}\vec{R}_{n-1}\bigr{)}-\vec{R}_{n-1}\\\
&=a_{0}\vec{r}_{n}+\bigl{(}b_{1}-1\bigr{)}\vec{R}_{n-1}\;,\end{split}$ (10)

and

$\vec{R}_{n}=\begin{cases}\vec{R}_{n-1}+(\alpha/u)d\hat{R}_{n}&,\ {\rm if}\
\bigl{|}d\vec{R}_{n}\bigr{|}>\alpha/u\;,\\\ \vec{R}_{n-1}+d\vec{R}_{n}&,\ {\rm
otherwise}\;,\end{cases}$ (11)

where $d\hat{R}_{n}=d\vec{R}_{n}/\bigl{|}d\vec{R}_{n}\bigr{|}$. A typical
choice of parameter is $\alpha=0.5$, ensuring that drift corrections are
smaller than the expected displacement resolution.

##  4 Results

We perform tests with synthetic data, where each synthetic image comprises a
Gaussian function integrated over each pixel, representing either a bright or
dark particle on a dark or bright background, with Poisson shot noise added
(see Fig. 1). To ensure that these synthetic images are experimentally
relevant, we match their maximum pixel intensity to the saturation capacity of
the CMOS sensor used in our experiments (10700 electrons for a Sony IMX250).
When adding dark noise to the images, we also match it to the dark noise of
the sensor (2.4 electrons). For each method and choice of parameters, we
analyze data sets consisting of 1000 images, whose displacements are randomly
distributed, but exactly known.

Figure 1: Gaussian particle tracking study: comparison of various image
analysis methods, both in the presence and absence of dark noise. Legend
indicates the particular method used: center of mass (CM), CM subtracting the
mean background (CM $-$ BG), uniform cross-correlation (CC-U), shot-noise-
weighted cross-correlation (CC-SN), and maximum likelihood estimation (MLE).
The Gaussian width as a fraction of the image width ($w_{0}/N$) is used as a
comparison metric. Insets show samples of the analyzed images: the top half
presents the original, and the bottom half adjusts the pixel values (through
rescaling, offsetting, and saturating) to highlight the presence or absence of
shot noise.

For the bright spot on a dark background we consider both the case without and
with dark noise, where a constant value is added to the entire image before
the Poisson noise is introduced. This means that subtracting the mean
background value from the image with dark noise does not yield, in general,
the image without dark noise. For the dark spot on the bright background we
only consider the case with added dark noise.

We compare the uniformly-weighted (CC-U) and shot-noise-weighted (CC-SN)
cross-correlation methods described in Sec. 2 with the standard center of mass
(CM) and maximum likelihood estimation (MLE) methods, using the real-time
approach described in Sec. 3. When using the MLE method, we provide the exact
displacement as the initial value for the iterative optimization algorithm. By
doing so, the MLE method produces results that represent the best achievable
outcomes from numerical methods, offering a benchmark against which other
methods can be compared.

Figure 1(a) illustrates the performance of each method in recovering
displacements for a bright spot in the absence of dark noise. In this case,
all methods exhibit roughly the same performance as the MLE method,
effectively reaching the shot noise limit for the standard deviation of the
position error, $w_{0}/(2\sqrt{N_{\gamma}})$, where $N_{\gamma}$ is the total
number of photons detected in the Gaussian [13]. This is valid as long as the
particle size (proportional to the Gaussian width $w_{0}$) does not constitute
a large fraction of the image. In this idealized context, it is evident that
the CM algorithm is the preferred choice, given its simplicity and low
computational cost.

Figure 1(b) shows the performance of each method in recovering displacements
for the more realistic case of a bright spot with dark noise. In this case,
the CM algorithm is used to analyze both the original images and the images
with the exact mean background subtracted (CM $-$ BG). While subtracting the
mean background substantially enhances the CM method’s performance, it still
falls short of reaching the shot noise limit, except for a small range of
$w_{0}$, a feature which is used in [5]. On the other hand, both the CC-U and
CC-SN methods yield results comparable to the ones obtained using the MLE
method, and are therefore preferred due to their significantly lower
computational cost. As in the previous case, the error for the CC-U method is
slightly larger than the ones for the CC-SN and MLE methods.

Figure 1(c) shows how each method performs in recovering displacements for a
dark spot in the presence of dark noise. When using the CM method, we subtract
the image from the exact mean background value to approximate a bright spot on
a dark background. Nevertheless, the displacement errors for the CM method are
still about an order of magnitude above the shot noise limit. While the CC-U
method exhibits larger errors than the MLE method, they are still notably
smaller than those associated with the CM method. The errors for the CC-SN
method remain comparable to the ones for the MLE method.

Figure 2: Experimental data: backlit levitated borosilicate glass microsphere.
Artifacts in the image background are from dust on the optics and window of
the vacuum chamber.

###  4.1 Experimental data

We show the CC-SN method can be used for real-time analysis by tracking a
borosilicate glass microsphere levitated in a magneto-gravitatational trap in
high vacuum. The microsphere is back-illuminated with collimated light from a
pulsed $660\text{\,}\mathrm{nm}$ LED and imaged on a CMOS camera (Basler
acA2440-$75\text{\,}\mathrm{\SIUnitSymbolMicro m}$), giving the particle the
appearance of a dark disk on a bright background with a bright spot in the
center (see Fig. 2; see also [12, 3]).

Additional force to suspend the approximately
$68\text{\,}\mathrm{\SIUnitSymbolMicro m}$ diameter microsphere is provided by
charge on the microsphere in an electric field produced by applying a voltage
across the pole pieces [3]. The recorded data consists of $256\times 128$
pixel images of the microsphere, corresponding to a field of view
approximately $310\times$150\text{\,}\mathrm{\SIUnitSymbolMicro m}$$, with 12
bits per pixel for ten minutes at a rate of $470\text{\,}\mathrm{Hz}$, which
is near the maximum frame rate supported by the camera’s USB3 interface.

Because the microsphere image displays complex structure with non-Gaussian
features such as diffraction rings on the edges and a bright spot at the
center, the CM algorithm is not suitable for tracking the particle. We instead
use the CC-SN algorithm to measure the displacement of the particle in real-
time using a GeForce RTX 3080 Ti GPU, while simultaneously recording the
images. For comparison, we perform an offline analysis using the CC-SN
algorithm and our iterative method described in Sec. 3 up to the fifth RI. In
both cases the upsampling factor is set to $256$.

We find excellent agreement between the real-time and offline methods, with
differences less than 1.5 pixels for the entire data set, as shown in Fig. 3.
We note that the real-time analysis uses the slow drift correction algorithm
described in Sec. 3, but the iterative method does not. The predictable
consequence is a systematic drift in both the horizontal and vertical
displacements, observed in the bottom panels of Fig. 3.

Figure 3: Comparison between real-time and offline algorithms at the fifth
reference image for experimental data. The calibration is
$1.2\text{\,}\mathrm{\SIUnitSymbolMicro m}$/pixel. Top: Displacements $r_{0}$.
The apparent beating in the detected vertical motion is due to sensitivity to
a combination of multiple degrees of freedom separated in frequency by
${\sim}$$3\text{\,}\mathrm{Hz}$. Bottom: residual (in pixels) of the
displacements from both methods.

##  5 Conclusions

In this paper we introduced a shot-noise weighted cross-correlation (CC-SN)
method, which can be used to determine image displacements approaching the
shot noise limit. The method was tested on several different images, including
a bright spot on a dark background, with and without dark noise, a dark spot
on a bright background with shot noise, and on complex experimental data. We
found that the CC-SN method is superior to the uniformly weighted cross-
correlation (CC-U) method, with errors comparable with those obtained using
the maximum likelihood (MLE) method, while being significantly cheaper than
MLE from a computation standpoint.

CC-SN is particularly useful for feedback control systems, as it can be used
for real-time image analysis. As a proof of principle, we used the CC-SN
method for real-time tracking of the position of a microsphere levitated in a
magneto-gravitational trap in high vacuum. Ultimately, we anticipate using the
method for feedback cooling of the motion of a particle to the shot-noise
limit as well as state preparation and tracking for a new measurement of the
Newtonian constant of Gravitation using a levitated particle.

Using this method for real-time analysis requires an implementation of the
algorithms discussed in this paper that is capable of processing images at a
rate at least the camera’s frame rate. To this end, we provide an open-source
toolkit for real-time image analysis called RETINAS [14]. The toolkit contains
the methods described in Sec. 2 implemented in CUDA, C, and Python, all of
which can be accessed via a user-friendly Python interface. Ref. [14] contains
both the code and extensive documentation on how to use the code and reproduce
all of the results shown in this paper.

The real-time CC-SN method described in this paper represents a significant
advancement in real-time image analysis, bridging a critical gap in current
particle tracking algorithms. Its application and further development offer a
way to streamline real-time image analysis across various disciplines,
offering a technique that is not only effective but also accessible and
adaptable.

Funding This material is based upon work supported by the National Science
Foundation under Grant No. 1806596, 1912083, 1950282, 2011783, and 2227079.

Acknowledgement BD gratefully acknowledges a Block Gift from the Coherent /
II-VI Foundation. This research made use of the resources of the High
Performance Computing Center at Idaho National Laboratory, which is supported
by the Office of Nuclear Energy of the U.S. Department of Energy and the
Nuclear Science User Facilities under Contract No. DE-AC07-05ID14517.

Disclosure The authors occasionally used ChatGPT for language and readability
improvements, and assume full responsibility for the contents of the
publication.

Data availability Data underlying the results presented in this paper are
available in Ref. [15].

## References

  * [1] A. R. Wade and F. W. Fitzke, Optics Express 3, 190 (1998).
  * [2] B. Davison _et al._ , Journal of Geophysical Research: Earth Surface 125, e2019JF005492 (2020).
  * [3] C. W. Lewandowski _et al._ , Phys. Rev. Applied 15, 014050 (2021).
  * [4] Y. Minowa _et al._ , arXiv preprint arXiv:2204.05506 (2022).
  * [5] A. J. Berglund _et al._ , Optics express 16, 14064 (2008).
  * [6] M. Guizar-Sicairos _et al._ , Optics letters 33, 156 (2008).
  * [7] S. van der Walt _et al._ , PeerJ 2, e453 (2014).
  * [8] A. Paintdakhi _et al._ , Molecular microbiology 99, 767 (2016).
  * [9] J. M. Graving _et al._ , Elife 8, e47994 (2019).
  * [10] B. H. Savitzky _et al._ , Microscopy and Microanalysis 27, 712 (2021).
  * [11] Argunşah _et al._ , Scientific Reports 12, 12405 (2022).
  * [12] B. R. Slezak _et al._ , New Journal of Physics 20, 063028 (2018).
  * [13] M. T. Hsu _et al._ , Journal of Optics B: Quantum and Semiclassical Optics 6, 495 (2004).
  * [14] L. R. Werneck, RETINAS GitHub (Visited on: October 12, 2023).
  * [15] L. R. Werneck _et al._ , Zenodo, doi:10.5281/zenodo.8336549 (2023).

references