# Estimating gravimetric effects for ordinary bodies

Pier Franco Nali1
1Independent scholar, Cagliari, Italy.
<EMAIL_ADDRESS>

(Submitted on dd-MMM-yyyy)

###### Abstract

The main topic of this article is a discussion about the best way to show
students that the proportionality of mass and weight, strictly true for point-
like particles, is an excellent approximation for objects of “normal” size.
The usual way of addressing this issue, although very simple, is not entirely
satisfactory. Our approach considers first and second order, coordinate
dependent, gravimetric effects, connected to the internal geometry of objects;
these effects, extremely small, are estimated through examples.

##  1 Introduction

Direct proportionality of mass and weight is a well-established principle,
proven as an experimental fact for all bodies in the same place. However,
apart from the special case of uniform gravitational field, this principle is
only valid locally, that is for point particles. When both the variability of
the gravitational field and the bodies’ internal structure cannot be ignored,
the point-particle approximation fails and the proportionality of mass and
weight cannot be regarded as strictly exact. Anyway, discrepancies involved
are generally tiny and can safely be overlooked in most situations we commonly
experience.

It should be deemed that teachers usually show students that for common bodies
the acceleration of gravity does not vary appreciably (or, otherwise said, the
gravitational field is uniform) within the size of the object, which can be
done very simply. At least at the college level, but also at the high-school
level, after having presented relative motions and introduced the “apparent
forces” in non-inertial reference frames - or “inertial forces” as we want to
call them - teachers specify that the weight force on the earth’s surface is
the resultant of the (true) gravitational force and of the apparent (inertia)
forces, in particular the centrifugal force if the body is stationary.

Limiting the problem to the gravitational component only (let’s say it
$\mathbf{G}$) teachers follow the usual simple path of differentiation $\Delta
G/G=\Delta R^{-2}/R^{-2}=-2\Delta R/R=-2h/R$, where $R$ is the radius of the
earth and $h$ is the height of the body (or, if students do not know the
differentials, simply calculate $G(R+h)/G(R)=R^{2}/(R+h)^{2}$). Since the
difference in the first order is already very small, $2h/R$ typically being of
the order of $10^{-8}$, that of the second order is obviously negligible for
any practical effect.

But we immediately realize that this method is flawed as the inverse
proportion of the gravitational force to the square of the distance is
strictly valid for point particles (and in the other special case of
centrally-symmetric homogeneous bodies). In this way we implicitly assume the
conclusion, as affirming uniformity of the gravitational field within point
objects is tautological: our argument contains a circularity and runs into a
logical fallacy (petitio principii). Furthermore, properly speaking, $h$ does
not correlate to the object size: it merely represents a (small) displacement
of the body (or rather, its center of mass) from the surface of the earth.

Thus the argument is ineffective for extended bodies of arbitrary shape as it
does not properly capture the effect on weight force (in either a rotating or
a non-rotating frame) of the variation of $\mathbf{g}$ within the object size;
it would be preferable to find a different approach, allowing us to address
the topic in broader generality and rigour.

##  2 Background

According to the “Declaration on the unit of mass and on the definition of
weight; conventional value of $g_{n}$”: ¡¡The word “weight” denotes a quantity
of the same nature as a “force”: the weight of a body is the product of its
mass and the acceleration due to gravity; in particular, the standard weight
of a body is the product of its mass and the standard acceleration due to
gravity.¿¿[1] And, contextually: ¡¡The kilogram is the unit of mass; it is
equal to the mass of the international prototype of the kilogram.¿¿111The
definition of the kilogram in terms of the international prototype is obsolete
and no longer in force since 20 May 2019; it has been redefined in terms of
the Planck constant. The value adopted for the standard acceleration due to
gravity (on earth) is
$g_{n}=$980.665\text{\,}\mathrm{cm}\text{\,}{\mathrm{s}}^{-2}$$. Thus, for the
weight force:

$\mathbf{w}=m\mathbf{g},$ (1)

where $m$ is the mass of the object and $\mathbf{g}$ is the acceleration
vector due to gravity (We denote vectors, like $\mathbf{g}$, $\mathbf{w}$,
$\mathbf{G}$, as bold letters and represent their magnitudes, like $g$, $w$,
$G$, as italic letters).

More generally, weight means the gravitational force (or this plus the
centrifugal force) on a small mass compared to that of the source (e.g. “the
weight of astronauts on the Moon”).

It is important to notice that “weight” and “gravitational force” are the same
force but the use of either term is contextual and it is good practice to
adhere to conventions on their use to avoid ambiguity. Calling the
gravitational force on a celestial body “weight” creates confusion, and
contradicts the convention that reserves this word for practical use (on
weight vs gravitational force see e.g.[2]).

Hence the use of “weight” should be reserved for the force experienced by an
object with mass in a gravitational field (e.g. light is bent by gravitational
fields, although it has not weight, because it is massless).

We are now ready to introduce the alternate differentiation pathway
$\Delta\mathbf{g}/g=1/g\left(\Delta\mathbf{g}/\Delta R\right)\Delta
R=R/g\left(\Delta\mathbf{g}/\Delta R\right)h/R$, where $h$ is now the linear
size of the object and $\mathbf{g}(R)$ (gravitational + centrifugal) does not
have a form given a priori. In this way we see that the less restrictive
condition that the gradient $\Delta\mathbf{g}/\Delta R$ is the same order of
$\mathbf{g}/R$ would fulfill the requirement of uniformity of the
gravitational field within the size of the object.

In a more formal way, if we think of the body made up of particles of masses
$m_{i}$, being $m=\sum_{i}m_{i}$ the mass of the entire body and
$\rho(\mathbf{r})=\sum_{i}m_{i}\delta(\mathbf{r}-\mathbf{r}_{i})$ its density,
where $\delta$ is the Dirac delta function, the weight force is the
generalization of eq. (1) by integration over the whole space:

$\mathbf{w}=\iiint\rho(\mathbf{r})\mathbf{g}(\mathbf{r})\,\differential\tau=\sum_{i}m_{i}\mathbf{g}(\mathbf{r}_{i}).$
(2)

For an uniform field ($\mathbf{g}(\mathbf{r})=g\hat{\mathbf{z}}$) the equation
(2) reduces to (1) and the internal size and geometry of the body are
irrelevant.

But the uniformity condition for $\mathbf{g}$ only holds approximately near
the earth’s surface. The earth’s gravitational field is not uniform even on a
small scale; modern gravimeters allow us to appreciate $g$ with eight or nine
significant digits (some $\mu$gals); this is how to say that variations in the
earth’s gravity between points even a few centimeters apart are detectable
instrumentally.

In this article we consider very small (some ppb) coordinate dependent
effects. At these scales there are several others effects, both instrumental
and environmental, which are not so weak. For example, $\mathbf{g}$ has a
dependence on time: the effect of terrestrial tides alone is two orders of
magnitude greater (a few hundreds $\mu$gals), to which are to be added the
effects of tides in the oceans, the hydrological and barometric components,
also variable over time, and so on; also the motion (if any) of the measuring
instrument has to be taken into account. Moreover, although independent of the
mass $m$ of the body under consideration, $\mathbf{g}$ generally depends on
“other” masses; it will be assumed that these external masses vary very
slowly.

Obviously we should not forget to mention the major non-gravitational
contribution to weight, that of the centrifugal force due to diurnal rotation
(which is a component of $\mathbf{g}$); the centrifugal force is some part per
thousand of the gravitational force and the dependence of the two forces on
the distance, from the center or from the axis, is different; the effect on
the weight of Archimedes’ thrust in the air is also significant. In the
following we shall leave all these effects just mentioned aside from the
present study and focus our attention on coordinate related ones, dependent on
size and geometric configuration of objects.

It is easy to see that for bodies for which experiments can be established,
such as for bodies near the earth’s surface, these coordinate-dependent
effects are far too small for the standard resolution of dynamometers and
scales. To this end we take Taylor series expansions truncated to the second-
order of the functions $\mathbf{g}(\mathbf{r}_{i})$ in eq. (2) centered about
$\mathbf{R}_{cm}=\frac{1}{m}\sum_{i}m_{i}\mathbf{r}_{i}$, that is the radius
vector conducted from the origin of the coordinates to the body center of
mass:222Under suitable analyticity conditions for the functions
$\mathbf{g}(\mathbf{r}_{i})$.

By fixing the origin of the coordinates in the center of the earth, the
equation (LABEL:eq3) is quite exact for all applications in which we need to
evaluate the weight force on bodies located on or near the earth’s surface. In
fact, under these conditions $|\mathbf{r}_{i}-\mathbf{R}_{cm}|\ll R$ with
$|\mathbf{R}_{cm}|\cong R$ (assuming
$R\simeq$6.371\text{\times}{10}^{6}\text{\,}\mathrm{m}$$ for the earth’s mean
radius) and the third and higher-order terms can be overlooked.

##  3 The first order effect. Implications for precision mass measurements

For a (small) displacement of a body from $P$ to a near point $P^{\prime}$ we
can express the variation of $\mathbf{g}$ within the body by means of the eq.
(LABEL:eq3) as:

$\mathbf{g}(\mathbf{r^{\prime}}_{i})=\mathbf{g(P^{\prime})}+\left(\mathbf{r^{\prime}}_{i}-\mathbf{P^{\prime}}\right)\dotproduct\gradient\mathbf{g}\bigr{\rvert}_{\mathbf{P^{\prime}}}$
(4)

(overlooking the small second order term). For a rigid body holds the distance
preserving condition
$\norm{\mathbf{r^{\prime}}_{i}-\mathbf{P^{\prime}}}\equiv\norm{\mathbf{r}_{i}-\mathbf{P}}$.
For the sake of simplicity we assume the special condition of a purely
translational displacement (preserving distance, angle, sense, and
orientation) such that
$\mathbf{r^{\prime}}_{i}-\mathbf{P^{\prime}}\equiv\mathbf{r}_{i}-\mathbf{P}$,
so that, if we take $P$ as the center of mass, the second term of the right-
hand side vanishes identically by introducing eq. (4) in eq. (2), and eq. (4)
reduces to

$\mathbf{g(P^{\prime})}=\mathbf{g}_{cm}+\left(\mathbf{P^{\prime}}-\mathbf{R}_{cm}\right)\dotproduct\gradient\mathbf{g}\bigr{\rvert}_{\mathbf{R}_{cm}}.$
(5)

In turn, eq. (2) reduces to

$\mathbf{w^{\prime}}=\mathbf{w}_{cm}+\mathbf{w}_{-1},$ (6)

where

$\mathbf{w}_{-1}=m\left(\mathbf{P^{\prime}}-\mathbf{R}_{cm}\right)\dotproduct\gradient\mathbf{g}\bigr{\rvert}_{\mathbf{R}_{cm}}.$
(7)

Applying the gradient criterion $\gradient\mathbf{g}\sim\mathbf{g}/R$ we
easily obtain $w_{-1}/w_{cm}\sim h/R$, so, for a body similar in size and mass
to the obsolete kilogram prototype ($h\sim$ a few centimeters) $w_{-1}\sim
10^{-7}-10^{-8}$ N.

In recent decades, in view of the redefinition of the SI units, in particular
the kilogram, the goal set by the CGPM was to achieve accuracy of the order of
$10^{-8}$. Laboratories such as the National Physical Laboratory in the UK and
numerous other metrology labs around the world have worked for years to
achieve the required accuracy with Kibble’s balances. At the meeting of 17th
May 2019 of the CCM, I. A. Robinson (NPL) stated: ¡¡Whilst, at present, it is
theoretically possible to measure the principal quantities to around 2–3 parts
in $\mathit{10^{9}}$ a number of other effects in the apparatus must be taken
into account.¿¿[3] This is a number of practical reasons, which limit
accuracy, such as alignments, vibrations, etc. The NPL also plans to develop
simpler Kibble’s balances, affordable and operable in laboratories not as
highly–specialized as NPL, capable of $10^{-8}$ accuracy. At this accuracy
level, a number of systematic effects has to be taken into account, including
gravimetric contributions (see e.g. [4]). An historical account of the
development of these sensitive balances in the context of the proposed reform
of the SI is outlined in [5].

Apparently, gravimetric effects such as those we are talking about were first
considered in the early 1970s in connection with the development at the
National Bureau of Standards of the “One Kilogram Balance” NBS No. 2, whose
standard deviation was approx. $4\text{\,}\mathrm{\SIUnitSymbolMicro g}$.[6]

Such kind of balances, used for comparing masses, compare the attractive
gravitational forces between weights (or loads) and the earth. It is assumed
(often implicitly) that these forces are exactly proportional to the masses of
the loads (in vacuum) and do not vary during the measurement. The force on a
standard weight used for the comparison of masses depends on the distance from
the center of the earth to the center of gravity of the weight.333We
anticipate here the notion of center of gravity that we will resume later. To
practical effects of the discussion carried out in this section we can
consider the center of gravity coincident with the center of mass, although
the two concepts, in general, are to be kept distinct. A second weight, of a
different configuration, may have its center of gravity at a different
distance from its base and thus the distance of the weight’s center of gravity
from the center of the earth will be different when the weight is placed on
the weighing pan (which operates with the bases of the weights to be compared
virtually on the same level). In this way, the constant of proportionality
between the gravitational forces and the masses of the weights on the pan will
be slightly altered, leading to a systematic error in the results of the
comparisons between the masses, the so-called “gravitational configuration
effect” introduced by Almer and Swift.[7]

If we consider a reference weight $w_{r}=m_{r}g(R)$ and a second equal weight
$w_{x}$, whose centers of gravity are spaced by a distance $\Delta h=d$ above
their bases, then, from eqs. (6) and (7):

$\displaystyle w_{x}=w_{x,cm}+w_{x,-1},$ $\displaystyle w_{x,cm}=m_{x}g(R),$
$\displaystyle
w_{x,-1}=w_{x,cm}\cdot\frac{1}{g}\partialderivative{g}{h}\left(\Delta
h\right),$

or, in the approximation of the gravitational component alone,

$w_{x,-1}=-w_{x,cm}\cdot\frac{2d}{R}.$ (8)

And, having imposed $w_{x}=w_{r}$,

$m_{x}=m_{r}+m_{-1}=m_{r}+\frac{2dm_{r}}{R}.$ (9)

The term $m_{-1}=2dm_{r}/R$ in eq. (9) is the (first order) corrective term
that must be applied to the mass of the second weight to take into account the
difference in the force of gravity on the weights placed on the weighing pan
of the balance whose centers of gravity are at different distances from their
bases. The corrective term can be evaluated independently of the equation (8),
valid in the approximation of the gravitational component alone, by directly
measuring the acceleration of the free fall $g$ and the gradient of the
gravitational field $\partialderivative*{g}{h}$ in the place in which the mass
calibration takes place.

In the case $d=$+1\text{\,}\mathrm{cm}$$ the correction for the comparison of
nominal weights of $1\text{\,}\mathrm{kg}$ calculated using the equation (9)
is approx. $+3\text{\,}\mathrm{\SIUnitSymbolMicro g}$. The old Pt-Ir kilogram
prototype (density $21.55\text{\,}\mathrm{kg}\text{/}{\mathrm{dm}}^{3}$) is a
right circular cylinder with a volume of approx.
$46.5\text{\,}{\mathrm{cm}}^{3}$ and approx. the same height
($39\text{\,}\mathrm{mm}$) as the diameter. Stainless steel samples (density
$8.00\text{\,}\mathrm{kg}\text{/}{\mathrm{dm}}^{3}$), having volume
($125\text{\,}{\mathrm{cm}}^{3}$) respecting the same proportions, have a
height of $54.2\text{\,}\mathrm{mm}$. The resulting distance of the samples’
centers of mass (/gravity) from their base is higher than that of the
prototype’s center of mass from its base by an amount of
$7.6\text{\,}\mathrm{mm}$, which leads to a correction of
$+2.4\text{\,}\mathrm{\SIUnitSymbolMicro g}$. For comparison, as Almer and
Swift stated: ¡¡Currently, mass comparisons at the 1-kg level can be carried
out with standard deviations as small as 1.5 parts in 109.¿¿[7]

This correction is far from being the most significant; the largest volume
($\approx$80\text{\,}{\mathrm{cm}}^{3}$$) of the stainless steel 1-kg samples
results in a correction for the aerostatic thrust of approx.
$+94\text{\,}\mathrm{mg}$ (assuming an air density of
$1.2\text{\,}\mathrm{kg}\text{/}{\mathrm{m}}^{3}$), that is about 40,000 times
the gravitational effect.[8] Nonetheless the gravitational correction becomes
significant for high precision mass measurements. In fact, accuracy is limited
not only by the achievable precision and uncertainty associated with the value
of the sample, but also by systematic errors. It can be said that the accuracy
of the results of the measurements is achieved only after all the relevant
systematic errors have been identified and evaluated. This implies that in the
design of an experiment all factors, even those that at first appear small,
must be estimated to establish their potential importance as systematic
factors affecting the measured results.

##  4 The second order effect

By introducing eq. (LABEL:eq3) into eq. (2) and noticing that the first order
term vanishes identically for the choice of $\mathbf{R}_{cm}$ , eq. (2)
reduces to

$\displaystyle\mathbf{w}=\mathbf{w}_{cm}+\mathbf{w}_{-2},$

where

$\displaystyle\mathbf{w}_{cm}=m\mathbf{g}_{cm}$

and

$\displaystyle\mathbf{w}_{-2}=\frac{1}{2}\sum_{i}m_{i}\left(\left(\mathbf{r}_{i}-\mathbf{R}_{cm}\right)\dotproduct\gradient\right)^{2}\mathbf{g}\bigr{\rvert}_{\mathbf{R}_{cm}}$

(or, in the continuous limit)

$\displaystyle=\frac{1}{2}\iiint\rho(\mathbf{r})\left(\left(\mathbf{r}-\mathbf{R_{cm}}\right)\dotproduct\gradient\right)^{2}\mathbf{g}\bigr{\rvert}_{\mathbf{R_{cm}}}\,\differential\tau.$

The term $\mathbf{w}_{cm}$ is the weight force acting on the material point to
which the body is reduced, having the mass of the body and located in its
center of mass.

The term $\mathbf{w}_{-2}$ is a second order gravimetric correction that takes
into account the effect of the internal geometry of the body, estimated as
follows:

$\displaystyle
w_{-2}=w_{cm}\cdot\frac{1}{2}\cdot\frac{1}{g}\partialderivative[2]{g}{z}\left(\Delta
z\right)^{2},$

or, in the purely gravitational component approximation,

$\displaystyle w_{-2}=w_{cm}\cdot\frac{3d^{2}}{R^{2}}$

(where $\Delta z=d$ is the linear size of the object). It represents the
difference due to mass distribution around the center of mass compared to the
situation in which all the mass is thought to be concentrated in one point.
More formally, it can be shown (see e.g.[9]) that the mass distribution
intervenes to second order through the inertia tensor of the body. For a right
circular cylinder of mass $1\text{\,}\mathrm{kg}$ a few centimeters high, like
a copy of the old Pt-Ir kilogram prototype, the order of magnitude of the
$w_{-2}$ term is $\sim 10^{-16}-10^{-17}$ N, the same of the weight of the
equivalent mass of 1 joule, just 1/10 of that of the mass of an Escherichia
coli bacterium and one hundred thousand times smaller than that of the mass of
a human cell.

Although fully negligible for bodies of ordinary mass near the surface of the
earth, similar but a bit more significant effects occur in various kinds of
problems, often faced with methods borrowed from celestial mechanics; in these
situations, all the possible contributions must be carefully evaluated both in
theoretical analyses and in the design of the experiments. A typical example
are tidal phenomena, whose effects depend on the gradient of the gravitational
field, rather than on intensity, and the variations of the gravitational force
from one part of the object to the other must be considered. Meanwhile, there
is no doubt that in these situations the bodies cannot be thought of as
material points; Newton had already noticed that the exact results obtained
for point-like particles are only approximate in presence of gravitational
force between extended bodies attracting at short distances. In celestial
mechanics it is usually satisfactory to stop calculations at the second order
of approximation.

Moreover, as the size of the objects under consideration are on a planetary or
sub-planetary scale, i.e. a significant fraction (say, from a few thousandths
to a few hundredths) of the earth’s radius (think, for example, of lithosphere
segments of which we want to study the isostatic conditions), or when the
bodies are very close to an attracting center (a situation encountered in
geophysical and astrophysical contexts), also the assumptions under which the
equation (LABEL:eq3) holds can fail and additional contributions should be
considered.

In addition, sometimes it is not even possible to set up experiments or carry
out direct measurements; when this occurs, the evaluation of gravitational
forces needs ad hoc modeling of objects, which may require, for example, the
computation of quadruple or sextuple integrals and numerical integration (see,
for example,[10]).

##  5 The elusive center of gravity. Near–uniform field

The slight variation of the gravitational field within the size of earthly
objects brings us to the interesting questions of the parallel field and the
center of gravity.

The earth’s gravitational field can be locally modeled by a field consisting
of parallel vectors of (slightly) non-uniform intensity. This picture is
useful because it allows us to introduce the “scalar weight” $w$ in a coherent
way,444“scalar” here does not mean invariant under rotation; here we intend
1-dimensional 1-component scalar field. providing a tool to face and clarify
the problem of determining a unique point (if any) where you can think applied
the total weight force acting on all the particles of the body, i.e. its
center of gravity.555The center of gravity is susceptible to other
definitions, which we will not deal with here. A definition different from
that of the weighted average can be given, for example, in the case of the
spherically symmetric field. A real gravitational field cannot be both
parallel and non-uniform at the same time. It is convenient to examine the
case of the near-uniform field, which, in addition to being simplistic,
reproduces the gravitational field near the earth’s surface with an excellent
degree of approximation. Furthermore, with this choice, the problem can be
dealt with in one dimension. For the usual central field

$\displaystyle\mathbf{g(r)}=-k\,\frac{\mathbf{r}}{\norm{r}^{3}}$

($k=GM$ for the earth’s gravitational field) $\divergence\mathbf{g}=0$
everywhere. In the near-uniform model we consider a small cylindrical region
where there is a field of vectors parallel to $\hat{\mathbf{z}}$, having non-
uniform modulus, so defined:

$\mathbf{g(r)}=g(z)\hat{\mathbf{z}}=-kz^{-2}\hat{\mathbf{z}},$ (10)

with $z\gtrsim R$.

The equation (10) does not represent a real Newtonian gravitational field as
$\mathbf{g}$ does not have zero divergence. However, for $z$ large enough,
i.e. far from the center of the field (e.g. near the earth’s surface), the
divergence is small and the eq. (10) is a very good approximation, locally
(far from the center of the earth), of a gravitational field generated by a
spherically symmetric mass distribution.666We assume the simplified picture of
spherical earth, uniform density, not rotating; we abstract from all possible
disturbing factors (assuming absence of air, no influence of celestial bodies,
etc.). The ininfluence of the body under examination on the central
gravitational field is also assumed (external field approximation). In this
framework, the center of gravity of a body can be defined through the
“equipollent” moment condition (see [11], p. 18). The moment of a single force
on a particle is perpendicular to the force and the vector radius from the
coordinate origin to the position of the particle. In general, however, this
is not true for a system of forces; the total moment of a system of forces
around a point $O$ (the pole, which we will also assume as the origin of the
coordinates) is generally not perpendicular to the total force vector acting
on the system.

The moment $\mathbf{T}_{eq}$ of the system of forces equipollent to a single
weight force $\mathbf{w}$ acting on the body satisfies the vector equation

$\mathbf{T}_{eq}=\mathbf{R}_{cg}\crossproduct\mathbf{w},$ (11)

where $\mathbf{w}$ is the total weight force acting on the body, defined by
eq. (2) and $\mathbf{R}_{cg}$ is the radius vector joining the pole with the
point of application of this force, i.e. with the body center of gravity. The
total moment of the forces acting on the system is by definition
$\mathbf{T}=\sum_{i}\mathbf{r}_{i}\crossproduct\mathbf{w}_{i}$ , and the total
weight force $\mathbf{w}=\sum_{i}\mathbf{w}_{i}$ . Imposing the
perpendicularity condition to these two vectors is equivalent to making the
equation (11) valid for $\mathbf{T}$, that we rewrite as

$\sum_{i}\left(\mathbf{r}_{i}-\mathbf{R}_{cg}\right)\crossproduct\mathbf{w}_{i}=\mathbf{0}.$
(12)

The equation (12) (torque equation) does not have solution if $\mathbf{T}$ and
$\mathbf{w}$ are not orthogonal (and neither is zero) and in this case the
center of gravity vector $\mathbf{R}_{cg}\equiv(X,Y,Z)$ cannot be determined
by this method. We do not examine here the existence conditions of the
solutions of the torque equation, whose detailed discussion can be found, for
example, in [12]. Fortunately, in the special case of parallel field the
orthogonality condition is met.777Another case in which this condition is met
is that of a planar system of forces. If we choose the $z$–axis in the
direction of the field, then $\mathbf{w}_{i}=w_{i}\,\hat{\mathbf{z}}$ and eq.
(12) reduces to the linear system

$\left\\{\begin{array}[]{lrcrcr}\sum_{i}&(x_{i}&-&X)\quad w_{i}&=&0\quad,\\\
\sum_{i}&(y_{i}&-&Y)\quad w_{i}&=&0\quad.\end{array}\right.$

The moment of total weight force will have only the $x$ and $y$ components
different from zero, from which the $X$ and $Y$ components of the vector
$\mathbf{R}_{cg}$ can be calculated; these define the line of action of the
total weight force. There remains the $z$ component to be determined (the
torque equation for the $z$ component is a null identity). We observe,
however, that under the assumptions made the equation (12) can be rewritten as

$\left(\sum_{i}w_{i}\mathbf{r}_{i}-w\mathbf{R}_{cg}\right)\crossproduct\hat{\mathbf{z}}=\mathbf{0}.$
(13)

Then, as the pole $O$ can be chosen arbitrarily and $\hat{\mathbf{z}}$ is a
fixed vector, the equation (13) can be satisfied by choosing the vector
$\mathbf{R}_{cg}$ defined as (see [11], p. 48)

$\mathbf{R}_{cg}=\frac{1}{w}\sum_{i}w_{i}\mathbf{r}_{i}=\frac{1}{w}\sum_{i}m_{i}g\left(\mathbf{r}_{i}\right)\mathbf{r}_{i}$
(14)

or, in the continuous limit,

$\mathbf{R}_{cg}=\frac{1}{w}\iiint\rho(\mathbf{r})g\left(\mathbf{r}\right)\mathbf{r}\,\differential\tau,$
(15)

which constitute the definition of the center of gravity in the case of
parallel field. For a uniform field the equation (15) becomes

$\mathbf{R}_{cm}=\frac{1}{m}\iiint\rho(\mathbf{r})\mathbf{r}\,\differential\tau,$
(16)

that is $\mathbf{R}_{cg}$ coincides with the center of mass radius vector. In
the equations (15) and (16) it is implied that

$\displaystyle w$
$\displaystyle=\iiint\rho(\mathbf{r})g(\mathbf{r})\,\differential\tau,$
$\displaystyle m$ $\displaystyle=\iiint\rho(\mathbf{r})\,\differential\tau.$

With series expansions of $w$ and $g\left(\mathbf{r}\right)$ in eq. (15)
around the center of mass,888See previous note 2. using the equations
(LABEL:eq3) and (10), and truncating after the first order, we have (we omit
the detailed steps):

$\displaystyle\mathbf{R}_{cg}=\mathbf{R}_{cm}+\mathbf{R}_{-1}+\cdots,$

If as an example we consider a solid in the shape of a right cylinder or a
rectangle parallelepiped, very elongated with respect to its basis, resting on
the earth’s surface so as to approach the situation of a parallel and near-
uniform field, we reduce the problem to one dimension. If $h$ is the height of
the solid, the $z$–coordinate of its center of mass will be given by
$Z_{cm}=R+h/2$; we also express the variable of integration as a function of
the coordinate in the system of the center of mass $\zeta=z-Z_{cm}$; finally,
for simplicity, suppose the solid of uniform density $\rho$. Then we can write
the $Z$ coordinate of the center of gravity as

The term $Z_{-1}=-h^{2}/6Z_{cm}\cong-h^{2}/6R$ in the first order of
approximation represents the displacement of the center of gravity apart from
the center of mass. This is a tiny difference: in the case of Dubai’s Burj
Khalifa, currently the tallest building in the world
($h=$829.80\text{\,}\mathrm{m}$$), the center of gravity is only about
$2\text{\,}\mathrm{cm}$ below the center of mass! The center of gravity is a
specially elusive concept. It identifies a defined point, but, unlike the
center of mass, it does not have a definite position. Its position depends, in
general, on the relative positions of the body under consideration and the
attractive mass. As can be seen from the equation (LABEL:eq20), when the
distance $Z_{cm}\cong R$ of the body from the center of the earth increases,
the center of gravity approaches the center of mass. This feature makes it
difficult to work with the center of gravity and in practice this concept is
seldom used. The detailed treatment of this and other interesting problems
related to the center of gravity is beyond our scope; an introductory
discussion on these topics can be found on the Wikipedia page “Centers of
gravity in non-uniform
fields”999https://en.wikipedia.org/wiki/Centers_of_gravity_in_non-
uniform_fields and related talk,
101010https://en.wikipedia.org/wiki/Talk:Centers_of_gravity_in_non-
uniform_fields to which the interested reader is referred.

##  6 Conclusions

We have established that the proportionality of mass and weight for ordinary
bodies can be taken as an excellent approximation in all cases of practical
interest. However, it is advisable for students to always clarify the limits
of validity of this approximation, both in their theoretical meaning and for
the aspects related to the sensitivity of the experiments. For this purpose,
the gradient criterion $\Delta\mathbf{g}/\Delta R\sim\mathbf{g}/R$ is suitable
for exploring the variation of the gravitational force within the size of the
body. While it is easy to show that this gravimetric effect is negligible for
ordinary bodies, special caution should be observed when, in investigating
certain areas, you go beyond the validity range of the point particle
approximation. In geophysics, hydrostatics and astrophysics various situations
are encountered of strongly inhomogeneous gravitational field and the
gravitational effects connected to the internal geometry of the bodies cannot
be neglected. Such effects must be carefully considered; for example: in
celestial mechanics and astrodynamics, in the calculation of the short-
distance interaction of non-spherical shaped bodies (see, e.g., [13, 14, 15]);
in geophysics, in the calculation of the gravimetric field of a polyhedral
plate (see, e.g., [16, 17, 18]); in hydrostatics, in the computation of the
thrust, where the pressure gradient is replaced by the product of the density
of the fluid and the gravitational field (see, e.g., [19]). These problems are
addressed on a case-by-case basis and often require the development of
specific solutions.

## Acknowledgments

I thank Santo Armenìa for drawing my attention to this topic. I also thank the
reviewer for helpful comments/suggestions.

## References

  * [1] CGPM,  Proceedings of the 3rd CGPM (1901), 1901, p. 70.
  * [2] I. Galili,  Int. J. Sci. Educ. 23, 1073 (2001).
  * [3] I. A. Robinson,  17th meeting of the CCM (17th May 2019).
  * [4] I. A. Robinson and S. Schlamminger,  Metrologia 53, A46 (2016).
  * [5] R. P. Crease,  Phys. World 24 (03), 39 (2011).
  * [6] H. E. Almer,  J. Res. Natl. Bur. Stand. (U. S.) 76C (01–02), 1 (1972).
  * [7] H. E. Almer and H. F. Swift,  Rev. Sci. Instrum. 46, 1174 (1975).
  * [8] Z. J. Jabbour and S. L. Yaniv,  J. Res. Natl. Inst. Stand. Technol. 106, 25 (2001).
  * [9] D. Hestenes,  New Foundations for Classical Mechanics (Kluwer, Dordrecht, 2002), p. 520.
  * [10] J. Stirling,  New J. Phys. 19, 073032 (2017).
  * [11] M. F. Beatty,  Principles of Engineering Mechanics, Volume 2 (Springer, New York, 2016).
  * [12] K. R. Symon,  Mechanics, 3rd edition (Addison-Wesley, Reading, 1971).
  * [13] J. Ashenberg,  Celestial Mech Dyn Astr 99, 149 (2007).
  * [14] Y. Shi, Y. Wang and S. Xu,  Celestial Mech Dyn Astr 129, 307 (2017).
  * [15] X. Hou, D. J. Scheeres and X. Xin,  Celestial Mech Dyn Astr 127, 369 (2017).
  * [16] B. Banerjee and S. P. Das Gupta,  Geophysics 42, 1053 (1977).
  * [17] D. Nagy,  Geophysics 31, 362 (1966).
  * [18] R. Karcol and R. Pašteka,  Pure Appl. Geophys. 176, 257 (2019).
  * [19] F. M. S. Lima and F. F. Monteiro,  Rev. Bras. de Ensino de Fis. 35, 3701 (2013).


# Observation of Localization Using a Noisy Quantum Computer

Kazue Kudo<EMAIL_ADDRESS>Department of Computer Science, Ochanomizu
University, Tokyo 112-8610, Japan  Graduate School of Information Sciences,
Tohoku University, Sendai 980-8579, Japan

###### Abstract

Quantum dynamics in a strongly disordered quantum many-body system show
localization properties. The initial state memory is maintained owing to slow
relaxation when the system is in the localized regime. This work demonstrates
how localization can be observed using a noisy quantum computer by evaluating
the magnetization and twist overlap in a quantum spin chain after short-time
evolution. The quantities obtained from quantum-circuit simulation and real-
device computation show their apparent dependence on disorder strength,
although real-device computation suffers from noise-induced errors
significantly. Using the exact diagonalization of the Hamiltonian, we analyze
how noise-induced errors affect those quantities. The analysis also suggests
how the twist overlap can reflect the information on the eigenstates of the
Hamiltonian.

##  I Introduction

Localization in quantum systems has been extensively investigated across a
wide variety of contexts. One particular area of interest is many-body
localization (MBL), which is eigenstate localization in strongly disordered
interacting systems. MBL can be regarded as an extended concept of Anderson
localization. Despite decades of investigation, MBL is still attracting
considerable attention [1, 2, 3, 4]. Many numerical studies of the MBL
transition using various quantities derived from eigenstates [5, 6, 7, 8, 9,
10, 11, 12, 13, 14, 15, 16, 17] and quantum dynamics [18, 19, 20] have
demonstrated features of a continuous transition from the thermal to localized
phases. However, the exponents of finite-size scaling analysis in quenched
disordered systems violate the Harris criterion [21, 22], indicating that the
system sizes used in those numerical studies are too small to properly analyze
the MBL transition. In other words, the transition observed in the numerical
works of finite-size systems is not a true phase transition but rather a
crossover. Although the existence of the true MBL phase is controversial, the
MBL transition point is located at a much stronger disorder if it exists [23,
24, 25, 26]. In this paper, we consider a finite-size system to be in the
localized regime where it displays certain MBL properties.

Numerous experimental studies have shown the crossover from the thermal to
localized regimes [27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39].
Methods based on the memory effect are commonly used to detect MBL in
experiments, as slow relaxation in the localized regime results in the initial
state remaining almost unchanged for a long time. Experimental investigations
often focus on simple quantities such as imbalance and magnetization, which
can be easily obtained by qubit measurements. However, these quantities do not
provide much information on eigenstates of the system Hamiltonian. Although
tomographic techniques can provide information on eigenstates, they require
significant effort. In contrast, twist overlap, which can also be easily
obtained from qubit measurements, exhibits distinct behavior beyond the memory
effect and can reflect the information on eigenstates [40]: The twist overlap
obtained from quantum dynamics behaves similarly to that evaluated using
eigenstates.

The objective of this study is to investigate the observation of MBL using a
noisy quantum computer. The effect of noise varies depending on the type of
observable. In this work, we aim to provide theoretical and numerical
explanations of why noise-induced errors have different effects on different
observables. In particular, we focus on localization detection based on
quantum dynamics and demonstrate the disorder dependence of both the
magnetization and twist overlap. A comparison of results obtained from
quantum-circuit simulation and real-device computation provides insights into
noise-induced errors. Our findings reveal that the twist overlap is more
vulnerable to noise-induced errors than the magnetization. Furthermore, we
analyze the results, especially in the localized regime, using the exact
diagonalization of the Hamiltonian. The analysis suggests why the twist
overlap is more susceptible to noise and how it can reflect the information on
eigenstates.

In this study, we employ the method of localization detection using a gate-
model quantum computer, inspired by experiments of localization detection
using a quantum annealer [37]. When using the quantum annealer, it is
difficult to adjust interaction strength and local field separately. In
contrast, a gate-model quantum device allows the independent control of the
parameters for interaction strength and local field, enabling the separate
control of disorder strength. However, the quantum-circuit computation of time
evolution suffers from approximation error owing to Suzuki–Trotter
decomposition [41, 42]. Increasing the Trotter number reduces the
approximation error, but it also increases the number of quantum gates,
thereby increasing gate errors during real-device computation. Thus, the
balance between the approximation and gate errors is critical in a noisy
quantum device. To reduce noise-induced errors, we focus on short-time
behavior instead of long-time evolution. Although approximate quantum circuits
can further reduce noise-induced errors [43, 44, 39], we employ a
straightforward method based on the Suzuki–Trotter decomposition in this work
because results obtained using this simple method are preferable for error
analysis.

##  II  Model and Methods

###  II.1 Hamiltonian and Observables

The Hamiltonian of a one-dimensional transverse Ising model with local random
fields consists of the Ising Hamiltonian $H_{\rm I}$ and the transverse-field
Hamiltonian $H_{\rm TF}$: $H=H_{\rm I}+H_{\rm TF}$, where

$\displaystyle H_{\rm I}$
$\displaystyle=\sum_{j=1}^{L-1}J\sigma^{z}_{j}\sigma^{z}_{j+1}+\sum_{j=1}^{L}h_{j}\sigma^{z}_{j},$
(1) $\displaystyle H_{\rm TF}$
$\displaystyle=-\sum_{j=1}^{L}\Gamma\sigma^{x}_{j}.$ (2)

Here, $L$ is the system size, and $\sigma_{j}^{x}$ and $\sigma_{j}^{z}$
represent the Pauli operators of components $x$ and $z$, respectively, at site
$j$. The interaction is antiferromagnetic; hence, the parameter $J$ is
positive. The strength of the transverse field is denoted by $\Gamma$ and is
also taken as positive. Throughout this study, these parameters are fixed to
$J=\Gamma=1$. The local fields $h_{j}$ are uniformly distributed random
numbers in the interval $[-w,w]$, where $w$ denotes the disorder strength.

We observe the magnetization and twist overlap at the end of the time
evolution, namely, at the final time $t=T_{\rm fin}$. Both observables are
easily obtained by measuring each qubit in a quantum device. These observables
averaged over different disorder realizations serve as probes for localization
[40]. In this study, the initial state is given as the all-spin-up state.

The magnetization in the $z$ direction is defined as
$M_{z}=\langle\psi|\hat{S}^{z}|\psi\rangle$, where $|\psi\rangle$ is a
wavefunction, and

$\displaystyle\hat{S}^{z}=\sum_{j=1}^{L}\sigma_{j}^{z}.$ (3)

Since $\hat{S}^{z}$ is diagonal in the computational basis, its eigenstates
are given by the computational basis states $|k\rangle$, which satisfy

$\displaystyle\hat{S}^{z}|k\rangle=m_{k}|k\rangle,\quad
m_{k}=\sum_{j=1}^{L}s_{j}^{(k)}.$ (4)

Here, $m_{k}$ is the eigenvalue corresponding to $|k\rangle$. The spin
configuration corresponding to $|k\rangle$ is described by
$\bm{s}^{(k)}\in\\{-1,1\\}^{L}$, whose $j$th element is $s_{j}^{(k)}$. In the
quantum-circuit simulation and real-device computation, the magnetization is
expressed as $M_{z}=\sum_{k}p_{k}m_{k}$, where $p_{k}$ is the experimental
probability associated with the spin configuration $\bm{s}^{(k)}$.

The magnetization, which is initially equal to the system size $L$, relaxes to
zero because of antiferromagnetic interactions in a weak-disorder case. In
contrast, the magnetization in a strong-disorder case is expected to retain
its initial value owing to the memory effect.

The twist overlap $z=\langle\psi|\hat{U}_{\rm twist}|\psi\rangle$ is the
overlap between a state $|\psi\rangle$ and its twisted state $\hat{U}_{\rm
twist}|\psi\rangle$. Here, $\hat{U}_{\rm twist}$ is the twist operator defined
by

$\hat{U}_{\rm
twist}=\exp\left[i\sum_{j=1}^{L}\theta_{j}\sigma^{z}_{j}/2\right],$ (5)

where $\theta_{j}=2\pi j/L$. Since the twist operator is also diagonal in the
computational basis, we have

$\displaystyle\hat{U}_{\rm twist}|k\rangle=\exp(iu_{k})|k\rangle,\quad
u_{k}=\frac{\pi}{L}\sum_{j=1}^{L}js_{j}^{(k)}.$ (6)

In the quantum-circuit simulation and real-device computation, the twist
overlap is given by $z=\sum_{k}p_{k}\exp(iu_{k})$. Consequently, both the
twist overlap and the magnetization are derived from an identical set of
measurement outcomes.

The twist operator generates a spin-wave-like excitation through the rotation
of spins around the $z$-axis at angles $\theta_{j}$ [45, 46]. When
$|\psi\rangle$ is a thermal eigenstate, the twist overlap vanishes because the
twisted state with a spin-wave-like excitation is orthogonal to the original
state. However, the long-wavelength perturbation given by the twist operator
has little effect on localized eigenstates, leading to a finite twist overlap.
The twist overlap can also be evaluated using a wavefunction that is a linear
combination of eigenstates.

###  II.2 Time Evolution

Time evolution is based on the Schrödinger equation

$i\frac{d}{dt}|\psi(t)\rangle=H|\psi(t)\rangle,$ (7)

whose solution is represented as

$|\psi(t)\rangle=e^{-iHt}|\psi_{0}\rangle,$ (8)

where $|\psi_{0}\rangle$ denotes the initial state. The wavefunction
$|\psi(t)\rangle$ at any time can be obtained through the exact
diagonalization of the Hamiltonian. Writing the initial state as a liner
combination $|\psi_{0}\rangle=\sum_{k}c_{k}|\phi_{k}\rangle$ of the
eigenstates $|\phi_{k}\rangle$ yields

$|\psi(t)\rangle=\sum_{k=1}^{2^{L}}c_{k}e^{-iE_{k}t}|\phi_{k}\rangle,$ (9)

where $E_{k}$ is the eigenenergy corresponding to $|\phi_{k}\rangle$. We here
refer to the time dependence of wavefunctions obtained from the exact
diagonalization as the exact time evolution.

For quantum simulations using quantum circuits, the time evolution operator
should be approximated using product formulas [41, 42]. Using the first-order
Suzuki–Trotter decomposition, we have

$\displaystyle e^{-iH\tau}=e^{-i(H_{\rm I}+H_{\rm TI})\tau}=\left(e^{-iH_{\rm
I}\delta}e^{-iH_{\rm TI}\delta}\right)^{m},$ (10)

where $\delta=\tau/m$ and $m$ is called the Trotter number or Trotter steps.
Similarly, the second-order Suzuki–Trotter decomposition leads to

$\displaystyle e^{-iH\tau}=\left(e^{-iH_{\rm TI}\delta/2}e^{-iH_{\rm
I}\delta}e^{-iH_{\rm TI}\delta/2}\right)^{m}.$ (11)

The Trotter error for the $p$th-order decomposition scales as
$O(\delta^{p+1})$. Single-qubit rotation gates applied to the $j$th qubit
denoted by $R^{X}_{j}(\theta)$ and $R^{Z}_{j}(\theta)$ implement the unitary
operators $e^{-i(\theta/2)\sigma_{j}^{x}}$ and
$e^{-i(\theta/2)\sigma_{j}^{z}}$, respectively. Similarly, two-qubit $ZZ$
rotation gates implement operators for interaction terms.

###  II.3 Real-Device Computation

We utilized ibmq_manila, one of the IBM Quantum Falcon Processors, to perform
computation on a real quantum computer. The machine consists of five qubits
arranged in a linear structure.

The straightforward execution of quantum circuits on noisy quantum hardware
causes errors. To reduce noise-induced errors without relying on fault
tolerance strategies, we utilized a measurement mitigation technique called
matrix-free measurement mitigation (M3), which is implemented through the
mthree Qiskit package [47]. M3 computes corrected measurement probabilities
within a reduced subspace, using either the lower-upper (LU) decomposition or
a preconditioned iterative method.

We evaluated observables by measurements with 1024 shots per circuit in both
quantum-circuit simulation and real-device computation.

##  III Disorder Dependence

We measured the magnetization and twist overlap at the end of time evolution,
which we fixed at $T_{\rm fin}=1.5$ in this study. This time is sufficiently
long, especially for cases with middle to strong disorders, as supported by
the time evolution results presented in Sect. IV. We set the system size to
$L=5$ for quantum-circuit simulation, real-device computation, and exact time
evolution to facilitate comparison.

Figure 1:  (Color online) Disorder dependence of (a,b) the magnetization in
the $z$-direction $M_{z}$ and (c,d) the absolute square of the twist overlap
$|z|^{2}$. Panels (a) and (c) show the average, and panels (b) and (d) show
the standard deviation. The data obtained from the exact time evolution are
labeled “exact”. The simulation results obtained from the $p$th-order
Suzuki–Trotter decomposition with the Trotter number $m$ are labeled $(p,m)$.
The data evaluated with eigenstates of the Hamiltonian are labeled “ES”.

The results of quantum-circuit simulations, computed using a classical
computer, are shown in Fig. 1. The simulation results labeled $(p,m)$ are
obtained using the $p$th-order Suzuki–Trotter decomposition with the Trotter
number $m$. The data are averaged over $10^{4}$ disorder realization samples.
The results labeled “exact” are obtained from the exact time evolution. The
disorder dependence of $M_{z}$ in the simulations is almost identical to that
in the exact computation.

The twist overlap simulation results are compared with those obtained from the
exact time evolution and the ones evaluated with eigenstates (labeled “ES”) of
the Hamiltonian. The twist overlap evaluated with eigenstates in each
realization is calculated and averaged over 16 eigenstates around the center
of the energy spectrum. The average of $|z|^{2}$ in the simulations shows a
trend similar to the results obtained from the exact time evolution and the
ones evaluated with eigenstates. However, the standard deviation
$\delta|z|^{2}$ for the results obtained from eigenstates is smaller than the
others. The simulations based on the first-order Suzuki–Trotter decomposition
yields a smaller $|z|^{2}$ than the exact time evolution, whereas the second-
order decomposition leads to a larger $|z|^{2}$.

Figure 2:  (Color online) Disorder dependence of (a,b) the magnetization in
the $z$-direction $M_{z}$ and (c,d) the absolute square of the twist overlap
$|z|^{2}$. Panels (a) and (c) show the average, and panels (b) and (d) show
the standard deviation. The real-device computation results for the $p$th-
order Suzuki–Trotter decomposition with the Trotter number $m$ are labeled
$(p,m)$. The data labeled “exact” and “ES” are the same as those in Fig. 1.

The real-device computation results shown in Fig. 2 indicate that noise-
induced errors are significant. The data labeled $(p,m)$ are computed using
the $p$th-order Suzuki–Trotter decomposition with the Trotter number $m$ and
averaged over $10^{3}$ disorder realization samples. The magnetization
obtained from real-device computation increases with disorder strength.
However, the average $M_{z}$ is smaller than the result obtained from the
exact time evolution, particularly in the large-$w$ region. The average and
standard deviation of $|z|^{2}$ obtained from real-device computation also
increase with disorder strength, although they are much smaller than those
obtained from the exact time evolution. The errors in $M_{z}$ and $|z|^{2}$
are more significant for $m=15$ than for $m=10$, whereas the difference
between $p=1$ and $p=2$ is negligible.

##  IV Time Evolution

Figure 3:  (Color online) Time evolution of (a,b) the magnetization in the
$z$-direction $M_{z}$ and (c,d) the absolute square of the twist overlap
$|z|^{2}$ for a system size of $L=5$. Each panel displays ten different
samples. Panels (a) and (c) share the same samples of wavefunctions, and the
same applies to panels (b) and (d).

Figure 3 shows the time evolution of the magnetization in the $z$-direction
($M_{z}$) and the absolute square of the twist overlap ($|z|^{2}$) for ten
different disorder realization samples, calculated using the exact time
evolution for a system with $L=5$. Figures 3(a) and 3(c) share the same
samples of wavefunctions, implying that each sample of $M_{z}$ is obtained
from the same realization as the corresponding sample of $|z|^{2}$. The same
applies to Figs. 3(b) and 3(d). In a weak disorder case ($w=1$), the initial
decay of $M_{z}$ ceases around $t\simeq 1$, whereas that of $|z|^{2}$ does
before $t=1$. These results suggest that the final time should be $T_{\rm
fin}>1$. To reduce noise-induced errors, a short $T_{\rm fin}$ is preferable.
Thus, we set $T_{\rm fin}=1.5$ for computation in Sect. III.

Although both the magnetization and the twist overlap tend to relax to zero in
the thermal regime ($w=1$), they exhibit different behaviors in the localized
regime ($w=10$). The magnetization in the localized regime fluctuates around
the initial value owing to the memory effect, whereas $|z|^{2}$ does not
maintain the initial value. The time average of $|z|^{2}$ is less than unity,
although it is still much larger than that in the thermal regime, and the
variation in time average is notable. Moreover, some samples of $|z|^{2}$ show
large-amplitude oscillations in the localized regime. These properties are
reflected in the large variance of the twist overlap shown in Fig. 1(c).

To understand the difference, we theoretically analyze the time evolution of
the system, focusing on the localized regime. Using Eq. (9), we describe the
time dependence of the magnetization and twist overlap as follows:

$\displaystyle M_{z}(t)$
$\displaystyle=\sum_{k}|c_{k}|^{2}\langle\phi_{k}|\hat{S}^{z}|\phi_{k}\rangle$
$\displaystyle\quad+\sum_{k\neq
l}c_{k}c_{l}^{*}\langle\phi_{l}|\hat{S}^{z}|\phi_{k}\rangle
e^{-i(E_{k}-E_{l})t},$ (12) $\displaystyle z(t)$
$\displaystyle=\sum_{k}|c_{k}|^{2}\langle\phi_{k}|\hat{U}_{\rm
twist}|\phi_{k}\rangle$ $\displaystyle\quad+\sum_{k\neq
l}c_{k}c_{l}^{*}\langle\phi_{l}|\hat{U}_{\rm twist}|\phi_{k}\rangle
e^{-i(E_{k}-E_{l})t}.$ (13)

The first terms on the right-hand side of Eqs. (12) and (13) represent long-
time averages, whereas the second terms correspond to time-dependent parts.

In the MBL phase, the eigenstates of the Hamiltonian are approximately
expressed as product states of local integrals of motion (LIOMs) [48, 6, 38].
Assuming that these states are also sufficiently close to the computational
basis in the localized regime considered here, we describe the eigenstate
$|\phi_{k}\rangle$ corresponding to the basis state $|k\rangle$ as a linear
combination of $|k\rangle$ and a state perpendicular to $|k\rangle$, i.e.,

$\displaystyle|\phi_{k}\rangle=b_{k}|k\rangle+|\alpha_{k}\rangle,$ (14)

where $\langle k|\alpha_{k}\rangle=0$. We also assume that
$\left|\langle\alpha_{k}|\alpha_{l}\rangle\right|\ll 1$ for $l\neq k$.

Using Eqs. (12) and (14) and assuming
$\left|\langle\alpha_{k}|\alpha_{l}\rangle\right|\ll 1$, we approximate the
time dependence of magnetization as

$\displaystyle M_{z}(t)$
$\displaystyle\simeq\sum_{k}|c_{k}|^{2}|b_{k}|^{2}m_{k}+2\sum_{k}m_{k}g_{k}(t),$
(15)

where

$\displaystyle g_{k}(t)$ $\displaystyle=\sum_{l(\neq
k)}|A_{kl}|\cos[(E_{k}-E_{l})t-\beta_{kl}],$ (16)

with $A_{kl}=c_{k}c_{l}^{*}b_{k}$. Here, $\beta_{kl}$ is defined through
$\tan\beta_{kl}=(\mathrm{Im}\;A_{kl})/(\mathrm{Re}\;A_{kl})$. Similarly, using
Eqs. (13) and (14) and $\left|\langle\alpha_{k}|\alpha_{l}\rangle\right|\ll
1$, we approximate $|z(t)|^{2}$ as

$\displaystyle|z(t)|^{2}$
$\displaystyle\simeq|z_{0}|^{2}+4|z_{0}|\sum_{k}\cos(u_{k}-\beta_{0})g_{k}(t)$
$\displaystyle\quad+\left|\sum_{k}e^{iu_{k}}g_{k}(t)\right|^{2},$ (17)

where $z_{0}=\sum_{k}|c_{k}|^{2}\langle\phi_{k}|\hat{U}_{\rm
twist}|\phi_{k}\rangle$ and $\beta_{0}$ is defined through
$\tan\beta_{0}=(\mathrm{Im}\;z_{0})/(\mathrm{Re}\;z_{0})$.

The long-time averages of $M_{z}$ and $|z|^{2}$ in the localized regime are
approximated as follows:

$\displaystyle\overline{M_{z}(t)}$
$\displaystyle\simeq\sum_{k}|c_{k}|^{2}|b_{k}|^{2}m_{k},$ (18)
$\displaystyle\overline{|z(t)|^{2}}$
$\displaystyle\simeq\sum_{k}|c_{k}|^{4}|b_{k}|^{4}.$ (19)

Here, we have neglected the terms related to $|A_{kl}|^{2}$, which are assumed
to be negligible for $k\neq l$. Since $\overline{M_{z}(t)}\simeq L$ in the
localized regime, as shown in Fig. 3(b), we can assume that
$|c_{k}|^{2}\approx 1$ and $|b_{k}|^{2}\approx 1$ at $k$, where $m_{k}=L$.
Additionally, since $|c_{k}|^{4}<|c_{k}|^{2}$ and $|b_{k}|^{4}<|b_{k}|^{2}$,
we can also understand why $\overline{|z(t)|^{2}}$ is slightly smaller than
unity.

Figure 4:  (Color online) (a,b) Plots of the overlap $|b_{k}|^{2}$ between the
basis state $|k\rangle$ and the corresponding eigenstate $|\phi_{k}\rangle$ of
the Hamiltonian. (c,d) Plots of the overlap $|c_{k}|^{2}$ between the
eigenstate $|\phi_{k}\rangle$ and the initial state $|\psi_{0}\rangle$. Each
panel displays ten samples that are obtained from the same samples of
wavefunctions as in Fig. 3.

Our assumptions in the localized regime are supported numerically. Figure 4
shows the overlap $|b_{k}|^{2}=|\langle k|\phi_{k}\rangle|^{2}$ between the
basis state $|k\rangle$ and the corresponding eigenstate $|\phi_{k}\rangle$ of
the Hamiltonian, as well as the overlap
$|c_{k}|^{2}=|\langle\phi_{k}|\psi_{0}\rangle|^{2}$ between the eigenstate
$|\phi_{k}\rangle$ and the initial state $|\psi_{0}\rangle$, for $w=1$ and
$w=10$. Here, we index the eigenstate $|\phi_{k}\rangle$ of the Hamiltonian as
$k=\mathrm{argmax}_{i}|\langle k|\phi_{i}\rangle|^{2}$. Note that
$|b_{0}|^{2}=|c_{0}|^{2}$ because the initial state is equivalent to the basis
state of $k=0$. In several samples for $w=10$, we observe that
$|b_{k}|^{2}\simeq 1$, which is consistent with the assumption that the
eigenstates are close to the computational basis in the localized regime.

The distributions of $|b_{k}|^{2}$ and $|c_{k}|^{2}$ also explain the behavior
of time averages in the localized regime. Whereas $|b_{k}|^{2}\simeq 1$ in
typical samples for $w=10$, $|b_{k}|^{2}\sim 0.5$ in some samples, although
$|b_{k}|^{2}$ is still mostly larger than that for $w=1$. The variation in
$|b_{k}|^{2}$ is related to that in $|z|^{2}$. Additionally, $|c_{k}|^{2}$ has
the largest peak at $k=0$, which corresponds to the initial state, and for
most samples with $w=10$, $|c_{0}|^{2}\approx 1$ and $|c_{k}|^{2}\ll 1$ for
$k\neq 0$. These observations explain why $\overline{M_{z}(t)}\simeq m_{0}=L$
and $\overline{|z(t)|^{2}}\simeq 1$ for typical samples in the localized
regime.

Here, we analyze the time-dependent parts of the magnetization and twist
overlap given by $f_{\rm mg}(t)$ and $f_{\rm tw}(t)$, respectively, which are
expressed as

$\displaystyle f_{\rm mg}(t)$ $\displaystyle=2\sum_{k}m_{k}g_{k}(t),$ (20)
$\displaystyle f_{\rm tw}(t)$
$\displaystyle=4|z_{0}|\sum_{k}\cos(u_{k}-\beta_{0})g_{k}(t).$ (21)

We have neglected the terms including $|A_{kl}|^{2}$, which are assumed to be
negligible for $k\neq l$. Since both $f_{\rm mg}(t)$ and $f_{\rm tw}(t)$
depend on the same time-dependent function $g_{k}(t)$, they have the same
period and phase of oscillation. However, as seen in Figs. 3(b) and 3(d),
their oscillation amplitudes differ considerably in the localized regime. The
difference originates from the factors $m_{k}$ and $\cos(u_{k}-\beta_{0})$. As
$m_{k}=0$ for many basis states, only a few elements of $g_{k}$ contribute to
$f_{\rm mg}$. The lack of contributions from the state with $m_{k}=0$ explains
the reason for the relatively small amplitudes of magnetization oscillation
seen in Fig. 3(b). In contrast, $\cos(u_{k}-\beta_{0})$ can take various
values within the range $[-1,1]$. Most elements of $g_{k}$ contribute to
$f_{\rm tw}$, interfering constructively or destructively. The variety in this
factor explains the various amplitudes of $|z(t)|^{2}$ in Fig. 3(d).

The above theoretical and numerical explanations are based on the exact
diagonalization of the Hamiltonian. They capture the crucial aspects of the
difference between the magnetization and the twist overlap.

##  V Discussion

###  V.1 Trotter Error and Noise-Induced Errors

The numbers of single- and two-qubit gates used in real-device computation
depend on the basis gates supported by the specific quantum device employed.
In the real device used in this work, an RZZ gate is implemented as an RZ gate
wrapped by a pair of CNOT gates, whereas an RX gate is substituted with a
combination of four single-qubit gates. In total, for $n$ qubits and $m$
Trotter steps, the numbers of single-qubit gates in the first- and second-
order decompositions are $(6n-1)m$ and $(6n-1)m+5n$, respectively. The number
of two-qubit gates is $2(n-1)m$ for both the first- and second-order
decompositions.

In the real device used in this work, single- and two-qubit gates are
calibrated frequently. The errors of CNOT and $\sqrt{X}$ gates were on the
order of $10^{-3}$ and $10^{-4}$, respectively, and the readout error was on
the order of $10^{-2}$.

In this study, we have chosen Trotter numbers so that the Trotter errors in
the quantum-circuit simulation are acceptable for evaluating the magnetization
and twist overlap. Increasing the Trotter number improves the approximation
accuracy, but also increases noise-induced errors in real-device computation.
Two-qubit gates contribute larger errors than single-qubit gates, and their
number increases almost linearly with the Trotter number, regardless of the
order of the Suzuki–Trotter decomposition used. Thus, noise-induced errors are
more sensitive to the Trotter number than the order of the Suzuki–Trotter
decomposition.

###  V.2 Difference in Effect of Errors

As shown in Fig. 2, the twist overlap suffers from noise-induced errors more
severely than the magnetization. In the localized regime, the real-device
computation results for $M_{z}$ are approximately $3/4$ of the exact time
evolution result, whereas those of $|z|^{2}$ are less than $1/3$ of the exact
time evolution result. The decay ratio of the average and standard deviation
is almost the same, suggesting that noise-induced errors are more prominent in
the time average part than in the oscillating parts. Considering that
$\overline{M_{z}(t)}$ is estimated as Eq. (18), we can suppose that noise-
induce errors reduce the dominant elements of $|c_{k}|^{2}|b_{k}|^{2}$ to
$3/4$ of the expected value and that these dominant elements correspond to
states with large $m_{k}$ values. On the other hand, the $\overline{z(t)}$ of
real-device computation results, which is estimated as Eq. (19), is much less
than $(3/4)^{2}$ of the expected value. The unexpectedly small values of
$\overline{z(t)}$ imply that errors in the states with small $m_{k}$ values
affect the twist overlap.

The error in $|z|^{2}$ is approximately three times larger than that in
$M_{z}$. As both the twist overlap and the magnetization are derived from
measurements of individual qubits, this relationship remains independent of
the noise strength. If we manage to reduce the error in $M_{z}$ to around 3%,
the error in $|z|^{2}$ becomes approximately 10%. The reduction in error makes
the twist overlap a potentially more informative quantity.

###  V.3 Information on Eigenstates

Figure 1(c) illustrates that the disorder dependence of $|z|^{2}$ evaluated
with the wavefunction after time evolution is similar to that evaluated with
the eigenstates of the Hamiltonian. The twist overlap $z_{e}$ evaluated with
the eigenstates is described as

$\displaystyle|z_{e}|^{2}=\frac{1}{|B|}\sum_{k\in
B}\left|\langle\phi_{k}|\hat{U}_{\rm
twist}|\phi_{k}\rangle\right|^{2}\simeq\frac{1}{|B|}\sum_{k\in B}|b_{k}|^{4},$
(22)

where $B$ denotes the index set of the eigenstates used for the evaluation and
$|B|$ is its size. Equation (22) indicates that $|z_{e}|^{2}$ approximately
represents the simple average of $|b_{k}|^{4}$ in set $B$ in the localized
regime. In contrast, $\overline{|z(t)|^{2}}$, which is estimated as Eq. (19)
in the localized regime, corresponds to an average weighted with the overlap
between eigenstates and the initial state. In the localized regime, where
$|b_{k}|^{2}\approx 1$ and $|c_{0}|^{2}\approx 1$ for a typical sample,
$|z_{e}|^{2}\simeq\overline{|z(t)|^{2}}$ on average.

In contrast, as shown in Fig. 1(d), the variance of $|z|^{2}$ is more
significant for the twist overlap evaluated after time evolution than for that
evaluated with eigenstates. This increase in variance is due to the
oscillating part of $|z(t)|^{2}$. Since the oscillation amplitude is
significant in some samples, the twist overlap can have various values at the
end of time evolution, leading to a significant variance.

In the quantum-circuit simulation and real-device computation, the evaluation
of the twist overlap is based on the measurement of each qubit. It is not a
straightforward summation of measured values ($s_{j}$); rather, it is
expressed as an infinite series in terms of $s_{j}$. An alternative method for
evaluating the twist overlap is to employ the Hadamard test between
$|\psi(t)\rangle$ and $\hat{U}_{\rm twist}|\psi(t)\rangle$, which includes the
measurement of only a single ancilla qubit. This approach is applicable in
quantum-circuit simulation as gate errors are excluded. However, the Hadamard
test involves control-$U$ operations, potentially resulting in notable two-
qubit gate errors in real-device computations.

##  VI Conclusions

In this study, we investigated the method to detect MBL based on quantum
dynamics, focusing on magnetization and twist overlap. We examined the
disorder dependence of these observables and analyzed their time evolution in
the localized regime both theoretically and numerically. By analyzing the time
evolution, we gained insights into how noise-induced errors affect the
observables.

The magnetization is simple and robust against noise. The disorder dependence
of the magnetization clearly shows the crossover between thermal and localized
regimes. The magnetization is evaluated from measured values of spin
configurations ($s_{j}$), which directly correspond to the eigenvalues
($\sum_{j}s_{j}$) of the magnetization operator. Furthermore, the
magnetization has a large expectation value for the initial state, leading to
a significant contribution of the initial state and small contributions of
other states. This type of observable is robust against noise because states
with small expectation values do not significantly contribute to the
observable.

On the other hand, the twist overlap is vulnerable to noise-induced errors,
although the disorder dependence exhibits a weak signature of the crossover in
real-device computation results. The evaluation of the twist overlap relies on
the measured values ($s_{j}$) through the phases ($u=\sum_{j}(\pi j/L)s_{j}$)
of the eigenvalues ($\exp(iu)$) of the twist operator. However, the
information of phases is lost in the absolute value of the twist overlap, and
that of eigenstates of the operator emerges. Furthermore, every eigenstate of
the Hamiltonian can contribute to the twist overlap, implying that noise-
induced errors are not systematically reduced, unlike the magnetization case.
Despite being vulnerable to noise, this type of observable can provide richer
information about the eigenstates of the Hamiltonian than other observables
such as magnetization under a noiseless condition.

In conclusion, quantum simulations using a real noisy device can detect the
crossover between thermal and localized regimes, provided that appropriate
observables are chosen. Some observables, such as the magnetization, are
robust against noise-induced errors, whereas others, such as the twist
overlap, are vulnerable. The results of this study provide helpful insights
into the selection of observables for probing MBL using a real quantum
computer.

###### Acknowledgements.

We acknowledge the use of IBM Quantum services for this work. The views
expressed are those of the author and do not reflect the official policy or
position of IBM or the IBM Quantum team. This work was partially supported by
the JSPS KAKENHI Grant Number JP23H04499.

## References

  * Parameswaran and Vasseur [2018] S. A. Parameswaran and R. Vasseur, Many-body localization, symmetry and topology, Rep. Prog. Phys. 81, 082501 (2018).
  * Abanin  _et al._ [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement, Rev. Mod. Phys. 91, 021001 (2019).
  * Gopalakrishnan and Parameswaran [2020] S. Gopalakrishnan and S. A. Parameswaran, Dynamics and transport at the threshold of many-body localization, Phys. Rep. 862, 1 (2020).
  * Tikhonov and Mirlin [2021] K. S. Tikhonov and A. D. Mirlin, From Anderson localization on random regular graphs to many-body localization, Ann. Phys. Special Issue on Localisation 2020, 435, 168525 (2021).
  * Pal and Huse [2010] A. Pal and D. A. Huse, Many-body localization phase transition, Phys. Rev. B 82, 174411 (2010).
  * Serbyn  _et al._ [2013a] M. Serbyn, Z. Papić, and D. A. Abanin, Local Conservation Laws and the Structure of the Many-Body Localized States, Phys. Rev. Lett. 111, 127201 (2013a).
  * Serbyn and Moore [2016] M. Serbyn and J. E. Moore, Spectral statistics across the many-body localization transition, Phys. Rev. B 93, 041424(R) (2016).
  * Khemani  _et al._ [2016] V. Khemani, F. Pollmann, and S. L. Sondhi, Obtaining Highly Excited Eigenstates of Many-Body Localized Hamiltonians by the Density Matrix Renormalization Group Approach, Phys. Rev. Lett. 116, 247204 (2016).
  * Khemani  _et al._ [2017a] V. Khemani, S. P. Lim, D. N. Sheng, and D. A. Huse, Critical Properties of the Many-Body Localization Transition, Phys. Rev. X 7, 021013 (2017a).
  * Khemani  _et al._ [2017b] V. Khemani, D. N. Sheng, and D. A. Huse, Two Universality Classes for the Many-Body Localization Transition, Phys. Rev. Lett. 119, 075702 (2017b).
  * Bera  _et al._ [2015] S. Bera, H. Schomerus, F. Heidrich-Meisner, and J. H. Bardarson, Many-Body Localization Characterized from a One-Particle Perspective, Phys. Rev. Lett. 115, 046603 (2015).
  * Bera and Lakshminarayan [2016] S. Bera and A. Lakshminarayan, Local entanglement structure across a many-body localization transition, Phys. Rev. B 93, 134204 (2016).
  * Kjäll  _et al._ [2014] J. A. Kjäll, J. H. Bardarson, and F. Pollmann, Many-Body Localization in a Disordered Quantum Ising Chain, Phys. Rev. Lett. 113, 107204 (2014).
  * Luitz  _et al._ [2015] D. J. Luitz, N. Laflorencie, and F. Alet, Many-body localization edge in the random-field Heisenberg chain, Phys. Rev. B 91, 081103(R) (2015).
  * Hopjan and Heidrich-Meisner [2020] M. Hopjan and F. Heidrich-Meisner, Many-body localization from a one-particle perspective in the disordered one-dimensional bose-hubbard model, Phys. Rev. A 101, 063617 (2020).
  * Gray  _et al._ [2018] J. Gray, S. Bose, and A. Bayat, Many-body localization transition: Schmidt gap, entanglement length, and scaling, Phys. Rev. B 97, 201105(R) (2018).
  * Kudo and Deguchi [2018] K. Kudo and T. Deguchi, Finite-size scaling with respect to interaction and disorder strength at the many-body localization transition, Phys. Rev. B 97, 220201(R) (2018).
  * Bardarson  _et al._ [2012] J. H. Bardarson, F. Pollmann, and J. E. Moore, Unbounded Growth of Entanglement in Models of Many-Body Localization, Phys. Rev. Lett. 109, 017202 (2012).
  * Serbyn  _et al._ [2013b] M. Serbyn, Z. Papić, and D. A. Abanin, Universal Slow Growth of Entanglement in Interacting Strongly Disordered Systems, Phys. Rev. Lett. 110, 260601 (2013b).
  * Enss  _et al._ [2017] T. Enss, F. Andraschko, and J. Sirker, Many-body localization in infinite chains, Phys. Rev. B 95, 045121 (2017).
  * Harris [1974] A. B. Harris, Effect of random defects on the critical behaviour of ising models, J. Phys. C: Solid State Phys. 7, 1671 (1974).
  * Chayes  _et al._ [1986] J. T. Chayes, L. Chayes, D. S. Fisher, and T. Spencer, Finite-size scaling and correlation lengths for disordered systems, Phys. Rev. Lett. 57, 2999 (1986).
  * Imbrie [2016] J. Z. Imbrie, On Many-Body Localization for Quantum Spin Chains, J. Stat. Phys. 163, 998 (2016).
  * Sierant and Zakrzewski [2022] P. Sierant and J. Zakrzewski, Challenges to observation of many-body localization, Phys. Rev. B 105, 224203 (2022).
  * Morningstar  _et al._ [2022] A. Morningstar, L. Colmenarez, V. Khemani, D. J. Luitz, and D. A. Huse, Avalanches and many-body resonances in many-body localized systems, Phys. Rev. B 105, 174205 (2022).
  * Sels [2022] D. Sels, Bath-induced delocalization in interacting disordered spin chains, Phys. Rev. B 106, L020202 (2022).
  * Schreiber  _et al._ [2015] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Observation of many-body localization of interacting fermions in a quasirandom optical lattice, Science 349, 842 (2015).
  * Smith  _et al._ [2016] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, Many-body localization in a quantum simulator with programmable random disorder, Nat. Phys. 12, 907 (2016).
  * Bordia  _et al._ [2016] P. Bordia, H. P. Lüschen, S. S. Hodgman, M. Schreiber, I. Bloch, and U. Schneider, Coupling identical one-dimensional many-body localized systems, Phys. Rev. Lett. 116, 140401 (2016).
  * Bordia  _et al._ [2017] P. Bordia, H. Lüschen, S. Scherg, S. Gopalakrishnan, M. Knap, U. Schneider, and I. Bloch, Probing Slow Relaxation and Many-Body Localization in Two-Dimensional Quasiperiodic Systems, Phys. Rev. X 7, 041047 (2017).
  * Lüschen  _et al._ [2017a] H. P. Lüschen, P. Bordia, S. S. Hodgman, M. Schreiber, S. Sarkar, A. J. Daley, M. H. Fischer, E. Altman, I. Bloch, and U. Schneider, Signatures of Many-Body Localization in a Controlled Open Quantum System, Phys. Rev. X 7, 011034 (2017a).
  * Lüschen  _et al._ [2017b] H. P. Lüschen, P. Bordia, S. Scherg, F. Alet, E. Altman, U. Schneider, and I. Bloch, Observation of Slow Dynamics near the Many-Body Localization Transition in One-Dimensional Quasiperiodic Systems, Phys. Rev. Lett. 119, 260401 (2017b).
  * Xu  _et al._ [2018] K. Xu, J.-J. Chen, Y. Zeng, Y.-R. Zhang, C. Song, W. Liu, Q. Guo, P. Zhang, D. Xu, H. Deng, K. Huang, H. Wang, X. Zhu, D. Zheng, and H. Fan, Emulating Many-Body Localization with a Superconducting Quantum Processor, Phys. Rev. Lett. 120, 050507 (2018).
  * Kohlert  _et al._ [2019] T. Kohlert, S. Scherg, X. Li, H. P. Lüschen, S. Das Sarma, I. Bloch, and M. Aidelsburger, Observation of Many-Body Localization in a One-Dimensional System with a Single-Particle Mobility Edge, Phys. Rev. Lett. 122, 170403 (2019).
  * Rubio-Abadal  _et al._ [2019] A. Rubio-Abadal, J.-y. Choi, J. Zeiher, S. Hollerith, J. Rui, I. Bloch, and C. Gross, Many-Body Delocalization in the Presence of a Quantum Bath, Phys. Rev. X 9, 041014 (2019).
  * Gong  _et al._ [2021] M. Gong, G. D. de Moraes Neto, C. Zha, Y. Wu, H. Rong, Y. Ye, S. Li, Q. Zhu, S. Wang, Y. Zhao, F. Liang, J. Lin, Y. Xu, C.-Z. Peng, H. Deng, A. Bayat, X. Zhu, and J.-W. Pan, Experimental characterization of the quantum many-body localization transition, Phys. Rev. Research 3, 033043 (2021).
  * Filho  _et al._ [2022] J. L. C. d. C. Filho, Z. G. Izquierdo, A. Saguia, T. Albash, I. Hen, and M. S. Sarandy, Localization transition induced by programmable disorder, Phys. Rev. B 105, 134201 (2022).
  * Chiaro  _et al._ [2022] B. Chiaro, C. Neill, A. Bohrdt, M. Filippone, F. Arute, K. Arya, R. Babbush, D. Bacon, J. Bardin, R. Barends, S. Boixo, D. Buell, B. Burkett, Y. Chen, Z. Chen, R. Collins, A. Dunsworth, E. Farhi, A. Fowler, B. Foxen, C. Gidney, M. Giustina, M. Harrigan, T. Huang, S. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, J. Kelly, P. Klimov, A. Korotkov, F. Kostritsa, D. Landhuis, E. Lucero, J. McClean, X. Mi, A. Megrant, M. Mohseni, J. Mutus, M. McEwen, O. Naaman, M. Neeley, M. Niu, A. Petukhov, C. Quintana, N. Rubin, D. Sank, K. Satzinger, T. White, Z. Yao, P. Yeh, A. Zalcman, V. Smelyanskiy, H. Neven, S. Gopalakrishnan, D. Abanin, M. Knap, J. Martinis, and P. Roushan, Direct measurement of nonlocal interactions in the many-body localized phase, Phys. Rev. Res. 4, 013148 (2022).
  * Liu  _et al._ [2023] S. Liu, S.-X. Zhang, C.-Y. Hsieh, S. Zhang, and H. Yao, Probing many-body localization by excited-state variational quantum eigensolver, Phys. Rev. B 107, 024204 (2023).
  * Kudo [2022] K. Kudo, Localization Detection Based on Quantum Dynamics, Entropy 24, 1085 (2022).
  * Hatano and Suzuki [2005] N. Hatano and M. Suzuki, Finding Exponential Product Formulas of Higher Orders, in  _Quantum Annealing and Other Optimization Methods_, Lecture Notes in Physics, edited by A. Das and B. K. Chakrabarti (Springer, Berlin, Heidelberg, 2005) pp. 37–68.
  * Childs  _et al._ [2021] A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, Theory of trotter error with commutator scaling, Phys. Rev. X 11, 011020 (2021).
  * Jaderberg  _et al._ [2020] B. Jaderberg, A. Agarwal, K. Leonhardt, M. Kiffner, and D. Jaksch, Minimum hardware requirements for hybrid quantumâclassical DMFT, Quantum Science and Technology 5, 034015 (2020).
  * Jaderberg  _et al._ [2022] B. Jaderberg, A. Eisfeld, D. Jaksch, and S. Mostame, Recompilation-enhanced simulation of electronâphonon dynamics on IBM quantum computers, New Journal of Physics 24, 093017 (2022).
  * Nakamura and Todo [2002] M. Nakamura and S. Todo, Order Parameter to Characterize Valence-Bond-Solid States in Quantum Spin Chains, Phys. Rev. Lett. 89, 077204 (2002).
  * Kutsuzawa and Todo [2022] T. Kutsuzawa and S. Todo, Nested Iterative Shift-invert Diagonalization for Many-body Localization in the Random-field Heisenberg Chain, arXiv:2203.09732 (2022).
  * [47] mthree: Matrix-free measurement mitigation (m3)., accessed: 2023-09-13.
  * Huse  _et al._ [2014] D. A. Huse, R. Nandkishore, and V. Oganesyan, Phenomenology of fully many-body-localized systems, Phys. Rev. B 90, 174202 (2014).