# State Expansion of a Levitated Nanoparticle in a Dark Harmonic Potential

Eric Bonvin  Louisiane Devaud  Massimiliano Rossi  Andrei Militaru  Photonics
Laboratory, ETH Zürich, 8093 Zürich, Switzerland  Quantum Center, ETH Zürich,
Zürich, Switzerland  Lorenzo Dania  Institut für Experimentalphysik,
Universität Innsbruck, Technikerstrasse 25, 6020 Innsbruck, Austria  Dmitry S.
Bykov  Institut für Experimentalphysik, Universität Innsbruck,
Technikerstrasse 25, 6020 Innsbruck, Austria  Oriol Romero-Isart  Institute
for Theoretical Physics, University of Innsbruck, Innsbruck, Austria
Institute for Quantum Optics and Quantum Information, Austrian Academy of
Sciences, Innsbruck, Austria  Tracy E. Northup  Institut für
Experimentalphysik, Universität Innsbruck, Technikerstrasse 25, 6020
Innsbruck, Austria  Lukas Novotny  Martin Frimmer  Photonics Laboratory, ETH
Zürich, 8093 Zürich, Switzerland  Quantum Center, ETH Zürich, Zürich,
Switzerland

###### Abstract

Levitated nanoparticles in vacuum are prime candidates for generating
macroscopic quantum superposition states of massive objects. Most protocols
for preparing these states necessitate coherent expansion beyond the scale of
the zero-point motion to produce sufficiently delocalized and pure phase-space
distributions. Here, we spatially expand and subsequently recontract the
thermal state of a levitated nanoparticle by modifying the stiffness of the
trap holding the particle. We achieve state-expansion factors of 25 in
standard deviation for a particle initially feedback-cooled to a center-of-
mass thermal state of $155\text{\,}\mathrm{mK}$. Our method relies on a hybrid
scheme combining an optical trap, for cooling and measuring the particle’s
motion, with a Paul trap for expanding its state. Consequently, state
expansion occurs devoid of measurement backaction from photon recoil, making
this approach suitable for coherent wavefunction expansion in future
experiments.

### Introduction.

Quantum physics predicts that a massive object can be prepared in a
macroscopic quantum superposition state—a state in which the center-of-mass
position is delocalized over length scales larger than the object’s size [1].
Matter-wave experiments have demonstrated this remarkable phenomenon with
macromolecules as large as $10^{4}$ atomic mass units [2]. Over the past
decade, levitated nanoparticles in vacuum [3, 4] have emerged as promising
candidates to increase the mass of such macroscopic quantum states by at least
four orders of magnitude [5, 6, 7, 8, 9, 10, 11, 12, 13]. The first important
merit of these levitated systems is that their motional state can be purified
by ground-state cooling [14, 15, 16, 17, 18, 19]. Their second important merit
is that levitation offers the opportunity to manipulate the trapping potential
and thereby control the motional state of a levitated object [20, 21, 22, 23].
This second merit sets levitated systems apart from clamped optomechanical
oscillators [24]. In those systems, impressive optomechanical quantum control
has been achieved [25, 26, 27, 28, 29, 30, 31, 32], but the length scale set
by the zero-point motion of the oscillator cannot be tuned in situ.

Control over the trapping potential is a key requirement in both approaches
pursued to expand the wave-function of a levitated nanoparticle. The first
approach relies on purely optical potentials together with free evolution in
the absence of a potential [12]. The second approach aims to exploit evolution
in a dedicated potential in the absence of decoherence due to measurement
back-action inevitably associated with optical interactions [33, 34, 35]. To
this end, this strategy proposes to use dark, non-optical potentials generated
by radio-frequency (RF) electric fields [13, 10]. On the experimental side,
hybrid approaches of optical traps combined with RF (Paul) traps have been
investigated [36, 37, 38, 39, 40, 41, 42] but remain to be exploited for
expansion protocols.

Generating and handling a quantum mechanically pure state of a nanoparticle in
a hybrid trap is a daunting experimental challenge. Nevertheless, expansion
protocols can already be tested on classical states to accelerate future
progress when moving to the quantum regime.

In this work, we expand the classical thermal state of motion of a levitated
nanoparticle by transferring it from an optical trap to a weakly confining
dark potential implemented with a linear Paul trap. After evolution in the
dark potential, the particle is re-trapped optically and its position is
measured, allowing for repetition of the protocol. We demonstrate a maximum
state expansion by a factor of 25 in standard deviation and subsequent
contraction to the starting size in addition to a contribution from heating
during the protocol. This work is an important step towards coherent
wavefunction expansion of levitated nanoparticles in the quantum regime.

((a))

((b))

((c))

Figure 1:  Phase-space representation of a particle in a harmonic potential
undergoing a frequency-jump protocol. (a) The particle is initially in a
thermal state, described by a uniform Gaussian distribution of size $\Delta
u_{0}$. (b) The state evolves in the weaker potential for increasing $t_{d}$.
(c) The maximum size $\Delta u_{\text{max}}$ is reached for $t_{d}=T_{p}/4$.

### Frequency-jump protocol.

We start by illustrating the dynamics of a harmonic oscillator with position
$u$ and velocity $\dot{u}$ under a so-called frequency-jump protocol [43, 44,
45, 46, 47]. In this protocol, we initialize the oscillator in a thermal state
in a (stiff) potential with frequency $\omega_{o}$ and period $T_{o}$ (in this
work implemented with an optical trap). With appropriately normalized axes,
this initial thermal state is described in phase space by a Gaussian
distribution with widths $\Delta u_{0}$ and $\Delta\dot{u}_{0}$ along the
position and velocity axes, respectively, schematically shown in Fig. 1(a). At
time $t=0$, we switch to a weaker potential with frequency
$\omega_{p}<\omega_{o}$ and period $T_{p}$ (implemented with a Paul trap in
the absence of laser light) for a time $t_{d}$. The evolution of the state in
the weaker potential as a function of $t_{d}$ is illustrated in Fig. 1(b). As
$t_{d}$ increases, the state expands spatially until reaching a maximum size
$\Delta u_{\text{max}}$ for $t_{d}=T_{p}/4$, as shown in Fig. 1(c). As $t_{d}$
increases further, the spatial state size returns to its initial value,
reached at $t_{d}=T_{p}/2$. To characterize the phase-space distribution
acquired at time $t_{d}$, we instantaneously switch the potential back to
frequency $\omega_{o}$. Now, the state rotates in phase space with preserved
shape. In this work, we experimentally realize this frequency-jump protocol in
a hybrid Paul-optical trap.

### Experimental setup.

Our setup overlaps an optical trap with a linear Paul trap, in which we
levitate a $177\text{\,}\mathrm{nm}$ diameter silica nanosphere. This
configuration, shown in Fig. 2, allows us to selectively trap the same
particle in either trap. What follows is a simplified description of the setup
(further details in Ref. [48]).

Figure 2:  Experimental setup. A laser beam is focused by a lens to form an
optical trap. The trapped particle’s center-of-mass motion is recorded on a
quadrant photodetector. The laser can be switched on and off with an acousto-
optic modulator (AOM). The optical trap is superimposed with a Paul trap. It
consists of a microfabricated chip held in the focal plane and endcap
electrodes (held at a potential $V_{\mathrm{cap}}$). The two pairs of on-chip
electrodes are driven with RF voltages ($V_{\mathrm{rf,+}}$ and
$V_{\mathrm{rf,-}}$, respectively) for in-plane confinement.

The optical trap is generated by a laser ($x$-polarized, wavelength
$1565\text{\,}\mathrm{nm}$, power $500\text{\,}\mathrm{mW}$) focused by a lens
($\text{NA}=0.77$, optical axis along $z$). The particle’s center-of-mass
oscillation frequencies are $2\pi\times$(44, 58, 10) $\text{\,}\mathrm{kHz}$
along $(x,y,z)$, respectively. The trap can be switched on and off with an
acousto-optic modulator on a timescale faster than $100$ ns. The forward
scattered light is collected by an identical lens, and sent onto a quadrant
photodiode for position detection in the $(x,y,z)$-basis.

The Paul trap consists of a glass substrate micromachined to a wheel-trap
geometry (FEMTOprint, Switzerland) [49, 50, 51, 52]. The metallized electrodes
are oriented along the orthogonal directions $u$ and $v$, which lie in the
focal plane and are tilted at $45^{\circ}$ relative to $x$ and $y$. The two
pairs of RF electrodes are driven with opposite phases at a frequency of
$33\text{\,}\mathrm{kHz}$ and a nominal peak-to-peak voltage of
$400\text{\,}\mathrm{V}$. Confinement along $z$ is achieved by applying a
typical endcap voltage of $70\text{\,}\mathrm{V}$ to the metallic lens
holders. This results in secular center-of-mass modes at $2\pi\times$(6, 6, 3)
$\text{\,}\mathrm{kHz}$ along $(u,v,z)$, respectively. The radial and axial
mode frequencies can be tuned from $0.16\text{\,}\mathrm{kHz}$ by adjusting
the applied voltages. The trap is mounted in a vacuum chamber and kept at a
pressure of $2.2\text{\times}{10}^{-6}\text{\,}\mathrm{mbar}$. Before starting
our experiments, we compensate for stray fields in the trapping region using
the shim electrodes of our trap [48, 53].

### Results.

We execute a frequency-jump protocol by rapidly transferring the particle from
the high-frequency optical trap to the lower frequency Paul trap and back.
Experimentally, this is achieved by switching the optical trap off and back
on, while keeping the Paul trap enabled at all times. More specifically, we
execute the following protocol: (i) for $t<0$, the particle is held in the
optical trap, and initialized to a feedback-cooled thermal state of motion;
(ii) at time $t=0$, feedback cooling is disabled, and the optical trap is
rapidly switched off; (iii) for $0<t<t_{d}$, the particle evolves in the Paul
trap until (iv) at $t=t_{d}$, the optical trap is switched back on to
recapture the particle optically; and (v) the final position of the particle
in phase space is characterized for $t>t_{d}$. In Fig. 3(a), we show the
measured particle position during one realization of this frequency-jump
protocol with $t_{d}=$72\text{\,}\mathrm{\SIUnitSymbolMicro s}$$ (black data-
points). The measurement of the particle’s position for times $t<0$ is largely
impacted by noise, as expected from a state under feedback cooling. In
contrast, the oscillations in the particle’s position are clearly visible at
the end of the protocol ($t>$150\text{\,}\mathrm{\SIUnitSymbolMicro s}$$). We
note that for times $0<t<t_{d}$, no measurement record exists, since the
optical field is off. Furthermore, between $t=72~{}\mu$s and $\approx
150~{}\mu$s, we observe a slow recovery of the recorded time trace after re-
enabling the optical trap due to AC-coupling of the data acquisition system.
To avoid unwanted artefacts in our data analysis, we ignore the first
$120\text{\,}\mathrm{\SIUnitSymbolMicro s}$ of signal upon recapture [open
data points in Fig. 3(a), see [54]].

We use the measurement record to retrodict the position and velocity after
recapture at $t=t_{d}$ in the optical trap (see Supplemental Material [54]).
The retrodicted trajectory is shown as the dashed blue line in Fig. 3(a). The
phase-space distribution before the frequency jump is inferred by retrodicting
the phase-space position of the particle for protocols executed with
$t_{d}=0$. Thus, all phase-space distributions are generated using the
identical data-processing protocol. Another challenge arises from the coupling
of the motion along $x$ and $y$ by the electric fields in the Paul trap, which
are aligned along the $u$ and $v$ axes. Therefore, in the following, we
combine the measurements of $x$ and $y$ to analyze motion along $u$.

((a))

((b))

((c))

((d))

Figure 3:  (a) Recorded position during a single realization of a frequency-
jump protocol (black data points). The optical trap is disabled from $t=0$ to
$t=t_{d}=$72\text{\,}\mathrm{\SIUnitSymbolMicro s}$$, allowing the particle to
evolve in the Paul trap. The oscillation after recapture is fit (dashed blue
line) to estimate the position and velocity of the particle at time of
recapture $t=t_{d}$ (red circle). (b) Initial state of the particle,
reconstructed from 500 repetitions of the protocol with $t_{d}=0$ (bin size
$$1\text{\,}\mathrm{nm}$\times$1\text{\,}\mathrm{nm}$$). (c) For
$t_{d}=$72\text{\,}\mathrm{\SIUnitSymbolMicro s}$\approx T_{p}/4$, the state
has visibly expanded along the position axis. (d) For
$t_{d}=$140\text{\,}\mathrm{\SIUnitSymbolMicro s}$\approx T_{p}/2$, the state
has contracted back to a size similar to the initial distribution.

A phase-space distribution of the particle before the frequency jump, i.e., at
$t_{d}=0$, generated from 500 repetitions of the experiment is shown in Fig.
3(b). Each repetition is started at the same time in the micromotion cycle.
Here, the particle is cooled to a state size of $1.5\pm
0.1\text{\,}\mathrm{nm}$ in the optical trap. This corresponds to a center-of-
mass temperature of $T_{\text{CoM}}=$155\pm 25\text{\,}\mathrm{mK}$$. The
state size $\Delta u$ is defined as the standard deviation of the position
$u$. We characterize the measurement noise independently and subtract it from
the signal when evaluating $\Delta u$. The error is evaluated using the
bootstrap method [55]. Figure 3(c) shows the distribution after an evolution
time $t_{d}=$72\text{\,}\mathrm{\SIUnitSymbolMicro s}$$, corresponding to a
quarter oscillation period in the Paul trap. In comparison to the initial
state, the distribution for $t_{d}=$72\text{\,}\mathrm{\SIUnitSymbolMicro s}$$
is stretched along the position axis, reaching a state size of $26.4\pm
1.2\text{\,}\mathrm{nm}$. This state expansion in a dark potential is the main
result of this Letter. In Fig. 3(d), we show the distribution for $t=140~{}\mu
s\approx T_{p}/2$, where the state has contracted to $\Delta u=$4.4\pm
0.2\text{\,}\mathrm{nm}$$, as expected. The increase in state size during
evolution from Fig. 3(b) to Fig. 3(d) is a signature of heating due to
collisions of the particle with gas molecules [33, 56, 57].

((a))

((b))

((c))

Figure 4:  State size $\Delta u$ as a function of evolution time $t_{d}$, for
different Paul trap stiffnesses, quantified by the frequency ratio $r$. Solid
lines: full model described in the Supplemental Material [54]. Dashed lines:
simplified model discussed in the main text. Dotted lines: heating
contribution $\Delta u_{h}$ under the simplified model.

We proceed by investigating the role of the frequency ratio
$r=\omega_{o}/\omega_{p}$ (tuned via the voltage of the Paul trap RF drive) in
our expansion protocol. Figure 4(a) shows the measured state size $\Delta u$
(blue circles) as a function of evolution time $t_{d}$, for a frequency ratio
$r=8.8$. We identify the following three features. First, the state size grows
and shrinks with a period $T_{p}/2$, as expected. Second, we observe a fast
modulation whose period matches that of the RF-field driving the Paul trap.
Third, the minimum of each expansion cycle trends upwards, suggesting that the
minimal state size continuously grows as a consequence of heating during the
expansion protocol.

We investigate these features for larger frequency ratios $r=14.5$ and
$r=24.3$ in Figs. 4(b) and 4(c), respectively. For weaker Paul traps, the main
oscillations are slower and of larger amplitude than for stiffer traps. In
Fig. 4(c), we achieve a state expansion by a factor 25. We observe that weaker
Paul traps are less heavily modulated by oscillations at the RF frequency,
which we attribute to the lower RF drive voltages. Finally, weaker traps are
more sensitive to heating, seen by the more rapid increase of the minima in
Fig. 4(c) as compared to Fig. 4(a).

### Discussion.

To understand our observations, we model the behavior of a particle’s state
subject to a frequency-jump protocol. To simplify the description, we assume
that the Paul trap is a harmonic potential. A complete model including the
effects of micromotion can be found in Ref. [54]. We model the particle in the
Paul trap as a harmonic oscillator of mass $m$ coupled to a thermal bath at
rate $\gamma$, following the equation of motion

$\ddot{u}+\omega_{p}^{2}u=\sqrt{2\gamma k_{\text{B}}T/m}~{}\xi_{\text{th}},$
(1)

where we neglect the damping term (since $\gamma\ll\omega_{p}$), we assume
$u^{2}\ll{k_{B}T/(m\omega_{p}^{2})}$, and each dot represents a time
derivative. The bath has temperature $T$ and is modeled using the stochastic
variable $\xi_{\text{th}}$, which satisfies
$\left\langle\xi_{\text{th}}(t)\xi_{\text{th}}(t^{\prime})\right\rangle=\delta(t-t^{\prime})$.
Here, $\left\langle\dots\right\rangle$ denotes an ensemble average. Solving
Eq. 1, we find that the spatial state size of the particle, given by the
standard deviation $\Delta u$ after the protocol, evolves according to $\Delta
u^{2}(t)=\Delta u_{c}^{2}(t)+\Delta u_{h}^{2}(t)$. The first contribution
stems from coherent evolution and is given by [44]

$\Delta u_{c}^{2}(t)=\Delta
u_{0}^{2}\cos^{2}(\omega_{p}t)+r^{2}\frac{\Delta\dot{u}_{0}^{2}}{\omega_{o}^{2}}\sin^{2}(\omega_{p}t),$
(2)

with $\Delta u_{0}$ the initial state size. The second term of Eq. (2) scales
with the ratio of the initial velocity state size $\Delta\dot{u}_{0}^{2}$ and
the (square of the) frequency of the Paul trap
$\omega_{p}^{2}=\omega_{o}^{2}/r^{2}$. When the protocol is started with a
thermal state of the optical potential, we have
$\Delta\dot{u}_{0}/\omega_{o}=\Delta u_{0}$ and the maximum spatial state size
from coherent evolution is given by $\Delta u_{c,\text{max}}=r\Delta u_{0}$,
such that the spatial expansion is simply dictated by $r$. The situation in
our experiment is slightly more complicated. The particle is initialized in
the optical trap in two uncorrelated thermal states along $x$ and $y$.
However, we perform our subsequent experiments in the rotated basis $uv$. The
distributions along $u$ and $v$ then become correlated, and the relation
$\Delta\dot{u}_{0}/\omega_{o}=\Delta u_{0}$ no longer holds. For more detail,
refer to the supplemental material [54]. The second contribution to $\Delta
u^{2}$ arises from heating and reads

$\Delta
u_{h}^{2}(t)=r^{2}\frac{\hbar\Gamma}{m\omega_{o}}\left(t-\frac{\sin(2\omega_{p}t)}{2\omega_{p}}\right).$
(3)

Here, to align with the literature [13], we have introduced the heating rate
$\Gamma$, which is expressed in units of phonons of the optical trapping
potential. In our experiment, where heating is dominated by coupling to the
residual gas in the vacuum chamber at temperature $T$, we have $\Gamma=\gamma
k_{B}T/(\hbar\omega_{o})$, which we measure to be $\Gamma=2\pi\times$926\pm
56\text{\,}\mathrm{kHz}$$.

Care must be taken when defining the optical trap frequency $\omega_{o}$ in
Eq. (2). Since we are interested in the particle motion along $u$, but the
optical trap’s eigenmodes are oriented along $x$, $y$, and $z$, we define the
effective oscillation frequency along $u$ as
$\omega_{u\text{,opt}}={[(\omega_{x\text{,opt}}^{2}+\omega_{y\text{,opt}}^{2})/2}]^{1/2}$.
We fit our model to the measured data using the Paul trap frequency
$\omega_{p}$ as the only fit parameter. All other parameters ($\Delta u_{0}$,
$\Delta\dot{u}_{0}$, $\omega_{o}$, $\gamma$, and $T$) are measured
independently.

The simplified model is shown in Fig. 4 as the dashed gray lines, while the
dotted gray lines show the heating contribution $\Delta u_{h}$. Our simplified
model reproduces all features of the data well, except the modulation at the
micromotion frequency. In contrast, the full model (solid line in Fig. 4)
accounts for all features observed during a frequency-jump protocol. Comparing
the two models in Fig. 4(a) reveals that, for stiff Paul traps, the presence
of the RF fields almost doubles the expanded state size as compared to a
simple frequency jump.

Another important observation is that our simple model explains the increase
of the heating contribution with increasing frequency ratio $r$,
experimentally observed in Fig. 4. Specifically, according to Eq. 3, at
compression time $T_{p}/2$, the total state size has no contribution due to
coherent evolution, and is given by

$\Delta u^{2}(T_{p}/2)=\Delta
u_{0}^{2}+r^{2}\frac{\hbar\Gamma}{m\omega_{o}}\frac{T_{p}}{2}.$ (4)

The only contribution to $\Delta u^{2}$ besides the initial value $\Delta
u_{0}^{2}$ is from heating, with a rate amplified by $r^{2}$.

Regarding the current limitations of our system, larger expansions could, in
principle, be achieved with the help of weaker potentials. However, the
associated increase in sensitivity to stray fields leads to frequent failures
to recapture the particle in the optical trap and must be overcome with more
accurate stray field and gravity compensation in the future.

### Conclusion.

We have constructed a hybrid Paul-optical trap with high-numerical-aperture
optical access, as used in systems achieving measurement-based ground-state
cooling of optically levitated nanoparticles [16, 15]. Using this platform, we
have exploited the large frequency difference between optical and Paul traps
to expand and subsequently recontract the thermal state size of an optically
levitated nanoparticle by a factor of 25. The expansion happens in the absence
of photon recoil, which is critical for extending this protocol to generate
macroscopic quantum states.

Having demonstrated expansion of a thermal state, a logical next step is to
move to ultra-high vacuum, unlocking access to ground-state cooling [15, 16,
17]. In the quantum regime, an exciting prospect is to observe quantum
squeezing of mechanical motion in the quadrature orthogonal to the one that is
coherently expanded. Indeed, in the absence of decoherence an expansion factor
of 25 would generate 28 dB of motional squeezing in variance. A challenge in
this endeavor is the increased sensitivity to heating during the protocol. On
the other hand, this increased sensitivity can potentially be exploited. Since
state expansion in the Paul trap happens in the absence of laser light, our
system may provide access to weak decoherence rates, such as from blackbody
emission, which are overwhelmed by photon recoil in optically levitated
systems [6, 8, 58]. Finally, our platform could also be used to generate
inverted quadratic potentials, by reversing the endcap voltages. This may
provide an avenue to reach large expansions in short times [59, 60, 13].

###### Acknowledgements.

This research was supported by the Swiss National Science Foundation (SNF)
through Grant No. 212599 and No. 217122, NCCR-QSIT program (Grant No.
51NF40-160591), the European Union’s Horizon 2020 research and innovation
program under Grant No. 951234 (Q-Xtreme), and by the Austrian Science Fund
(FWF) Project No. Y951 and Grant No. I5540.

## References

  * Arndt and Hornberger [2014] M. Arndt and K. Hornberger, Testing the limits of quantum mechanical superpositions, Nature Physics 10, 271 (2014).
  * Fein  _et al._ [2019] Y. Y. Fein, P. Geyer, P. Zwick, F. Kiałka, S. Pedalino, M. Mayor, S. Gerlich, and M. Arndt, Quantum superposition of molecules beyond 25 kDa, Nature Physics 15, 1242 (2019).
  * Millen  _et al._ [2020] J. Millen, T. S. Monteiro, R. Pettit, and A. N. Vamivakas, Optomechanics with levitated particles, Reports on Progress in Physics 83, 026401 (2020).
  * Gonzalez-Ballestero  _et al._ [2021] C. Gonzalez-Ballestero, M. Aspelmeyer, L. Novotny, R. Quidant, and O. Romero-Isart, Levitodynamics: Levitation and control of microscopic objects in vacuum, Science 374, eabg3027 (2021).
  * Romero-Isart  _et al._ [2011] O. Romero-Isart, A. C. Pflanzer, F. Blaser, R. Kaltenbaek, N. Kiesel, M. Aspelmeyer, and J. I. Cirac, Large Quantum Superpositions and Interference of Massive Nanometer-Sized Objects, Physical Review Letters 107, 020405 (2011).
  * Romero-Isart [2011] O. Romero-Isart, Quantum superposition of massive objects and collapse models, Physical Review A 84, 052121 (2011).
  * Scala  _et al._ [2013] M. Scala, M. S. Kim, G. W. Morley, P. F. Barker, and S. Bose, Matter-Wave Interferometry of a Levitated Thermal Nano-Oscillator Induced and Probed by a Spin, Physical Review Letters 111, 180403 (2013).
  * Bateman  _et al._ [2014] J. Bateman, S. Nimmrichter, K. Hornberger, and H. Ulbricht, Near-field interferometry of a free-falling nanoparticle from a point-like source, Nature Communications 5, 4788 (2014).
  * Wan  _et al._ [2016] C. Wan, M. Scala, G. W. Morley, ATM. A. Rahman, H. Ulbricht, J. Bateman, P. F. Barker, S. Bose, and M. S. Kim, Free Nano-Object Ramsey Interferometry for Large Quantum Superpositions, Physical Review Letters 117, 143003 (2016).
  * Pino  _et al._ [2018] H. Pino, J. Prat-Camps, K. Sinha, B. P. Venkatesh, and O. Romero-Isart, On-chip quantum interference of a superconducting microsphere, Quantum Science and Technology 3, 025001 (2018).
  * Stickler  _et al._ [2018] B. A. Stickler, B. Papendell, S. Kuhn, B. Schrinski, J. Millen, M. Arndt, and K. Hornberger, Probing macroscopic quantum superpositions with nanorotors, New Journal of Physics 20, 122001 (2018).
  * Neumeier  _et al._ [2022] L. Neumeier, M. A. Ciampini, O. Romero-Isart, M. Aspelmeyer, and N. Kiesel, Fast quantum interference of a nanoparticle via optical potential control (2022), arXiv:2207.12539 .
  * Roda-Llordes  _et al._ [2023] M. Roda-Llordes, A. Riera-Campeny, D. Candoli, P. T. Grochowski, and O. Romero-Isart, Macroscopic quantum superpositions via dynamics in a wide double-well potential (2023), arXiv:2303.07959 .
  * Delić  _et al._ [2020] U. Delić, M. Reisenbauer, K. Dare, D. Grass, V. Vuletić, N. Kiesel, and M. Aspelmeyer, Cooling of a levitated nanoparticle to the motional quantum ground state, Science 367, 892 (2020).
  * Magrini  _et al._ [2021] L. Magrini, P. Rosenzweig, C. Bach, A. Deutschmann-Olek, S. G. Hofer, S. Hong, N. Kiesel, A. Kugi, and M. Aspelmeyer, Real-time optimal quantum control of mechanical motion at room temperature, Nature 595, 373 (2021).
  * Tebbenjohanns  _et al._ [2021] F. Tebbenjohanns, M. L. Mattana, M. Rossi, M. Frimmer, and L. Novotny, Quantum control of a nanoparticle optically levitated in cryogenic free space, Nature 595, 378 (2021).
  * Kamba  _et al._ [2022] M. Kamba, R. Shimizu, and K. Aikawa, Optical cold damping of neutral nanoparticles near the ground state in an optical lattice, Optics Express 30, 26716 (2022).
  * Ranfagni  _et al._ [2022] A. Ranfagni, K. Børkje, F. Marino, and F. Marin, Two-dimensional quantum motion of a levitated nanosphere, Phys. Rev. Res. 4, 033051 (2022).
  * Piotrowski  _et al._ [2023] J. Piotrowski, D. Windey, J. Vijayan, C. Gonzalez-Ballestero, A. de los Ríos Sommer, N. Meyer, R. Quidant, O. Romero-Isart, R. Reimann, and L. Novotny, Simultaneous ground-state cooling of two mechanical modes of a levitated nanoparticle, Nature Physics 19, 1009 (2023).
  * Rashid  _et al._ [2016] M. Rashid, T. Tufarelli, J. Bateman, J. Vovrosh, D. Hempston, M. S. Kim, and H. Ulbricht, Experimental Realization of a Thermal Squeezed State of Levitated Optomechanics, Physical Review Letters 117, 273601 (2016).
  * Hebestreit  _et al._ [2018a] E. Hebestreit, M. Frimmer, R. Reimann, and L. Novotny, Sensing Static Forces with Free-Falling Nanoparticles, Physical Review Letters 121, 063602 (2018a).
  * Frimmer  _et al._ [2019] M. Frimmer, T. L. Heugel, Ž. Nosan, F. Tebbenjohanns, D. Hälg, A. Akin, C. L. Degen, L. Novotny, R. Chitra, O. Zilberberg, and A. Eichler, Rapid Flipping of Parametric Phase States, Physical Review Letters 123, 254102 (2019).
  * Ciampini  _et al._ [2021] M. A. Ciampini, T. Wenzl, M. Konopik, G. Thalhammer, M. Aspelmeyer, E. Lutz, and N. Kiesel, Experimental nonequilibrium memory erasure beyond Landauer’s bound (2021), arxiv:2107.04429 [cond-mat, physics:physics, physics:quant-ph] .
  * Aspelmeyer  _et al._ [2014] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Cavity optomechanics, Reviews of Modern Physics 86, 1391 (2014).
  * Rossi  _et al._ [2018] M. Rossi, D. Mason, J. Chen, Y. Tsaturyan, and A. Schliesser, Measurement-based quantum control of mechanical motion, Nature 563, 53 (2018).
  * O’Connell  _et al._ [2010] A. D. O’Connell, M. Hofheinz, M. Ansmann, R. C. Bialczak, M. Lenander, E. Lucero, M. Neeley, D. Sank, H. Wang, M. Weides, J. Wenner, J. M. Martinis, and A. N. Cleland, Quantum ground state and single-phonon control of a mechanical resonator, Nature 464, 697 (2010).
  * Lecocq  _et al._ [2015] F. Lecocq, J. B. Clark, R. W. Simmonds, J. Aumentado, and J. D. Teufel, Quantum Nondemolition Measurement of a Nonclassical State of a Massive Object, Physical Review X 5, 041037 (2015).
  * Riedinger  _et al._ [2016] R. Riedinger, S. Hong, R. A. Norte, J. A. Slater, J. Shang, A. G. Krause, V. Anant, M. Aspelmeyer, and S. Gröblacher, Non-classical correlations between single photons and phonons from a mechanical oscillator, Nature 530, 313 (2016).
  * Møller  _et al._ [2017] C. B. Møller, R. A. Thomas, G. Vasilakis, E. Zeuthen, Y. Tsaturyan, M. Balabas, K. Jensen, A. Schliesser, K. Hammerer, and E. S. Polzik, Quantum back-action-evading measurement of motion in a negative mass reference frame, Nature 547, 191 (2017).
  * Chu  _et al._ [2018] Y. Chu, P. Kharel, T. Yoon, L. Frunzio, P. T. Rakich, and R. J. Schoelkopf, Creation and control of multi-phonon Fock states in a bulk acoustic-wave resonator, Nature 563, 666 (2018).
  * Sletten  _et al._ [2019] L. R. Sletten, B. A. Moores, J. J. Viennot, and K. W. Lehnert, Resolving Phonon Fock States in a Multimode Cavity with a Double-Slit Qubit, Physical Review X 9, 021056 (2019).
  * Bild  _et al._ [2023] M. Bild, M. Fadel, Y. Yang, U. von Lüpke, P. Martin, A. Bruno, and Y. Chu, Schrödinger cat states of a 16-microgram mechanical oscillator, Science 380, 274 (2023).
  * Jain  _et al._ [2016] V. Jain, J. Gieseler, C. Moritz, C. Dellago, R. Quidant, and L. Novotny, Direct Measurement of Photon Recoil from a Levitated Nanoparticle, Physical Review Letters 116, 243601 (2016).
  * Maurer  _et al._ [2023] P. Maurer, C. Gonzalez-Ballestero, and O. Romero-Isart, Quantum theory of light interaction with a lorenz-mie particle: Optical detection and three-dimensional ground-state cooling, Phys. Rev. A 108, 033714 (2023).
  * Dania  _et al._ [2023] L. Dania, D. S. Bykov, F. Goschin, M. Teller, and T. E. Northup, Ultra-high quality factor of a levitated nanomechanical oscillator (2023), arxiv:2304.02408 [quant-ph] .
  * Millen  _et al._ [2015] J. Millen, P. Z. G. Fonseca, T. Mavrogordatos, T. S. Monteiro, and P. F. Barker, Cavity Cooling a Single Charged Levitated Nanosphere, Physical Review Letters 114, 123602 (2015).
  * Nagornykh  _et al._ [2015] P. Nagornykh, J. E. Coppock, and B. E. Kane, Cooling of levitated graphene nanoplatelets in high vacuum, Applied Physics Letters 106, 244102 (2015).
  * Fonseca  _et al._ [2016] P. Z. G. Fonseca, E. B. Aranas, J. Millen, T. S. Monteiro, and P. F. Barker, Nonlinear Dynamics and Strong Cavity Cooling of Levitated Nanoparticles, Physical Review Letters 117, 173602 (2016).
  * Delord  _et al._ [2017] T. Delord, L. Nicolas, M. Bodini, and G. Hétet, Diamonds levitating in a Paul trap under vacuum: Measurements of laser-induced heating via NV center thermometry, Applied Physics Letters 111, 013101 (2017).
  * Conangla  _et al._ [2018] G. P. Conangla, A. W. Schell, R. A. Rica, and R. Quidant, Motion Control and Optical Interrogation of a Levitating Single Nitrogen Vacancy in Vacuum, Nano Letters 18, 3956 (2018).
  * Conangla  _et al._ [2020] G. P. Conangla, R. A. Rica, and R. Quidant, Extending Vacuum Trapping to Absorbing Objects with Hybrid Paul-Optical Traps, Nano Letters 20, 6018 (2020).
  * Bykov  _et al._ [2022] D. S. Bykov, M. Meusburger, L. Dania, and T. E. Northup, Hybrid electro-optical trap for experiments with levitated particles in vacuum, Review of Scientific Instruments 93, 073201 (2022).
  * Janszky and Yushin [1986] J. Janszky and Y. Yushin, Squeezing via frequency jump, Optics Communications 59, 151 (1986).
  * Graham [1987] R. Graham, Squeezing and Frequency Changes in Harmonic Oscillations, Journal of Modern Optics 34, 873 (1987).
  * Lo [1990] C. F. Lo, Squeezing by tuning the oscillator frequency, Journal of Physics A: Mathematical and General 23, 1155 (1990).
  * Cosco  _et al._ [2021] F. Cosco, J. S. Pedernales, and M. B. Plenio, Enhanced force sensitivity and entanglement in periodically driven optomechanics, Physical Review A 103, L061501 (2021).
  * Xin  _et al._ [2021] M. Xin, W. S. Leong, Z. Chen, Y. Wang, and S.-Y. Lan, Rapid quantum squeezing by jumping the harmonic oscillator frequency, Phys. Rev. Lett. 127, 183602 (2021).
  * Bonvin  _et al._ [2023] E. Bonvin, L. Devaud, M. Rossi, A. Militaru, L. Dania, D. S. Bykov, M. Teller, T. E. Northup, L. Novotny, and M. Frimmer, Hybrid paul-optical trap with large optical access for levitated optomechanics (2023), arXiv:2312.10131 .
  * Chen [2017] J.-S. Chen,  _Ticking near the Zero-Point Energy: Towards 1e-18 Accuracy in Al + Optical Clocks_, Ph.D. thesis, University of Colorado (2017).
  * Chen  _et al._ [2017] J.-S. Chen, S. M. Brewer, C. W. Chou, D. J. Wineland, D. R. Leibrandt, and D. B. Hume, Sympathetic ground state cooling and time-dilation shifts in an ${{}^{27}\mathrm{Al}}^{+}$ optical clock, Phys. Rev. Lett. 118, 053002 (2017).
  * Brewer  _et al._ [2019] S. M. Brewer, J.-S. Chen, A. M. Hankin, E. R. Clements, C. W. Chou, D. J. Wineland, D. B. Hume, and D. R. Leibrandt, ${{}^{27}{\mathrm{Al}}^{+}}$ quantum-logic clock with a systematic uncertainty below ${10}^{-18}$, Phys. Rev. Lett. 123, 033201 (2019).
  * Teller  _et al._ [2023] M. Teller, V. Messerer, K. Schüppert, Y. Zou, D. A. Fioretto, M. Galli, P. C. Holz, J. Reichel, and T. E. Northup, Integrating a fiber cavity into a wheel trap for strong ion–cavity coupling, AVS Quantum Science 5, 012001 (2023).
  * Bonvin [2023] E. Bonvin,  _Hybrid Traps for Coherent State Expansion of a Levitated Nanoparticle_ , Doctoral thesis, ETH Zurich, Zurich (2023).
  * Bonvin and et al. [2023] E. Bonvin and et al., Supplemental Material (2023).
  * Efron [1979] B. Efron, Bootstrap Methods: Another Look at the Jackknife, The Annals of Statistics 7, 10.1214/aos/1176344552 (1979).
  * Tebbenjohanns  _et al._ [2019] F. Tebbenjohanns, M. Frimmer, A. Militaru, V. Jain, and L. Novotny, Cold Damping of an Optically Levitated Nanoparticle to Microkelvin Temperatures, Physical Review Letters 122, 223601 (2019).
  * Conangla  _et al._ [2019] G. P. Conangla, F. Ricci, M. T. Cuairan, A. W. Schell, N. Meyer, and R. Quidant, Optimal Feedback Cooling of a Charged Levitated Nanoparticle with Adaptive Control, Physical Review Letters 122, 223602 (2019).
  * Agrenius  _et al._ [2023] T. Agrenius, C. Gonzalez-Ballestero, P. Maurer, and O. Romero-Isart, Interaction between an Optically Levitated Nanoparticle and Its Thermal Image: Internal Thermometry via Displacement Sensing, Physical Review Letters 130, 093601 (2023).
  * Romero-Isart [2017] O. Romero-Isart, Coherent inflation for large quantum superpositions of levitated microspheres, New Journal of Physics 19, 123029 (2017).
  * Weiss  _et al._ [2021] T. Weiss, M. Roda-Llordes, E. Torrontegui, M. Aspelmeyer, and O. Romero-Isart, Large Quantum Delocalization of a Levitated Nanoparticle Using Optimal Control: Applications for Force Sensing and Entangling via Weak Forces, Physical Review Letters 127, 023601 (2021).
  * Hauer  _et al._ [2013] B. Hauer, C. Doolin, K. Beach, and J. Davis, A general procedure for thermomechanical calibration of nano/micro-mechanical resonators, Annals of Physics 339, 181 (2013).
  * Hebestreit  _et al._ [2018b] E. Hebestreit, M. Frimmer, R. Reimann, C. Dellago, F. Ricci, and L. Novotny, Calibration and energy measurement of optically levitated nanoparticle sensors, Review of Scientific Instruments 89, 033111 (2018b).
  * Joos and Lindner [1989] E. Joos and A. Lindner, Langevin equation for the parametric oscillator: A model for ion confinement in a radio frequency trap, Zeitschrift für Physik D Atoms, Molecules and Clusters 11, 295 (1989).
  * McLachlan [1964] N. W. McLachlan,  _Theory and Application of Mathieu Functions_ (Dover Publications, New York, 1964).
  * Meixner and Schäfke [1954] J. Meixner and F. W. Schäfke,  _Mathieusche Funktionen und Sphäroidfunktionen_ (Springer Berlin Heidelberg, 1954).
  * Abramowitz and Stegun [1964] M. Abramowitz and IA. Stegun,  _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_ (Dover New York, 1964).

Supplemental Materials

##  S1 Change of basis

In this section, we analyze the state-space distribution of a harmonic
oscillator under a coordinate transformation. Our goal is to find a relation
between the state size in the coordinate system $uv$ of the Paul trap and that
in the $xy$ system of the optical trap.

We consider two orthogonal modes $x$ and $y$ of a particle trapped in a
harmonic potential. We express the velocities of the particle as $\dot{x}$ and
$\dot{y}$, and the modes have resonance frequencies $\omega_{x}$ and
$\omega_{y}$, respectively. The particle is undergoing linear feedback cooling
along both axes, such that we consider two effective temperatures $T_{x}$ and
$T_{y}$. Each mode is then in a thermal state, with uncorrelated positions and
velocities, and the phase space probability distribution $\rho$ is given by

$\rho(x,y,\dot{x},\dot{y})=\rho_{q}(x,y)\rho_{p}(\dot{x},\dot{y}),$ (S1)

with

$\displaystyle\rho_{q}(x,y)$ $\displaystyle=\frac{m\omega_{x}\omega_{y}}{2\pi
k_{\text{B}}\sqrt{T_{x}T_{y}}}\text{exp}\left[{-m\left(\frac{\omega_{x}^{2}x^{2}}{2k_{\text{B}}T_{x}}+\frac{\omega_{y}^{2}y^{2}}{2k_{\text{B}}T_{y}}\right)}\right],$
(S2a) $\displaystyle\rho_{p}(\dot{x},\dot{y})$ $\displaystyle=\frac{m}{2\pi
k_{\text{B}}\sqrt{T_{x}T_{y}}}\text{exp}\left[{-m\left(\frac{\dot{x}^{2}}{2k_{\text{B}}T_{x}}+\frac{\dot{y}^{2}}{2k_{\text{B}}T_{y}}\right)}\right],$
(S2b)

where we have separated the total probability distribution into its position
contribution $\rho_{q}(x,y)$ and velocity contribution
$\rho_{p}(\dot{x},\dot{y})$. Both distributions in Eqs. (S2) have diagonal
covariance matrices, which can be seen from the lack of cross terms.

We now introduce two orthogonal directions $u$ and $v$ which lie in the
$xy$-plane and are tilted at $45^{\circ}$ with respect to $x$ and $y$. We
perform a change of basis from $xy$ to $uv$ by applying

$x=\frac{u+v}{\sqrt{2}},\qquad
y=\frac{u-v}{\sqrt{2}},\qquad\dot{x}=\frac{\dot{u}+\dot{v}}{\sqrt{2}},\qquad\dot{y}=\frac{\dot{u}-\dot{v}}{\sqrt{2}}.$
(S3)

We can perform this transformation on $\rho_{q}(x,y)$ and
$\rho_{p}(\dot{x},\dot{y})$ independently. We have

$\displaystyle\tilde{\rho}_{q}(u,v)$
$\displaystyle=\frac{m\omega_{x}\omega_{y}}{2\pi
k_{\text{B}}\sqrt{T_{x}T_{y}}}\text{exp}\left[{-\frac{1}{2}\frac{m}{2k_{\text{B}}}\begin{pmatrix}{u+v}&{u-v}\end{pmatrix}\begin{pmatrix}\omega_{x}^{2}/T_{x}&0\\\
0&\omega_{y}^{2}/T_{y}\end{pmatrix}\begin{pmatrix}{u+v}\\\
{u-v}\end{pmatrix}}\right]$ (S4a)
$\displaystyle=\frac{m\omega_{x}\omega_{y}}{2\pi
k_{\text{B}}\sqrt{T_{x}T_{y}}}\text{exp}\left[{-\frac{1}{2}\begin{pmatrix}u&v\end{pmatrix}\underbrace{\frac{m}{2k_{\text{B}}}\begin{pmatrix}\omega_{x}^{2}/T_{x}+\omega_{y}^{2}/T_{y}&\omega_{x}^{2}/T_{x}-\omega_{y}^{2}/T_{y}\\\
\omega_{x}^{2}/T_{x}-\omega_{y}^{2}/T_{y}&\omega_{x}^{2}/T_{x}+\omega_{y}^{2}/T_{y}\end{pmatrix}}_{\Sigma_{q}^{-1}}\begin{pmatrix}u\\\
v\end{pmatrix}}\right],$

where we explicitly highlight the covariance matrix $\Sigma_{q}$ in the new
basis. With analogous steps, we derive

$\tilde{\rho}_{p}(\dot{u},\dot{v})=\frac{m}{2\pi
k_{\text{B}}\sqrt{T_{x}T_{y}}}\text{exp}\left[{-\frac{1}{2}\begin{pmatrix}\dot{u}&\dot{v}\end{pmatrix}\underbrace{\frac{m}{2k_{\text{B}}}\begin{pmatrix}1/T_{x}+1/T_{y}&1/T_{x}-1/T_{y}\\\
1/T_{x}-1/T_{y}&1/T_{x}+1/T_{y}\end{pmatrix}}_{\Sigma_{p}^{-1}}\begin{pmatrix}\dot{u}\\\
\dot{v}\end{pmatrix}}\right],$ (S4b)

from which the analogous definition of the new covariance matrix $\Sigma_{p}$
for the velocities follows. Finally, we find the expressions for the new
covariance matrices:

$\displaystyle\Sigma_{q}$
$\displaystyle=\frac{k_{\text{B}}}{2m}\begin{pmatrix}T_{x}/\omega_{x}^{2}+T_{y}/\omega_{y}^{2}&&T_{x}/\omega_{x}^{2}-T_{y}/\omega_{y}^{2}\\\
T_{x}/\omega_{x}^{2}-T_{y}/\omega_{y}^{2}&&T_{x}/\omega_{x}^{2}+T_{y}/\omega_{y}^{2}\end{pmatrix},$
(S5a) $\displaystyle\Sigma_{p}$
$\displaystyle=\frac{k_{\text{B}}}{2m}\begin{pmatrix}T_{x}+T_{y}&&T_{x}-T_{y}\\\
T_{x}-T_{y}&&T_{x}+T_{y}\end{pmatrix}.$ (S5b)

We can appreciate from Eqs. (S4) that the original Gaussian distribution
remains Gaussian in the new frame of reference. However, we now observe
correlations between $u$ and $v$, as well as between $\dot{u}$ and $\dot{v}$,
whose magnitude depends on the difference between the temperatures and
resonance frequencies. Interestingly, if the original axes have different
resonance frequencies $\omega_{x}$ and $\omega_{y}$, it is not possible to set
all cross correlation terms to zero simultaneously in Eqs. (S5): either the
new positions or the new velocities will remain correlated. The position and
velocity variances in the new basis can directly be read off from the
covariance matrix:

$\Delta u^{2}=\Delta
v^{2}=\frac{k_{\text{B}}}{2m}\left(\frac{T_{x}}{\omega_{x}^{2}}+\frac{T_{y}}{\omega_{y}^{2}}\right),\qquad\Delta\dot{u}^{2}=\Delta\dot{v}^{2}=\frac{k_{\text{B}}}{2m}\left({T_{x}}+{T_{y}}\right).$
(S6)

From these relations, we conclude that the initial state does not satisfy
$\Delta\dot{u}_{0}/\omega_{o}=\Delta u_{0}$. Instead, for our system
parameters, we find $\Delta\dot{u}_{0}/\omega_{o}\simeq 0.85\Delta u_{0}$,
indicating that the maximum achievable state size is $\Delta
u_{c,\text{max}}=0.85r\Delta u_{0}$.

##  S2 Position and velocity measurement

### Detector calibration.

The particle’s position is recorded in the form of electric signals generated
by a quadrant photo-diode (QPD). We convert these signals into length units
for a more meaningful quantitative analysis by following the process outlined
in [61, 62]. A measurement record of the particle in the optical trap is
acquired at higher pressures (typically $5\text{\,}\mathrm{mbar}$), under
which conditions the particle is in thermal equilibrium with a bath at ambient
temperature. This knowledge allows us to convert the measured signal (in
volts) to the expected particle displacement (in meters). Calibration is
performed separately for the $x$ and $y$ modes of the particle’s motion.

### Processing of frequency-jump protocol measurement records.

Next, let us describe the data analysis process used to extract the position
and velocity of a particle after a frequency-jump protocol. Each measurement
record consists of two time traces, measuring the particle position along $x$
and $y$, respectively. For each individual time trace, we perform the
following three steps. First we truncate the time trace, only keeping the data
$120\text{\,}\mathrm{\SIUnitSymbolMicro s}$ after particle recapture, i.e.,
during phase (v) of the protocol. The $100\text{\,}\mathrm{\SIUnitSymbolMicro
s}$ buffer serves to reject effects related to the AC-coupling of our data
acquisition system from the data analysis. Second, we demodulate this data at
the optical trap frequency. Third, we adjust the phase of the demodulated
components to account for timing delays in the system. The in-phase and
quadrature components of the final demodulated signal describe the particle’s
position and velocity upon recapture, respectively.

The filter stage of the demodulation process has a $2\text{\,}\mathrm{kHz}$
bandwidth and is applied backwards in time, effectively retrodicting the
signal. Special care must be taken when choosing the demodulator frequency.
Across a series of experiments, the laser power, and hence the trap frequency,
is subject to slight fluctuations. For each realization, we adjust the
demodulator frequency to precisely match that of the measured signal,
preventing the filtering stage from introducing unwanted phase shifts.

Any timing delays during the measurement sequence manifest themselves as
angled phase space distributions. These have several different origins. First,
the trigger that launches the frequency-jump protocol is not precisely aligned
to the sampling grid of the data acquisition system. This introduces a delay
of $02\text{\,}\mathrm{\SIUnitSymbolMicro s}$, which we record individually
for each iteration. Next, the response time of the AOM is limited by the
velocity of the acoustic wave in the AOM crystal. This introduces a fixed
delay measured to be $0.6\text{\,}\mathrm{\SIUnitSymbolMicro s}$. Finally, we
observe an additional delay of $0.9\text{\,}\mathrm{\SIUnitSymbolMicro s}$ of
unknown origin. An example of the same phase space before and after accounting
for timing delays is shown in Fig. S1. The variable trigger delay causes the
phase space distribution to fan out, while the fixed delay adds an angle to
the distribution.

Figure S1: Example phase space distribution (a) before and (b) after delay
adjustment. The variable trigger delay causes the phase space distribution to
fan out, while the fixed delay adds an angle to the distribution.

### Characterization of measurement noise.

In order to characterize the measurement noise of our system, we perform an
identical data analysis process to the one described above, with the only
difference being a shifted spectral window during the demodulation phase.
Instead of demodulating our measured signals in a $2\text{\,}\mathrm{kHz}$
band around the frequency of our mode of interest, we shift this window to a
spectral region with no signal. This resulting phase space distribution is the
product of noise in our system. We refer to the position variance of this
distribution as our measurement noise, which we measure to be
$2.5\text{\,}\mathrm{n}\mathrm{m}^{2}$. This noise is then subtracted from the
measured variance of our particle during a frequency-jump protocol.

##  S3 Frequency-jump protocol in a Paul trap

In this section, we describe the state size of a particle undergoing a
frequency-jump protocol in a Paul trap. The derivation below closely follows
[63]. We assume a particle of mass $m$, with $n_{e}$ elementary charges, in a
Paul trap. It is coupled with a rate $\gamma$ to a bath at temperature $T$.
Neglecting friction, we assume that it follows the Mathieu equation [64]

$\frac{d^{2}u}{d\tau^{2}}+\left(a_{u}-2q_{u}\cos(2\tau)\right)u=\sqrt{\frac{2\gamma
k_{\text{B}}T}{m}}~{}\frac{4\xi_{\text{th}}(\tau)}{\Omega^{2}}$ (S7)

where $\Omega$ is the angular frequency of the RF driving tone, and we have
introduced the dimensionless time variable $\tau=\Omega t/2$. Here, $a_{u}$
and $q_{u}$ are the conventional Mathieu parameters along $u$, given by

$a_{u}=\frac{4n_{e}e}{m\Omega^{2}}\frac{\partial^{2}\phi_{\text{dc}}}{\partial
u^{2}},\quad
q_{u}=\frac{2n_{e}e}{m\Omega^{2}}\frac{\partial^{2}\phi_{\text{rf}}}{\partial
u^{2}},$ (S8)

with $\phi_{\text{dc}}$ and $\phi_{\text{rf}}$ the DC and RF potentials of the
Paul trap, respectively. Floquet theory [65, 66] provides the following
linearly independent solutions to the homogeneous equation of Eq. S7

$\displaystyle\lambda_{1}(\tau)$
$\displaystyle=\sum_{n=-\infty}^{\infty}C_{2n}\cos((2n+\beta)\tau),$ (S9)
$\displaystyle\lambda_{2}(\tau)$
$\displaystyle=\sum_{n=-\infty}^{\infty}C_{2n}\sin((2n+\beta)\tau),$

with $\beta$ the characteristic exponent of the Mathieu equation. From this,
we can describe the full solution to the equation of motion Eq. S7

$u(\tau)=\tilde{u}_{0}\lambda_{1}(\tau)+(\tilde{v}_{u,0}/{\omega_{p}})\lambda_{2}(\tau)+\int_{0}^{\tau}g(\tau,\tau^{\prime})\sqrt{\frac{2\gamma
k_{\text{B}}T}{m}}~{}\frac{4\xi_{\text{th}}(\tau^{\prime})}{\Omega^{2}}d\tau^{\prime}$
(S10)

where $\tilde{u}_{0}$ and $\tilde{v}_{u,0}$ are related to the initial
position and velocity by

$u_{0}=\tilde{u}_{0}\sum_{n}C_{2n},\quad
v_{u,0}=\tilde{v}_{u,0}\sum_{n}(2n+\beta)C_{2n}/\beta,$ (S11)

and $\omega_{p}=\beta\Omega/2$ is the Paul trap’s secular frequency. The
Green’s function $g(\tau,\tau^{\prime})$ has the form

$g(\tau,\tau^{\prime})=\frac{\left(\lambda_{1}(\tau^{\prime})\lambda_{2}(\tau)-\lambda_{1}(\tau)\lambda_{2}(\tau^{\prime})\right)}{W}$
(S12)

where $W=\lambda_{1}\lambda^{\prime}_{2}-\lambda^{\prime}_{1}\lambda_{2}$ is
the Wronskian of the system. Substituting Eq. S9, we get

$g(\tau,\tau^{\prime})=\frac{\sum_{n,n^{\prime}}C_{2n}C_{2n^{\prime}}\sin\left((2n+\beta)\tau-(2n^{\prime}+\beta)\tau^{\prime}\right)}{\sum_{n}(2n+\beta)C_{2n}^{2}}.$
(S13)

We are now in a position to evaluate the state size of a particle evolving in
a Paul trap for a duration $t$ using Eq. S10. We break it down into its
coherent and heating contributions according to $\Delta u^{2}(t)=\Delta
u_{c}^{2}(t)+\Delta u_{h}^{2}(t)$. The contribution from coherent evolution is
given by

$\Delta
u_{c}^{2}(t)=\Delta\tilde{u}_{0}^{2}\lambda^{2}_{1}(t)+r^{2}\frac{\Delta\tilde{u}_{0}^{2}}{\omega_{\text{o}}^{2}}\lambda_{2}^{2}(t)$
(S14)

and the contribution from heating is described by

$\Delta u_{h}^{2}(t)=r^{2}\frac{2\gamma
k_{\text{B}}T}{m\omega_{o}^{2}}\int_{0}^{t}g^{2}(t,t^{\prime})dt^{\prime}$
(S15)

where we have introduced the frequency ratio $r=\omega_{o}/\omega_{p}$.

Both the coherent and heating contributions resemble their counterparts in the
simplified model Eqs. 2 and 3, with the main difference being that the
trigonometric functions are replaced with infinite sums of trigonometric
functions, capturing the effect of micromotion. For the heating contribution
Eq. S15, an added consequence is that the integral of the Green’s function can
no longer be analytically evaluated. If we choose to apply the pseudo-
potential approximation to Eqs. S14 and S15 by keeping only the lowest order
term, i.e., taking $C_{2n}=\delta_{n0}$, we recover Eqs. 2 and 3. The full
model used in the main text is numerically evaluated by keeping terms up to
$n=\pm 1$, using the expression $C_{\pm 2}={-q_{u}}/{(4\pm 4\beta)}$ [66].


# Two Parallel PageRank Algorithms via Improving Forward Push

Qi Zhang  Rongxia Tang  Zhengan Yao  Jun Liang  Department of Mathematics, Sun
Yat-sen University, China  School of Software, South China Normal University,
China

###### Abstract

Initially used to rank web pages, PageRank has now been applied in many
fields. With the growing scale of graph, accelerating PageRank computing is
urged and designing parallel algorithm is a feasible solution. In this paper,
two parallel PageRank algorithms IFP1 and IFP2 are proposed via improving the
state-of-the-art Personalized PageRank algorithm, i.e., Forward Push.
Theoretical analysis indicates that, IFP1 can take advantage of the DAG
structure of the graph, where the dangling vertices improves the convergence
rate and the unreferenced vertices decreases the computation amount. As an
improvement of IFP1, IFP2 pushes mass to the dangling vertices only once but
rather many times, and thus decreases the computation amount further.
Experiments on six data sets illustrate that both IFP1 and IFP2 outperform
Power method, where IFP2 with 38 parallelism can be at most 50 times as fast
as the Power method.

###### keywords:

PageRank, Forward Push, Parallel

##  1 Introduction

S.Brin and L.Page[1, 2] proposed PageRank while dealing with the problem of
ranking web pages retrieved by Google. PageRank measures the importance of web
pages according to the network structure rather than contents. Generally, one
web page has higher PageRank value if it is linked by more web pages, or the
web pages with a link to it are with higher PageRank value themselves.
Nowadays, PageRank’s application goes far more beyond Internet[3]. We can find
it in many fields such as social network analysis, chemistry, molecular
biology, sports[4] and social sciences[5].

In the past decades, plenty of PageRank algorithms were proposed and among
which, Power method is the most basic one. Some improvements based on Power
method were presented. Kamvar[6] proposed the adaptive method which excluded
the vertices had converged from the computation immediately. Haveliwala[7] and
Kamvar[8] proposed the extrapolation method which focused on the second
largest eigenvalue and made the best use of the previous iteration result.
Kamvar[9] proposed the block method by utilizing the network’s block
structure. Gang[10] proposed the POWER-ARNOLDI method via introducing Arnoldi-
type algorithm into PageRank computation.

Monte Carlo(MC) method [11] is another important PageRank algorithm. MC method
simulates random walks on graph and approximates PageRank value by the
probability that random walks terminate at each vertex. Some improvements
based on MC method for undirected graph were proposed. Sarma[12] accelerated
the walk by stitching short paths into a long one. Luo[13, 14] proposed Radar
Push and thus obtained a result with lower variance.

In addition, distributed PageRank algorithms such as Sankaralingam[15],
Zhu[16], Stergiou[17], Ishii[18] and [19, 20, 21, 22, 23] were proposed. [24,
25, 26] took advantage of the DAG structure by linear algebra. [27, 28] tried
computing PageRank on GPU.

Algorithms mentioned above have respective advantages. With the explosive
growth in the scale of graph, accelerating PageRank computing is urged.
Parallelizing existing algorithms such as Power method[29] is a feasible
solution, however, there stand barriers as below.

  1. (1)

Algorithms based on the Power method converge slowly and can be partially
parallel only, since the dependency between different iterations can not be
eliminated.

  2. (2)

Algorithms based on the MC method perform poorly on real application, since
large memory space and bandwidth are required.

Forward Push, which is the state-of-the-art Personalized PageRank(PPR)
algorithm, has been attracting more attention recently. With initial mass
distribution $\bm{p}$, each vertex reserves $1-c$ proportion of the mass
receiving from its source vertices, then evenly pushes the remaining $c$
proportion to its target vertices. The mass each vertex obtains in the ending
is just the PPR value. Since each vertex executes reserving and pushing
operations independently, Forward Push is easy to parallelizing. Designing
parallel PageRank algorithm via improving Forward Push is an natural idea. We
can obtain the correct PageRank vector on strongly connected graph, however,
problems arise when there exist dangling vertices in the graph. Specifically,
based on the random walk model, PageRank requires random walk arriving at
dangling vertex randomly jump to any vertex, while Forward Push terminates
this walk. The slight difference leads that, utilizing Forward Push for
PageRank computation forcibly and counting PPR vector as PageRank vector
simply, we can either terminate pushing operation at dangling vertices and
thus obtain wrong result, or continue pushing mass at dangling vertices to
every vertex and thus generate plenty of computation.

In this paper, we propose parallel PageRank algorithms via improving Forward
Push which addressed the problems mentioned above. The contributions are as
follows.

  1. $(1)$

We reveal that PageRank vector is essentially the probability distribution of
mass and can be obtained by Forward Push which terminates pushing mass at
dangling vertices.

  2. $(2)$

We propose two parallel PageRank algorithms, IFP1 and IFP2, via improving
Forward Push. Compared with the Power method, both IFP1 and IFP2 can take
advantage of DAG structure of the graph, has higher convergence rate and
generates less computation amount.

  3. $(3)$

Experimental results on six data sets demonstrate that both IFP1 and IFP2
outperform the Power method, where IFP2 with 38 parallelism can be at most 50
times as fast as.

The remaining of this paper are as follows. In section 2, we introduce
PageRank and Forward Push. In section 3, we firstly give the solution of
computing PageRank via Forward Push which terminates pushing mass at dangling
vertices, then present IFP1 and corresponding theoretical analysis, and
propose IFP2 at last. Some numerical experiments are performed in section 4.
We summarize this paper in section 5.

##  2 Preliminary

In this section, PageRank and Forward Push will be briefly introduced. Then
the problem arose by dangling vertices when computing PageRank by Forward Push
will be detailed.

###  2.1 PageRank

Given graph $G(V,E)$, where $V=\\{v_{1},v_{2},\cdots,v_{n}\\}$,
$E=\\{(v_{i},v_{j}):i,j=1,2,\cdots,n\\}$ and $|E|=m$. Denote by
$\bm{A}=(a_{ij})_{n\times n}$ the adjacency matrix, where $a_{ij}$ is 1 if
$(v_{j},v_{i})\in E$ and 0 else. Denote by $\bm{P}=(p_{ij})_{n\times n}$ the
probability transition matrix, where

$p_{ij}=\left\\{\begin{array}[]{ll}a_{ij}/\sum\limits_{i=1}^{n}{a_{ij}},&\textrm{if
$\sum\limits_{i=1}^{n}{a_{ij}}\neq 0$},\\\
0,&\textrm{else}.\end{array}\right.$

We call vertex without in-link the unreferenced vertex, vertex without out-
link the dangling vertex. Delete unreferenced (dangling) vertex, the newly
generating unreferenced (dangling) vertex is called the weak unreferenced
(dangling) vertex. Let $\bm{d}=(d_{1},d_{2},...,d_{n})^{T}$, where $d_{i}$ is
1 if $v_{i}$ is dangling vertex and 0 else. Let
$\bm{p}=(p_{1},p_{2},\cdots,p_{n})^{T}$ denote the $n$-dimensional probability
distribution. Let $\bm{P}^{{}^{\prime}}=\bm{P}+\bm{p}\bm{d}^{T}$ and
$\bm{P}^{{}^{\prime\prime}}=c\bm{P}^{{}^{\prime}}+(1-c)\bm{p}\bm{e}^{T}$,
where $c\in(0,1)$ is damping factor and $\bm{e}=(1,1,...,1)^{T}$. Denote by
$\bm{\pi}$ the PageRank vector. According to[30], let
$\bm{p}=\frac{\bm{e}}{n}$, one of PageRank’s definitions is

$\bm{\pi}=\bm{P}^{{}^{\prime\prime}}\bm{\pi},\pi_{i}>0,\sum\limits_{i=1}^{n}\pi_{i}=1.$
(1)

PageRank can be interpreted from many perspectives. Based on random walk
model, a random walk starts at random vertex, with probability $c$ walks
according to the graph, and with probability $1-c$ terminates, then $\bm{\pi}$
is the probability distribution of random walk terminating at each vertex.
When $\bm{p}\neq\frac{\bm{e}}{n}$, PageRank converts to Personalized
PageRank(PPR), and $\bm{p}$ is called the personalized vector. PageRank is a
special case of PPR.

It should be noted that, the graph which rules the random walk is actually the
graph corresponding to $\bm{P}^{{}^{\prime}}$ but rather the original one. The
difference between these two is mainly on the dangling vertices, specifically,
the former artificially links each dangling vertex to every vertices while the
latter does nothing. As a consequence, when a random walk arriving at the
dangling vertices, the former requires with probability $c$ randomly choosing
a target vertex and continuing the random walk, while the latter means
terminating. Graph corresponding to $\bm{P}^{{}^{\prime}}$ contains more extra
edges. From the perspective of computing, the existing of dangling vertices
always generates more computation when walking according to
$\bm{P}^{{}^{\prime}}$, and that is why directly utilizing Forward Push in
computing PageRank is costly.

###  2.2 Forward Push

Forward Push is the state-of-the-art PPR algorithm. As described in Algorithm
1, with initial mass distribution $\bm{p}$, each vertex does the following
two: (1)reserve $1-c$ proportion of mass it has received; (2)evenly push the
remaining $c$ proportion to its target vertices. While Forward Push running,
the pushing mass, i.e., the mass needing to push to the target vertices
decreases and the reserved mass increases. The algorithm finishes when there’s
no vertex holds more than the pre-defined $\xi$ pushing mass. The reserved
mass of each vertex is just its PPR value.

1:

2:$G(V,E)$:The graph;

3:$c$:The damping factor;

4:$\bm{p}$:The personalized vector;

5:$\xi$:The tolerance of error.

6:

7:$\bm{\pi}$:Personalized PageRank vector.

8:

9:Each vertex $v_{i}$ maintains a data structure
$\langle\overline{\pi}_{i},h_{i}\rangle$.

10:Initially set $\overline{\pi}_{i}=0$, $h_{i}=p_{i}$.

11:

12:while There exists vertex $v_{i}$ satisfying $h_{i}>\xi$ do

13: $\overline{\pi}_{i}+=(1-c)h_{i}$;

14: for $v_{j}\in D(v_{i})$ do$\triangleright$ [$D(v_{i})$ are target vertices
of $v_{i}$.]

15: $h_{j}+=\frac{ch_{i}}{deg(v_{i})}$;

16: end for

17: $h_{i}=0$;

18:end while

19:$\pi_{i}=\overline{\pi}_{i}$.

Algorithm 1 Forward Push

Forward Push is convenient to parallelizing since each vertex executes
reserving and pushing operation independently. Previous works seldom detailed
how to address mass on the dangling vertices. As a special case of PPR,
PageRank can be obtained by Forward Push when the graph is strongly connected.
However, for graph containing dangling vertices, Forward Push can either

  1. (1)

terminate pushing the mass on dangling vertices, and thus obtain an incorrect
result;

  2. (2)

continue pushing mass on dangling vertices evenly to every vertex, and thus
generate plenty of computation.

Moreover, in real computing environment, if the scale of graph is large
enough, the mass each vertex getting from the dangling vertices will be so
small that the result may loss precision. Graphs abstracted from reality
always contain large proportion of dangling vertices. To design parallel
PageRank algorithm based on Forward Push, the problem arose by dangling
vertices needs to be addressed firstly.

##  3 Algorithm and convergence analysis

In this section, we firstly demonstrate that PageRank vector is essentially a
distribution of mass, and it is feasible to compute PageRank via Forward Push
that needs not push mass on dangling vertices. Then we propose the parallel
PageRank algorithm IFP1 by improving Forward Push and present theoretical
analysis as well. At last, we propose IFP2, the improvement of IFP1.

###  3.1 Computing PageRank vector via Forward Push

By expanding Formula (1), we have

$(\bm{I}-c\bm{P}^{{}^{\prime}})\bm{\pi}=(1-c)\bm{p}.$ (2)

Since $\rho(c\bm{P}^{{}^{\prime}})<1$,
$(\bm{I}-c\bm{P}^{{}^{\prime}})^{-1}=\sum\limits_{r=0}^{\infty}(c\bm{P}^{{}^{\prime}})^{r}$,
it follows that

$\bm{\pi}=(1-c)\sum\limits_{r=0}^{\infty}(c\bm{P}^{{}^{\prime}})^{r}\bm{p}.$
(3)

Formula (3) is just the algebraic form of Forward Push, the masses on dangling
vertices are pushed to every vertex according to $\bm{p}$.

By expanding (2), we have

$(\bm{I}-c\bm{P})\bm{\pi}=(c\bm{d^{T}\pi}+1-c)\bm{p}$.

Since $\bm{I}-c\bm{P}$ is invertible , it follows that

$\bm{\pi}=(c\bm{d^{T}\pi}+1-c)(\bm{I}-c\bm{P})^{-1}\bm{p}$.

Let $\gamma=c\bm{d^{T}\pi}+1-c$, it follows that

$\bm{\pi}=\gamma(\bm{I}-c\bm{P})^{-1}\bm{p}$.

Since $1=\bm{e^{T}\pi}=\gamma\bm{e^{T}}(\bm{I}-c\bm{P})^{-1}\bm{p}$, we have

$\gamma=\frac{1}{\bm{e^{T}}(\bm{I}-c\bm{P})^{-1}\bm{p}}$.

Since $\pi_{i}=\gamma\bm{e_{i}^{T}}(\bm{I}-c\bm{P})^{-1}\bm{p}$, where
$\bm{e}_{i}$ is $n$-dimensional vector with the $i_{th}$ element is 1 and 0
others, it follows that

$\pi_{i}=\frac{\bm{e_{i}^{T}}(\bm{I}-c\bm{P})^{-1}\bm{p}}{\bm{e^{T}}(\bm{I}-c\bm{P})^{-1}\bm{p}}$.

Since $(\bm{I}-c\bm{P})^{-1}=\sum\limits_{r=0}^{\infty}(c\bm{P})^{r}$, we have

$\pi_{i}=\frac{\bm{e_{i}^{T}}\sum\limits_{r=0}^{\infty}(c\bm{P})^{r}\bm{p}}{\bm{e^{T}}\sum\limits_{r=0}^{\infty}(c\bm{P})^{r}\bm{p}}.$
(4)

$\sum\limits_{r=0}^{\infty}(c\bm{P})^{r}\bm{p}$ is the algebraic form of
Forward Push as well, it is different from Formula (3) on two aspects: (1)it
terminates pushing mass on the dangling vertices; (2) it reserves 100% but not
$1-c$ proportion of the mass.

Formula (4) demonstrates that PageRank vector is essentially a distribution of
reserving mass. The proportion each vertex reserves has no effect on the final
result, however, if each vertex reserves 100% proportion of mass, there’s no
need for dangling vertices executing reserving operation. Moreover, it is
$\bm{P}$ but rather $\bm{P}^{{}^{\prime}}$ that ruled the process of Forward
Push, the mass on dangling vertices need not to be pushed any more. Formula
(4) implies a solution of addressing the problem arose by dangling vertices.

###  3.2 IFP1

The remaining issue is parallelizing. Restricted by the computing resource,
assigning exclusive thread for every vertex is infeasible. We can generate
some threads and assign vertices to them. Then IFP1 is proposed as Algorithm
2.

1:

2:$K$:The number of threads;

3:$\xi$:The lower bound of mass.

4:

5:$\bm{\pi}$:PageRank vector.

6:Each vertex $v_{i}$ maintains a data structure
$\langle\overline{\pi}_{i},h_{i}\rangle$.

7:Assign non-dangling vertices to $K$ threads, denote by $S_{j}$ the set of
vertices belonging to thread $j$.

8:Initially set $\overline{\pi}_{i}=0$, $h_{i}=1$.

9:Invoke $K$ Calculations and Management;$\triangleright$ [The $K$
Calculations and Management do in parallel.]

10:Calculate $\pi$ following
$\pi_{i}=\frac{\overline{\pi}_{i}+h_{i}}{\sum\limits_{i=1}^{n}(\overline{\pi}_{i}+h_{i})}$
while the Management terminates.

11:

12:function Calculation($j$)

13: while 1 do

14: for $v_{i}\in S_{j}$ do

15: if $h_{i}>\xi$ then

16: $\overline{\pi}_{i}=\overline{\pi}_{i}+h_{i}$;

17: for $u\in D(v_{i})$ do$\triangleright$ [$D(v_{i})$ is the set of target
vertices of $v_{i}$.]

18: $h_{u}=h_{u}+\frac{ch_{i}}{deg(v_{i})}$;

19: end for

20: $h_{i}=0;$

21: end if

22: end for

23: end while

24:end function

25:

26:function Management

27: while There exists non-dangling vertex satisfying $h_{i}>\xi$ do

28: end while

29: Terminate all the $K$ Calculations.

30:end function

Algorithm 2 IFP1

While IFP1 running, thread $j$ circularly scans $S_{j}$, any vertex $v_{i}$
satisfying $h_{i}>\xi$ will be processed. The mass received from its source
vertices are reserved by 100% proportion, and then pushed to its target
vertices by $c$ proportion. On the whole, the reserved mass increases and
pushing mass decreases. The reserved mass of unreferenced vertices and weak
unreferenced vertices stay steadily after several iterations, i.e., these
vertices get converged. If there exists no non-dangling vertex holds more than
the predefined $\xi$ pushing mass, IFP1 gets converged. IFP1 is similar to
Forward Push except the following:

  1. (1)

IFP1 is parallel while Forward Push is serial;

  2. (2)

IFP1 processes non-dangling vertices only while Forward Push addresses all of
them;

  3. (3)

IFP1 reverses 100% proportion of the mass while Forward Push reserves $1-c$
proportion.

It should be noted that, in multi-thread environment, the addition on $h_{i}$
must be atomic and thus data race may occur.

###  3.3 Algorithm analysis

We analyse IFP1 from three aspects, the convergence rate, error and
computation amount.

####  3.3.1 Convergence rate

Convergence rate relates to iteration rounds, however, as a parallel algorithm
IFP1 in fact has no iteration. We define one iteration of IFP1 as that the
threads finish scanning through the whole vertices. Since whether processes
the dangling vertices, the proportion of reserving mass and the parallelism
have no effect on the convergence rate, for convenience, we assume that
(1)IFP1 addresses all of the vertices; (2)the proportion of reserving mass is
still $1-c$; (3)the parallelism is 1.

Denote by $\overline{\bm{\pi}}^{I}(t)$ the mass having been reserved at the
beginning of the $t_{th}$ iteration, by $\overline{\bm{\pi}}^{R}(t)$ the mass
needing to push at the beginning of the $t_{th}$ iteration, then it follows
that

$||\overline{\bm{\pi}}^{I}(t)+\overline{\bm{\pi}}^{R}(t)||_{1}\leq n$,

where $\overline{\bm{\pi}}^{I}(0)=\bm{0}$,
$\overline{\bm{\pi}}^{R}(0)=\bm{e}$,
$\lim\limits_{t\to\infty}\frac{\overline{\bm{\pi}}^{I}(t)}{||\overline{\bm{\pi}}^{I}(t)||_{1}}=\bm{\pi}$,
$\lim\limits_{t\to\infty}\overline{\bm{\pi}}^{R}(t)=\bm{0}$,
$||\overline{\bm{\pi}}^{I}(t)||_{1}$ increases monotonously on $t$, and
$||\overline{\bm{\pi}}^{R}(t)||_{1}$ decreases monotonously on $t$. The
convergence rate can be measured by
$\frac{||\overline{\bm{\pi}}(t)||_{1}}{||\overline{\bm{\pi}}(t+1)||_{1}}$.

At the beginning of the $t_{th}$ iteration, some vertices having been
converged and denote them by $V_{U}(t)$. Denote by $V_{D}$ the set of dangling
vertices, it follows that

$||\overline{\bm{\pi}}^{R}(t)||_{1}=\sum\limits_{v_{i}\in
V_{U}(t)}\overline{\pi}^{R}_{i}(t)+\sum\limits_{v_{i}\in
V_{1}(t)}\overline{\pi}^{R}_{i}(t)+\sum\limits_{v_{i}\in
V_{2}(t)}\overline{\pi}^{R}_{i}(t)$,

where $V_{1}(t)=V_{D}-V_{U}(t)$ and $V_{2}(t)=V-V_{U}(t)-V_{D}$.

During the $t_{th}$ iteration, the reserved mass of $V_{U}(t)$ stay constant;
the mass needing to push of $V_{1}(t)$ decrease to 0; the sum of mass needing
to push of $V_{2}(t)$ decrease by $1-c$ proportion, thus it follows that

$||\overline{\bm{\pi}}^{R}(t+1)||_{1}=\sum\limits_{v_{i}\in
V_{U}(t)}\overline{\pi}^{R}_{i}(t)+c\sum\limits_{v_{i}\in
V_{2}(t)}\overline{\pi}^{R}_{i}(t)$.

Vertices belonging to $V_{U}(t)$ have been converged, we have

$\sum\limits_{v_{i}\in V_{U}(t)}\overline{\pi}^{R}_{i}(t)<|V_{U}(t)|\xi$.

Assuming that $\xi$ is sufficiently small that $|V_{U}(t)|\xi\to 0$, then

$\frac{||\overline{\bm{\pi}}^{R}(t+1)||_{1}}{||\overline{\bm{\pi}}^{R}(t)||_{1}}=\frac{c\sum_{v_{i}\in
V_{2}(t)}\overline{\pi}^{R}_{i}(t)}{\sum_{v_{i}\in
V_{1}(t)}\overline{\pi}^{R}_{i}(t)+\sum_{v_{i}\in
V_{2}(t)}\overline{\pi}^{R}_{i}(t)}$.

Let $\alpha(t)=\frac{\sum_{v_{i}\in
V_{2}(t)}\overline{\pi}^{R}_{i}(t)}{\sum_{v_{i}\in
V_{1}(t)}\overline{\pi}^{R}_{i}(t)+\sum_{v_{i}\in
V_{2}(t)}\overline{\pi}^{R}_{i}(t)}$, then

$\frac{||\overline{\bm{\pi}}^{R}(t+1)||_{1}}{||\overline{\bm{\pi}}^{R}(t)||_{1}}=c\alpha(t)$.

Since $V_{1}(t)\cup V_{2}=V-V_{U}(t)$, we have

$\alpha(t)=\frac{\sum_{v_{i}\in
V-V_{U}(t)-V_{D}}\overline{\pi}^{R}_{i}(t)}{\sum_{v_{i}\in
V-V_{U}(t)}\overline{\pi}^{R}_{i}(t)}$.

$\alpha(t)$ is the proportion of pushing mass on non-dangling vertices among
the total pushing mass at the beginning of the $t_{th}$ iteration.

Generally, $\alpha(t)$ negatively correlates with the proportion of dangling
vertices. While IFP1 running, the amount of non-converged vertices decreases,
then the proportion of dangling vertices among non-converged vertices
increases, and thus $\alpha(t)$ decreases. Let $\alpha=\max\limits_{t\geq
1}\left\\{\alpha(t)\right\\}$ and $\lambda=\alpha c$, then

$\min\limits_{t\geq
1}\left\\{\frac{||\overline{\bm{\pi}}^{R}(t)||_{1}}{||\overline{\bm{\pi}}^{R}(t+1)||_{1}}\right\\}=\frac{1}{\max\limits_{t\geq
1}\left\\{\alpha(t)\right\\}c}=\frac{1}{\alpha c}=\lambda^{-1}$.

$\lambda^{-1}$ can be viewed as IFP1’s convergence rate. Generally, the larger
proportion of dangling vertices a graph has, the higher convergence rate IFP1
is. The iterations IFP1 takes to get converged under the predefined threshold
$\xi$ can be estimated by

$||\overline{\bm{\pi}}^{R}(T)||_{1}=n\prod\limits_{t=0}^{T-1}c\alpha(t)=\sum\limits_{v_{i}\in
V_{1}(t)}{\overline{\pi}^{R}_{i}(T)}+\sum\limits_{v_{i}\in
V_{2}(t)}{\overline{\pi}^{R}_{i}(T)}$.

Since IFP1 gets converged when $||\bm{\pi}_{R}(t)||_{\infty}<\xi$ and mass on
$V_{1}(t)$ decrease to 0 during the $t_{th}$ iteration, one of the sufficient
conditions is

$n\prod\limits_{t=0}^{T-1}c\alpha(t)\leq
n(\lambda)^{T}<|V-V_{U}(t)-V_{D}|\xi$,

and thus

$T>\log_{\lambda}{\xi}+\log_{\lambda}{\frac{|V-V_{U}(T)-V_{D}|}{n}}$.

Generally, $|V-V_{U}(T)-V_{D}|=O(n)$, we have

$T=O(\log_{\lambda}{\xi}).$ (5)

####  3.3.2 Error

It is infeasible to estimate the error of each vertex since the complication
of graph structure. We discuss the effect of pre-defined threshold $\xi$ on
relative error $ERR(\xi)$ from the whole. Assuming
$||\overline{\bm{\pi}}^{R}(t)||_{1}=n\lambda^{t}$, then

$||\overline{\bm{\pi}}^{I}(T+t)-\overline{\bm{\pi}}^{I}(T)||_{1}\leq
n|\lambda^{T}-\lambda^{T+t}|=(1-\lambda^{t})n\lambda^{T}$.

Let $t\rightarrow\infty$, then we can obtain that

$ERR(\xi)=\lim\limits_{t\rightarrow\infty}\frac{||\bm{\pi}^{I}(T+t)-\bm{\pi}^{I}(T)||_{1}}{||\bm{\pi}^{I}(T+t)||_{1}}<\lim\limits_{t\rightarrow\infty}\frac{|\lambda^{T}-\lambda^{T+t}|}{1-\lambda^{T+t}}=\xi.$
(6)

Formula (6) can be viewed as an estimation of relative error under the
predefined threshold $\xi$. It is to be noted that the error analysis
mentioned above is rough, it considers no any factor of graph structure.

####  3.3.3 Computation amount

The computation of IFP1 mainly consists by additions and productions. Denote
by $m(t)$ the computation amount of the $t_{th}$ iteration, and by $M$ the
total computation amount for IFP1 getting converged. Since only vertices
belonging to $V-V_{U}(t)$ generate arithmetic operations during the $t_{th}$
iteration, it follows that

$m(t)=\sum_{v_{i}\in V-V_{U}(t)}(deg(v_{i})+1)$.

Let $\beta(t)=\frac{m(t+1)}{m(t)}$. $\beta(t)$ indicates the decreasing rate
of computation amount. Unreferenced vertices and weak unreferenced vertices
gradually exit computing while IFP1 running, thus $\beta(t)\leq 1$ and

$m+n=m(0)\geq m(1)\geq m(2)\geq\cdots\geq m(T)$.

Let $\beta=\max\limits_{1\leq t\leq T}\left\\{\beta(t)\right\\}$. Assuming
that $\beta<1$, i.e., there exist vertices getting converged at each
iteration, we have

$M=\sum\limits_{t=0}^{T}m(t)=\sum\limits_{t=0}^{T}\sum_{v_{i}\in
V-V_{U}(t)}(deg(v_{i})+1)<(m+n)\frac{1-\beta^{T}}{1-\beta}$,

where $T=O(\log_{\lambda}{\xi})$.

As a summary, it is clear that IFP1 could take advantage of DAG structure,
where the dangling vertices improve the convergence rate and the unreferenced
vertices lower the computation amount. Compared with Power method, on graph
containing DAG structure, IFP1 requires less iterations to get converged and
generates less computations at each iteration, thus is faster.

###  3.4 IFP2

IFP1 requires no the dangling vertices executing reserving operation, however,
that pushing mass to dangling vertex from its source vertices is still needed
and sometimes will be executed many times. Specifically, denote by $v_{d}$ a
dangling vertex and by $S(v_{d})$ the set of $v_{d}$’s source vertices, each
time $v_{i}$ satisfies $h_{i}>\xi$, the pushing operation from $v_{i}$ to
$v_{d}$ will be executed. On the other hand, we have

$\sum\limits_{t=0}^{T}\overline{\pi}_{d}^{R}(t)=\sum\limits_{v_{i}\in
S(v_{d})}\frac{c}{deg(v_{i})}\sum\limits_{t=0}^{T}\overline{\pi}_{i}^{R}(t)$,

where $\sum\limits_{t=0}^{T}\overline{\pi}_{i}^{R}(t)$ is the reserved mass of
$v_{i}$. That implies the mass on dangling vertices is completely determined
by their source vertices. If we do not push mass to the the dangling vertices
from their source vertices initially, but execute these operations after all
the non-dangling vertices getting converged, then only once pushing operation
is sufficient. Motivated by this, IFP2, an improvement of IFP1, is proposed as
Algorithm 3.

1:

2:$K$:The number of threads;

3:$\xi$:The lower bound of mass.

4:

5:$\bm{\pi}$:PageRank vector.

6:

7:Preprocess the graph data and invalid edges between dangling vertices and
their source vertices.

8:Each vertex $v_{i}$ maintains a data structure
$\langle\overline{\pi}_{i},h_{i}\rangle$.

9:Assign vertices to $K$ threads, denote by $S^{1}_{j}$ and $S^{2}_{j}$ the
set of non-dangling vertices and dangling vertices belonging to thread $j$
respectively.

10:Initially set $\overline{\pi}_{i}=0$, $h_{i}=1$.

11:Invoke $K$ Calculations and Management;$\triangleright$ [The $K$
Calculations and Management do in parallel.]

12:Calculate $\pi$ following
$\pi_{i}=\frac{\overline{\pi}_{i}}{\sum\limits_{i=1}^{n}\overline{\pi}_{i}}$
while the Management terminates.

13:

14:function Calculation($j$)

15: while 1 do

16: if $CTRL$ then $\triangleright$ [$CTRL$ is bool and initially true.]

17: Phase1:

18: for $v_{i}\in S^{1}_{j}$ do

19: if $h_{i}>\xi$ then

20: $\overline{\pi}_{i}=\overline{\pi}_{i}+h_{i}$;

21: for $u\in D(v_{i})$ do$\triangleright$ [$D(v_{i})$ is the set of target
vertices of $v_{i}$.]

22: $h_{u}=h_{u}+\frac{ch_{i}}{deg(v_{i})}$;

23: end for

24: $h_{i}=0;$

25: end if

26: end for

27: else

28: Phase2:

29: if $STS$ then $\triangleright$ [$STS$ is bool and initially true.]

30: for $v_{i}\in S^{2}_{j}$ do

31: $\overline{\pi}_{i}=h_{i}$;

32: for $u\in S(v_{i})$ do $\triangleright$ [$S(v_{i})$ is the set of source
vertices of $v_{i}$.]

33:
$\overline{\pi}_{i}=\overline{\pi}_{i}+\frac{c\overline{\pi}_{u}}{deg(u)}$;

34: end for

35: $h_{i}=0$;

36: end for

37: $STS=false$;

38: end if

39: end if

40: end while

41:end function

42:

43:function Management

44: while There exists non-dangling vertex satisfying $h_{i}>\xi$ do

45: end while

46: $CTRL=false$;

47: while There exists dangling vertex satisfying $h_{i}>\xi$ do

48: end while

49: Terminate all the $K$ Calculations.

50:end function

Algorithm 3 IFP2

Different from IFP1, IFP2 needs to address all of the vertices, both non-
dangling vertices and dangling vertices are assigned to the $K$ threads. IFP2
contains two phases, where phase1 addresses the non-dangling vertices and
phase2 addresses the dangling vertices. Phase1 is similar to IFP1 except that
no mass is pushed to the dangling vertices. While phase1 finished, all of the
non-dangling vertices get converged, phase2 starts and pushes mass to the
dangling vertices. Compared with IFP1, phase2 executes only once and thus some
pushing operation corresponding to dangling vertices are saved. Specifically,
during the $t_{th}$ iteration, there are $\sum\limits_{v\in
V_{D}-V_{U}(t)}|S(v)|$ mass pushing operations saved, thus the total saved
operations is $\sum\limits_{t=0}^{T}\sum\limits_{v\in V_{D}-V_{U}(t)}|S(v)|$.
Generally, the higher proportion of edges corresponding to dangling vertices
the graph contains, the less computation IFP2 generates.

It should be noted that, IFP2 needs some preprocessing, before phase1 running,
edges corresponding to dangling vertices should be invalidated. These work may
increase the time consumption, however, based on appropriate data structure,
the preprocessing cost is at most $o(\sum\limits_{v\in V_{D}}|S(v)|)$.
Experiments in the following will show that, the invaliding operations almost
add no extra CPU time consumption.

##  4 Experiment

In this section, both IFP1 and IFP2 are demonstrated experimentally. Firstly,
we introduce the computing environment and indexes to be used. Then, the
convergence of IFP1 and IFP2 are illustrated, and the comparison with Power
method as well. At last, the results of IFP2 with different threads is
elaborated.

###  4.1 Experiment Setting

All algorithms of this experiment are implemented with C++ on serves with
Intel(R) Xeon(R) Silver 4210R CPU 2.40GHz 40 processors and 190GB memory. The
operation system is Ubuntu 18.04.5 LTS. Six data sets are illustrated in table
1, where $n$, $m$, $n_{d}$, $m_{d}$ and $deg=\frac{m}{n}$ represent the number
of vertices, the number of edges, the number of dangling vertices, the number
of edges corresponding to the dangling vertices and the average degree
respectively. The CPU time consumption $T$ and max relative error
$ERR=\max\limits_{v_{i}\in V}\frac{|\overline{\pi}_{i}-\pi_{i}|}{\pi_{i}}$ are
taken to estimate the algorithms, where The true PageRank value $\pi_{i}$ is
obtained by Power method at the $210_{th}$ iteration. The damping factor
$c=0.85$.

Data Sets | $n$ | $m$ | $n_{d}$ | $m_{d}$ | $deg$
---|---|---|---|---|---
web-Stanford | 281903 | 2312497 | 172 | 410 | 8.21
Stanford-Berkeley | 683446 | 7583376 | 68062 | 994368 | 11.1
web-Google | 916428 | 5105039 | 176974 | 325725 | 5.57
in-2004 | 1382908 | 16917053 | 282306 | 1006484 | 12.23
soc-LiveJournal1 | 4847571 | 68993773 | 539119 | 2092984 | 14.23
uk-2002 | 18520343 | 298113762 | 2760973 | 16029357 | 16.10
Table 1: Data Sets

###  4.2 Convergence

The relation of $ERR$, $T$ and $\xi$ on six data sets is illustrated in Figure
1.

  1. (1)

The blue lines show that $ERR$ has a positive linear relation with $\xi$,
which is consistent with Formula (6). The red lines show that $T$ has a
negative exponential relation with $\xi$, which is consistent with Formula
(5). That red lines marked with triangles are lower than red lines marked with
squares show that IFP2 is faster than IFP1 on all six data sets.

  2. (2)

The blue lines show that $ERR$ scarcely changes when $\xi<10^{-15}$. It seems
both IFP1 and IFP2 have limitation in precision. We believe it is not a
algorithm flaw, but caused by the insufficiency of C++’s DOUBLE data type,
whose significant digit number is 15. $h_{i}<10^{-15}$ can not change
$\overline{\pi}_{i}$, thus $ERR$ stays constantly.

Figure 1: $\xi$ Versus $RES$ and $T$

###  4.3 Comparison with other algorithms

We compare IFP1 and IFP2 with Power method (SPI) [1, 2] and parallel Power
method (MPI)[27]. Both MPI, IFP1 and IFP2 are executed with 38 parallelism.
The relation of $ERR$ and $T$ is illustrated in Figures 2. The CPU time
consumption of preprocessing is illustrated in Table 2. Table 3 illustrates
the CPU time consumption $T$ when $ERR<0.001$.

  1. (1)

The magenta, red and blue lines are lower than green lines in Figure 2, it
shows that both MPI, IFP1 and IFP2 are faster than SPI, parallelizing is an
effective solution to accelerate PageRank computing.

  2. (2)

The magenta lines and red lines are lower than green lines in Figure 2, it
shows that both IFP1 and IFP2 outperform MPI. Since IFP1, IFP2 and MPI are
executed with the same parallelism, we believe the advantages owe to the
higher convergence rate of IFP1 and IFP2.

  3. (3)

For MPI and IFP1, the preprocessing is assigning vertices to threads, for
IFP2, the preprocessing includes assigning and invalidating the edges
corresponding to the dangling vertices. MPI assigns all of the vertices to
threads, thus consumes more time. The time consumption of invalidating is so
few that it can be done in 0.155s even for uk-2002.

  4. (4)

Table 3 shows that IFP2 with 38 parallelism can be at most 50 times as fast as
SPI, and at most 3 times faster than MPI, when $ERR<0.001$.

Figure 2: $T$ Versus $ERR$ Data sets | MPI | IFP1 | IFP2
---|---|---|---
web-Stanford | 0.014 | 0.012 | 0.012
Stanford-Berkeley | 0.041 | 0.030 | 0.031
web-Google | 0.043 | 0.033 | 0.036
in-2004 | 0.070 | 0.046 | 0.054
soc-LiveJournal1 | 0.236 | 0.178 | 0.194
uk-2002 | 0.926 | 0.629 | 0.784
Table 2: The time consumption of preprocessing Data sets | SPI | MPI | IFP1 | IFP2
---|---|---|---|---
web-Stanford | 14.445 | 0.949 | 0.303 | 0.285
Stanford-Berkeley | 27.525 | 1.709 | 1.038 | 0.995
web-Google | 32.182 | 1.936 | 0.718 | 0.671
in-2004 | 64.665 | 4.728 | 1.497 | 1.481
soc-LiveJournal1 | 423.183 | 21.415 | 12.702 | 12.617
uk-2002 | 1108.711 | 85.455 | 24.823 | 22.479
Table 3: The time consumption when $ERR<0.001$

###  4.4 Performance under different parallelism

We execute IFP2 with parallelism 4, 8, 16, 32 and 38. The relation of $ERR$
and $T$ with different parallelism are illustrated in Figure 3. With the
increasing of parallelism, the CPU time consumption $T$ decreases. Generally,
the more the parallelism the faster the computing. IFP2 with 16 parallelism is
faster than MPI with 38 parallelism, which proves IFP2 has higher convergence
rate than Power method again.

Figure 3: $T$ Versus $ERR$ under parallelism $4$,$8$,$16$,$32$ and $38$

##  5 Conclusion

PageRank is a basic problem of graph computation and parallelization is a
feasible solution of accelerating computing. In this paper, we firstly reveal
that PageRank vector is essentially the mass distribution, and base on which
two parallel PageRank algorithms IFP1 and IFP2 are proposed. The most
prominent feature of IFP1 is that it can make full use of the DAG structure,
generally, the larger proportion of dangling vertices, the higher convergence
rate, the larger proportion of unreferenced vertices and weak unreferenced
vertices, the less computation amount. IFP2 pushes mass to the dangling
vertices only once, thus, the more edges corresponding to the dangling
vertices, the more computation decreases. Experiments on six data sets
demonstrate that both IFP1 and IFP2 are faster than the Power method. However,
every coin has two sides, on graph containing no DAG structure such as
undirected graph, neither IFP1 nor IFP2 can outperform the Power method since
the atomic operation is costly. In the future, PageRank on dynamic graph can
be studied.

## References

  * [1] S. Brin, L. Page, Reprint of: The anatomy of a large-scale hypertextual web search engine, Computer networks 56 (18) (2012) 3825–3833.
  * [2] L. Page, S. Brin, R. Motwani, T. Winograd, The pagerank citation ranking: Bringing order to the web., Tech. rep., Stanford InfoLab (1999).
  * [3] D. F. Gleich, Pagerank beyond the web, Siam Review 57 (3) (2015) 321–363.
  * [4] S. Brown, A pagerank model for player performance assessment in basketball, soccer and hockey, arXiv preprint arXiv:1704.00583 (2017).
  * [5] Ying, Hou, Zhou, Dou, Wang, Shao, The changes of central cultural cities based on the analysis of the agglomeration of literati’s footprints in tang and song dynasties, Journal of Geo-information Science 22 (5) (2020) 945–953.
  * [6] S. Kamvar, T. Haveliwala, G. Golub, Adaptive methods for the computation of pagerank, Linear Algebra and its Applications 386 (2004) 51–65.
  * [7] T. Haveliwala, S. Kamvar, D. Klein, C. Manning, G. Golub, Computing pagerank using power extrapolation, Tech. rep., Stanford (2003).
  * [8] S. D. Kamvar, T. H. Haveliwala, C. D. Manning, G. H. Golub, Extrapolation methods for accelerating pagerank computations, in: Proceedings of the 12th international conference on World Wide Web, 2003, pp. 261–270.
  * [9] S. Kamvar, T. Haveliwala, C. Manning, G. Golub, Exploiting the block structure of the web for computing pagerank, Tech. rep., Stanford (2003).
  * [10] G. Wu, Y. Wei, A power-arnoldi algorithm for computing pagerank, Numerical Linear Algebra with Applications 14 (7) (2007) 521–546.
  * [11] K. Avrachenkov, N. Litvak, D. Nemirovsky, N. Osipova, Monte carlo methods in pagerank computation: When one iteration is sufficient, SIAM Journal on Numerical Analysis 45 (2) (2007) 890–904.
  * [12] A. D. Sarma, A. R. Molla, G. Pandurangan, E. Upfal, Fast distributed pagerank computation, in: International Conference on Distributed Computing and Networking, Springer, 2013, pp. 11–26.
  * [13] S. Luo, Distributed pagerank computation: An improved theoretical study, Proceedings of the AAAI Conference on Artificial Intelligence 33 (2019) 4496–4503.
  * [14] S. Luo, Improved communication cost in distributed pagerank computation–a theoretical study, in: International Conference on Machine Learning, PMLR, 2020, pp. 6459–6467.
  * [15] K. Sankaralingam, S. Sethumadhavan, J. Browne, Distributed pagerank for p2p systems, in: High Performance Distributed Computing, 2003. Proceedings. 12th IEEE International Symposium on, 2003, pp. 58–68.  doi:10.1109/HPDC.2003.1210016.
  * [16] Y. Zhu, S. Ye, X. Li, Distributed pagerank computation based on iterative aggregation-disaggregation methods, in: Proceedings of the 14th ACM international conference on Information and knowledge management, 2005, pp. 578–585.
  * [17] S. Stergiou, Scaling pagerank to 100 billion pages, in: Proceedings of The Web Conference 2020, 2020, pp. 2761–2767.
  * [18] H. Ishii, R. Tempo, Distributed randomized algorithms for the pagerank computation, IEEE Transactions on Automatic Control 55 (9) (2010) 1987–2002.
  * [19] H. Ishii, R. Tempo, E.-W. Bai, A web aggregation approach for distributed randomized pagerank algorithms, IEEE Transactions on automatic control 57 (11) (2012) 2703–2717.
  * [20] A. Suzuki, H. Ishii, Efficient pagerank computation via distributed algorithms with web clustering, arXiv preprint arXiv:1907.09979 (2019).
  * [21] T. Charalambous, C. N. Hadjicostis, M. G. Rabbat, M. Johansson, Totally asynchronous distributed estimation of eigenvector centrality in digraphs with application to the pagerank problem, in: 2016 IEEE 55th Conference on Decision and Control (CDC), IEEE, 2016, pp. 25–30.
  * [22] A. Suzuki, H. Ishii, Distributed randomized algorithms for pagerank based on a novel interpretation, in: 2018 Annual American Control Conference (ACC), IEEE, 2018, pp. 472–477.
  * [23] L. Dai, N. M. Freris, Fully distributed pagerank computation with exponential convergence, arXiv preprint arXiv:1705.09927 (2017).
  * [24] I. C. F. Ipsen, T. M. Selee, Pagerank computation, with special attention to dangling nodes, SIAM Journal on Matrix Analysis and Applications 29 (4) (2008) 1281–1296.  arXiv:https://doi.org/10.1137/060664331, doi:10.1137/060664331.
URL https://doi.org/10.1137/060664331

  * [25] Y. Lin, X. Shi, Y. Wei, On computing pagerank via lumping the google matrix, Journal of Computational and Applied Mathematics 224 (2) (2009) 702–708.
  * [26] Z. Zhu, Q. Peng, Z. Li, X. Guan, O. Muhammad, Fast pagerank computation based on network decomposition and dag structure, IEEE Access 6 (2018) 41760–41770.
  * [27] N. T. Duong, Q. A. P. Nguyen, A. T. Nguyen, H.-D. Nguyen, Parallel pagerank computation using gpus, in: Proceedings of the Third Symposium on Information and Communication Technology, 2012, pp. 223–230.
  * [28] Lai, The computations of pagerank scores based on gpgpu environment, Ph.D. thesis, SUN YAT-SEN University (2017).
  * [29] H. Migallón, V. Migallón, J. Penadés, Parallel two-stage algorithms for solving the pagerank problem, Advances in Engineering Software 125 (2018) 188–199.
  * [30] Berkhin, Pavel, A survey on pagerank computing, Internet Mathematics 2 (1) (2005) 73–120.