# Monte Carlo Spin Simulations of Magnetic Noise - The Search for Pivoting

D. L. Mickelsen  Department of Physics and Astronomy, University of
California, Irvine, California, 92697  Ruqian Wu  Department of Physics and
Astronomy, University of California, Irvine, California, 92697  Clare C. Yu
Department of Physics and Astronomy, University of California, Irvine,
California, 92697

###### Abstract

Superconducting quantum interference devices (SQUIDs) show great promise as
quantum bits (qubits) but continue to be hindered by flux noise. The flux
noise power spectra of SQUIDs go as $1/f^{\alpha}$, where $\alpha$ is the
temperature-dependent noise exponent. Experiments find
$0.5\lesssim\alpha\lesssim 1$. Furthermore, experiments find that the noise
power spectra versus frequency at different temperatures pivot about or cross
at a common point for each SQUID. To try to better understand the results and
motivated by experimental evidence that magnetic moments on the surface of
SQUIDS produce flux noise, we present the results of our Monte Carlo
simulations of various spin systems on 2D lattices. We find that only spin
glasses produce $\alpha\sim 1$ at low temperature. We find that aliasing of
the noise power spectra at high frequencies can lead to spectral pivoting if
it is in proximity to a knee at a slightly lower frequency. We show that the
pivot frequency depends on the method of site selection and how often the
magnetization is recorded. The spectral pivoting that occurs in our
simulations is due to aliasing and does not explain the spectral pivoting of
experiments.

##  I Introduction

Superconducting quantum interference devices (SQUIDs) can be used as quantum
bits (qubits). While SQUIDs hold great potential for quantum computing, they
suffer from noise and decoherence. One of the main sources of decoherence is
$1/f$ flux noise [1]. $1/f$ noise is characterized by noise power spectra that
go as $1/f^{\alpha}$, where $f$ is frequency and $\alpha$ is the noise
exponent. The noise exponents for SQUIDs lie in the range
$0.5\lesssim\alpha\lesssim 1$ for $1\text{ K}\gtrsim T\gtrsim 20\text{ mK}$
[2, 3, 4], where the noise exponent increases as temperature decreases.
Fluctuating magnetic moments were proposed as a source of flux noise in SQUIDs
[2]. This is consistent with experimental evidence of surface spins on normal
metals [5] and superconductors [6]. Sendelbach et al. measured a $1/T$
temperature-dependent flux through SQUIDs which is indicative of paramagnetic
spins [6]. Fluctuations of these surface spins cause flux noise because SQUIDs
are highly sensitive magnetometers. Furthermore, cross correlations between
fluctuations in the flux and the inductance in DC SQUIDs indicate that the
spins interact ferromagnetically [7].

Additional evidence for surface spins comes from density functional theory
(DFT) calculations of the oxide layer on the aluminum surface of SQUIDs. Since
SQUIDs are exposed to air, it is reasonable to expect oxygen to be adsorbed on
the surface below approximately $40\text{ K}$. $\text{O}_{2}$ is paramagnetic,
and DFT calculations show that oxygen retains a magnetic moment of $1.8\text{
}\mu_{B}$ after being adsorbed on the surface of sapphire ($\alpha$-Al2O3)
[8]. DFT also finds a low barrier ($\sim 10\text{ mK}$) for spin reorientation
so that spins rotate easily in the easy plane that is perpendicular to the
$\text{O}_{2}$ molecular bond [8, 9]. In addition, these calculations indicate
that the coupling between adsorbed oxygen molecules is ferromagnetic [8].
Monte Carlo simulations of ferromagnetically coupled $\text{O}_{2}$ on the
surface of sapphire are able to produce $1/f$ noise consistent with experiment
at higher temperatures [8]. Evidence of paramagnetic oxygen spins on SQUIDs
come from X-ray absorption spectroscopy (XAS) and x-ray magnetic circular
dichroism (XMCD) experiments that were carried out on thin films of aluminum
and niobium which are typical SQUID materials [9]. The experiments [9] confirm
the DFT predictions [8] that the bond axis of $\text{O}_{2}$ is tilted
$55^{\circ}$ away from the surface normal.

Surface treatments that remove or prevent oxygen adsorption on SQUIDs reduce
flux noise by a factor of four or five [9]. A protective coating of
nonmagnetic ammonia ($\text{NH}_{3}$) prevents the adsorption of
$\text{O}_{2}$ since $\text{NH}_{3}$ has a higher binding energy to
$\text{Al}_{2}\text{O}_{3}$. Ultraviolet irradiation of SQUIDs in an ultrahigh
vacuum removes adsorbed oxygen [9]. Although flux noise is reduced with these
treatments, it is not eliminated.

There is still the question as to how surface spins produce $1/f$ noise. The
Dutta-Horn model of $1/f$ noise [10] is the most common explanation for
$1/f$-type noise. The model assumes that independent, thermally-activated
processes exist, where each process individually produces a Lorentzian power
spectra. A distribution of barrier heights that is slowly varying on the order
of $k_{B}T$ leads to a range of relaxation times and hence, $1/f$ noise. In
the case of spins, spin flips can result from spins hopping in and out of
traps with different preferential spin orientations [11] or hopping between a
local moment state and the conduction band at the interface between a metal
and an insulator [12]. De Sousa proposed that magnetic noise arises from spin
flips of paramagnetic dangling bonds at the amorphous semiconductor/oxide
interface [13]. A range of relaxation times can also come from interactions
between spins [14, 15, 16, 17]. Faoro and Ioffe invoked spin diffusion via
RKKY interactions between spins. Chen and Yu suggested that interacting
surface spins can be modeled as a spin glass [15] while De proposed a model
with small spin clusters in which the spins have ferromagnetic RKKY
interactions [16, 17]. These theoretical models of spins on the surface of
SQUIDs have indeed been able to produce $1/f$ flux noise [11, 13, 14, 12, 15,
16, 17].

While it is well established that surface spins are a source of $1/f$ flux
noise, there are still mysteries associated with this noise. For example, why
is the noise exponent close to 1? While Ising spin glass models yield
$\alpha\lesssim 1$ [15], ferromagnetic spin system simulations find $\alpha>1$
at low temperatures [8].

Another puzzle comes from experiments using Nb SQUIDs. In this case, when the
noise power spectra are plotted as a function of frequency on a log-log plot,
the curves for different temperatures cross or intersect in the vicinity of a
single “pivot” frequency [3]. As a result, as the temperature decreases, the
low frequency noise increases and the high frequency noise decreases. The
noise near the pivot frequency is rather constant as a function of
temperature. For SQUIDs with different geometries, the power spectra for each
device still pivot. The pivot frequency is different for each SQUID, and there
is no clear relation between geometry and pivot frequency. Other experiments
involving SQUIDs of various materials also find spectral pivoting [4].

Futhermore, Anton et al. [3] integrated the (extrapolated) flux noise power
spectra $S_{\Phi}(f)$ over the frequency $f$ from $f_{1}=10^{-4}$ Hz to
$f_{2}=10^{9}$ Hz to obtain the mean square flux noise
$\langle\Phi^{2}\rangle$. They found that the mean square flux noise increased
by two to three orders of magnitude as the temperature increased from 0.1 K to
4 K.

Several models have been proposed to explain spectral pivoting. Spin diffusion
can explain pivoting if the system is close to a phase transition where the
diffusion coefficient is dependent on temperature [18]. The result is a range
of frequencies where power spectra cross as in experiment. Monte Carlo
simulations of Heisenberg spins in a cluster model can also produce spectral
pivoting [19]. This model bases the probability of spin flips on changes in
free energy instead of internal energy as in the standard Monte Carlo method.
This has the effect that lower entropy spin configurations are more favored.
This model produces $1/f^{\alpha}$ power spectra at low frequency and
$1/f^{2}$ power spectra at high frequency. The power spectra do not
intrinsically pivot. To get crossing, the $1/f^{\alpha}$ parts of the power
spectra are extended via extrapolation into the high-frequency range.

In an effort to understand spectral pivoting and what types of interactions
lead to the noise exponent $\alpha\sim 1$, we performed Monte Carlo
simulations of spins on 2D lattices since both DFT simulations [8] and
experiments [9] indicate spins can produce flux noise. Since oxygen spins have
an easy-plane anisotropy perpendicular to the $\text{O}_{2}$ bond [8], we test
various interacting spin models with different anisotropies.

Although we find that spectral pivoting can occur in Monte Carlo simulations
of classical XY and Heisenberg spins, the pivoting in our simulations is an
artifact of the simulations. In plots of noise power versus frequency, a high
frequency pivot occurs because the low-frequency knees are close in frequency
to the high-frequency aliasing of the power spectra, and thus the simulations
do not explain the experimentally observed pivoting. The low-frequency knee
refers to the crossover from a flat noise spectrum at low frequencies to
$1/f^{\alpha}$ at higher frequencies with the noise exponent $\alpha>0$.

Aliasing refers to a well-known phenomenon in signal processing. Suppose a
signal is sampled at equal time intervals, $\Delta$. Then we can define a
Nyquist frequency $f_{c}=1/(2\Delta)$. If we Fourier transform the time
series, the result has components with frequencies outside the frequency range
$-f_{c}<f<f_{c}$ that spuriously contribute to the Fourier transformed signal
in the range between $-f_{c}$ and $f_{c}$ [20]. We say that the true power
spectrum from frequencies outside that range are “folded” into the range. As a
simple example, consider two waves: $g_{1}(t)=\exp(2\pi if_{1}t)$ and
$g_{2}(t)=\exp(2\pi if_{2}t)$ where $f_{2}=f_{1}+(m/\Delta)$ where $m$ is a
positive integer so that $f_{2}>f_{1}$. It is easy to show that $g_{1}=g_{2}$
at time intervals of $\Delta$. So the higher frequency signal $g_{2}(t)$ will
contribute to the Fourier transformed signal at lower frequencies such as
$f_{1}$. In our simulations, aliasing appears as a minimum or flattening of
the noise power spectra at $f_{c}$.

In addition to trying to better understand spectral pivoting, we also use this
paper to explain some of the challenges and pitfalls of Monte Carlo
simulations and what happens when there are deviations from the standard
procedure. The paper is structured in the following manner. In Sec. II, we
describe the Hamiltonians of the spin models, how the simulations are
performed/equilibrated, and how the noise power is analyzed. In Sec. III, we
present our results on noise exponents, noise amplitudes, and why pivoting
occurs in our simulations.

##  II Method

###  II.1 Spin Models

The Ising, XY, and Heisenberg Hamiltonians of the 2D spin systems are given by

$H_{\text{Ising}}=-\sum_{\langle i,j\rangle}J_{ij}s_{i}\cdot s_{j}$ (1)
$H_{\text{XY}}=-\sum_{\langle
i,j\rangle}J_{ij}\mathbf{s}_{i}\cdot\mathbf{s}_{j}$ (2)
$H_{\text{Heis.}}=-\sum_{\langle
i,j\rangle}J_{ij}\mathbf{s}_{i}\cdot\mathbf{s}_{j}-A\sum_{i}(\mathbf{n}_{i}\cdot\mathbf{s}_{i})^{2}$
(3)

where $\mathbf{s}_{i}$ and $\mathbf{s}_{j}$ are classical spins of length 1 on
nearest-neighbor sites $i$ and $j$, respectively. $J_{ij}$ is the spin
exchange coupling. A positive value for $J_{ij}$ indicates a ferromagnetic
interaction. The second term in Eq. (3) is a spin anisotropy term. For each
site $i$, there is a local anisotropy axis $\mathbf{n}_{i}$. To model the
disorder of the SQUID surface, $\mathbf{n}_{i}$ points in a direction that
varies randomly from site to site. The random-axis anisotropic model was
proposed by Harris to describe magnetism in an amorphous material [21].

If the anisotropy $A$ is positive, then $\mathbf{n}_{i}$ is the easy axis for
spins; if the anisotropy $A$ is negative, then $\mathbf{n}_{i}$ is normal to
the easy plane for spin orientation.

Six spin models were simulated: noninteracting spins ($J_{ij}=0$, $A=-1$),
ferromagnet ($J_{ij}=1$), Poisson ferromagnet ($\langle J_{ij}\rangle=1$,
$\sigma^{2}_{J_{ij}}=0.2$), antiferromagnet ($J_{ij}=-1$), spin glass
($\langle J_{ij}\rangle=0,\sigma_{J_{ij}}=1$), and spin glass ferromagnet
($\langle J_{ij}\rangle=0.5,\sigma_{J_{ij}}=1$). While one would think that
the problem of a noninteracting anisotropic Heisenberg spin could be solved
exactly, this is not the case. However, one can make various approximations to
the solution [22]. For the Poisson ferromagnet, the couplings $J_{ij}$ are
chosen in the following way [8]. First random integers $C_{ij}$ are drawn from
a Poisson distribution with a mean of 5. Then the couplings are given by
$J_{ij}$ = 0.2 $C_{ij}$ so that for the Poisson ferromagnet, $\langle
J_{ij}\rangle=1$ and $\sigma^{2}_{J_{ij}}=0.2$. The temperature for this
system is measured in units of $\langle J_{ij}\rangle$.

For the spin glass and spin glass ferromagnet, $J_{ij}$ is chosen from a
Gaussian distribution with a variance $\sigma^{2}_{J_{ij}}$ and the
temperature is measured in units of $\sigma_{J_{ij}}$. One way to obtain the
random interactions associated with a spin glass is via RKKY interactions
[14]. Spins trapped in local moment states at the interface between a metal
and an insulator [12] could interact via RKKY oscillations in the conduction
electron gas of the metal. However, spins associated with adsorbates such as
oxygen molecules on the surface of the native metal oxide are too far away
from the conduction electrons of the metal to interact with each other via
RKKY interactions [8]. The interactions and the corresponding simulated spin
models are summarized in Table 1.

Interaction |  Spin Models
---|---
Noninteracting |  Heisenberg, $A=-1$
$J_{ij}=0$
Spin Glass $\langle J_{ij}\rangle=0$, $\sigma_{J_{ij}}=1$ |  Heisenberg, $A=0$
Heisenberg, $A=-10$
XY
Ising
Antiferromagnet $J_{ij}=-1$ |  Heisenberg, $A=0$
Heisenberg, $A=-10$
XY
Poisson Ferromagnet $\langle J_{ij}\rangle=1$, $\sigma^{2}_{J_{ij}}=0.2$ |  Heisenberg, $A=0$
Heisenberg, $A=-10$
XY
Ferromagnet $J_{ij}=1$ |  Heisenberg, $A=0$
Heisenberg, $A=-10$
XY
Spin Glass Ferromagnet $\langle J_{ij}\rangle=0.5$, $\sigma_{J_{ij}}=1$ |  Heisenberg, $A=0$
Heisenberg, $A=-10$
XY
Table 1: All six types of spin exchange couplings $J_{ij}$ and the
corresponding spins models that were simulated.

###  II.2 Simulation Details

We perform simulations with spins occupying every site of a 16 x 16 square
lattice. The system is initialized with randomly oriented spins. For
anisotropic systems, each run has a unique, but random, set of anisotropic
axes. A spin is allowed to reorient itself according to the Metropolis
algorithm [23]. In this algorithm, a trial move consists of first choosing a
spin on the lattice at random. The initial energy of this site $E_{i}$ is
calculated from the local field produced by its nearest neighbors and the
local anisotropy. A new orientation of the spin is chosen from a random
distribution for the Heisenberg and XY systems. For the Heisenberg systems,
the distribution is random on the unit sphere (the distribution is uniform in
$\phi$ and $\cos\theta$). For XY systems, the distribution is random on the
unit circle (uniform distribution in $\phi$). In the case of Ising systems,
the spin is flipped. The final energy $E_{f}$ of this site with the new spin
orientation is calculated. If the final energy is less than the initial
energy, then the new spin orientation is accepted. However if the final energy
is greater than the initial energy, then the flip is accepted with probability
$\exp[-(E_{f}-E_{i})/(k_{B}T)]$. A random number is generated from a uniform
distribution between 0 and 1; if it is less than the Boltzmann factor, then
the new orientation is accepted. This process continues for the remaining
sites within the lattice. The time it takes for one sweep through the lattice
is one Monte Carlo step (MCS).

The system is equilibrated as described later in section II.5. The total
magnetic moment of this system is obtained by summing over all the spins in
the lattice and is given by $\textbf{m}=\sum_{i=1}^{N}\textbf{s}_{i}$, where
$N$ is the number of sites. After equilibration, the total magnetic moment of
the lattice is recorded at every Monte Carlo step. The system is cooled from
its initial random spin configuration at $T=10$ to $T=0.5$.

###  II.3 Site Selection

For every sweep through the lattice of $N$ spins, $N$ spins are offered the
chance for reorientation, but only $M$ spins are at different sites. There are
two typical methods of selecting sites.

The “every-site method” involves giving each site in the lattice one
opportunity to reorient, i.e., $M=N$. This is the method used in this paper.
In the “random-site method,” $M$ sites are chosen for reorientation at random.
With this method, it is possible to select several sites more than once and
some not at all.

We would expect simulations of magnetic noise to produce white noise power
spectra at high temperatures. For the “every-site method,” this is true. The
power spectrum of Heisenberg spins in the high-temperature limit using the
“random-site method” goes as $1/f$, and the $1/f$ noise power spectra is
entirely due to the site selection method. For a given model, the power
spectra at different temperatures with the “random-site method” pivot at a
lower frequency compared to the “every-site method.”

###  II.4 Time Steps

One time step in a standard Monte Carlo simulation is one Monte Carlo step
(MCS). To simulate a slower recording rate, the total magnetic moment time
series can be recorded after several MCS. A comparison of recording every time
step and every 10 time steps is shown in Fig. 1. Although the maximum
frequency of the power spectra when recording every 10 time steps is smaller
than for recording every time step, the total noise power for both cases is
equal. At a given temperature, the total noise power is equal to the variance
of the total magnetic moment time series which is dimensionless. This means
that the noise power $S_{x}(f)$ has units of MCS. Recording every 10 time
steps causes the aliasing to occur at a lower frequency, but the location of
the low-frequency knee remains unchanged. The effect is that the pivoting
frequency is decreased. Unless otherwise noted, the time steps are recorded at
1 MCS intervals.

(One might be tempted to record the time series more often than one MCS in an
effort to extend to frequencies higher than $0.5\text{ MCS}^{-1}$. However, in
real physical systems, all spins evolve simulataneously so it would be
unphysical to record the total magnetic moment before all spins are given the
opportunity to reorient.)

Figure 1: Noise power $S_{x}(f)$ of the x-component of the total magnetic
moment versus frequency. Noise power spectra averaged over 10 segments for the
2D (a) isotropic ($A=0$) Heisenberg ferromagnet ($J_{ij}=1$) and (b)
noninteracting ($J_{ij}=0$), anisotropic ($A=-1$) Heisenberg models from
recording the total magnetic moment at every time step and every ten time
steps at $T=1$. The units of noise power are MCS, because the total noise
power is dimensionless.

###  II.5 Equilibration

The test for equilibration follows Bhatt and Young’s procedure for the
equilibration of Ising spin glasses [24]. Two independent replicas of each
system with the same exchange couplings and orientation of anisotropic axes
are created and run in parallel. The initial spin configurations for the two
replicas are different and random. For the set of spins {$\textbf{s}_{i}$}
with $N$ lattice sites, the spin autocorrelation function for the replica $n$,
after an equilibration time $t_{0}$, is

$Q^{(n)}(t)=\frac{1}{N}\sum_{i=1}^{N}\textbf{s}_{i}^{(n)}(t_{0})\cdot\textbf{s}_{i}^{(n)}(t_{0}+t),$
(4)

where the summation is over all lattice sites. The spin glass susceptibility
for replica $n$ is calculated as the second moment of this overlap and then
averaged over 200 different realizations of bonds and anisotropy axes. This
disorder average is denoted by $\left[\ldots\right]_{\text{av}}$:

$\chi_{\text{SG}}^{(n)}(t)=\frac{1}{N}\left[\left(\sum_{i=1}^{N}\textbf{s}_{i}^{(n)}(t_{0})\cdot\textbf{s}_{i}^{(n)}(t_{0}+t)\right)^{2}\right]_{\text{av}}.$
(5)

The equilibration time $t_{0}$ is chosen from the sequence
$1,3,10,30,100,300,\ldots{}$, etc. The idea is to compare
$\textbf{s}_{i}^{(n)}(t_{0})$ to $\textbf{s}_{i}^{(n)}(t_{0}+t)$ as
${t\to\infty}$ to see whether $\textbf{s}_{i}^{(n)}(t_{0}+t)$ has lost its
“memory” of $\textbf{s}_{i}^{(n)}(t_{0})$. In practice, the comparison is done
as ${t\to 3t_{0}}$. The spin glass susceptibility in Eq. (5) is averaged over
a length of time $t_{0}$:

$\chi_{\text{SG}}^{(n)}=\frac{1}{Nt_{0}}\left[\sum_{t=t_{0}}^{2t_{0}-1}\left(\sum_{i=1}^{N}\textbf{s}_{i}^{(n)}(t_{0})\cdot\textbf{s}_{i}^{(n)}(t_{0}+t)\right)^{2}\right]_{\text{av}}.$
(6)

The summation over $t$ starts at $t_{0}$ so that the distribution of
$Q^{n}(t)$ is Gaussian. The correlation of the spins at shorter times makes
the distribution deviate from a Gaussian.

For small values of $t_{0}$ and when the system is at low temperatures, there
are few spin fluctuations, so $Q^{(n)}(t)\sim 1$ and
$\chi_{\text{SG}}^{(n)}(t)\sim N$. This is in agreement with simulations.

We can also calculate $\chi^{(n)}_{\text{SG}}$ in the high-temperature limit.
We start with two Heisenberg spins $\textbf{s}_{1}$ and $\textbf{s}_{2}$ of
length 1 that represent $\textbf{s}_{i}^{(n)}(t_{0})$ and
$\textbf{s}_{i}^{(n)}(t_{0}+t)$, respectively, in Eq. (6). Without loss of
generality, we choose one spin to be along the z-axis. The average square of
the dot product is calculated as

$\begin{split}\langle\left(\textbf{s}_{1}\cdot\textbf{s}_{2}\right)^{2}\rangle_{\text{Heis.}}&=\langle
s_{1}^{2}s_{2}^{2}\cos^{2}(\theta)\rangle\\\
&=\langle\cos^{2}(\theta)\rangle,\end{split}$ (7)

where $\theta$ is the angle between $\textbf{s}_{1}$ and $\textbf{s}_{2}$. The
angle $\theta$ is also the polar angle of the spin. In our simulations,
$\cos(\theta)$ is chosen from a uniform distribution, so Eq. (7) can be
simplified:

$\begin{split}\langle\left(\textbf{s}_{1}\cdot\textbf{s}_{2}\right)^{2}\rangle_{\text{Heis.}}&=\langle\cos^{2}(\theta)\rangle\\\
&=\frac{1}{2}\int_{-1}^{1}\cos^{2}(\theta)d(\cos(\theta))\\\
&=\frac{1}{3}\end{split}$ (8)

We repeat the process for an XY spin with an angle $\phi$ between spins that
is chosen from a uniform distribution:

$\begin{split}\langle\left(\textbf{s}_{1}\cdot\textbf{s}_{2}\right)^{2}\rangle_{\text{XY}}&=\langle
s_{1}^{2}s_{2}^{2}\cos^{2}(\phi)\rangle\\\ &=\langle\cos^{2}(\phi)\rangle\\\
&=\frac{1}{2\pi}\int_{0}^{2\pi}\cos^{2}(\phi)d\phi\\\
&=\frac{1}{2}.\end{split}$ (9)

For Ising spins,

$\begin{split}\langle\left(\textbf{s}_{1}\cdot\textbf{s}_{2}\right)^{2}\rangle_{\text{Ising}}&=\langle\left(\pm
s_{1}s_{2}\right)^{2}\rangle\\\ &=\langle s_{1}^{2}s_{2}^{2}\rangle\\\
&=1.\end{split}$ (10)

From Eqs. (8), (9), and (10), we can see that for spins with $m$ components,

$\langle\left(\textbf{s}_{1}\cdot\textbf{s}_{2}\right)^{2}\rangle=\frac{1}{m}.$
(11)

Combining Eq. (11) with Eq. (6), for high temperatures, we get
$\chi^{(n)}_{\text{SG}}=\frac{1}{m}$ which is seen in simulations.

We then define the average of the two single replica susceptibilities

$\overline{\chi}_{\text{SG}}=\frac{1}{2}\left(\chi_{\text{SG}}^{(1)}+\chi_{\text{SG}}^{(2)}\right)$
(12)

as the two times spin glass susceptibility.

The spin glass susceptibility may also be calculated from the spin overlap of
the two different replicas. The mutual overlap between the spins
$\textbf{s}_{i}^{(1)}$ and $\textbf{s}_{i}^{(2)}$ of the two replicas is

$Q(t)=\frac{1}{N}\sum_{i=1}^{N}\textbf{s}_{i}^{(1)}(t_{0}+t)\cdot\textbf{s}_{i}^{(2)}(t_{0}+t).$
(13)

The spin glass susceptibility is calculated from the spin overlap:

$\chi_{\text{SG}}=\frac{1}{Nt_{0}}\left[\sum_{t=t_{0}}^{2t_{0}-1}\left(\sum_{i=1}^{N}\textbf{s}_{i}^{(1)}(t_{0}+t)\cdot\textbf{s}_{i}^{(2)}(t_{0}+t)\right)^{2}\right]_{\text{av}}.$
(14)

For all temperatures, as the equilibration time is approached, the spin glass
susceptibilities converge; the two times susceptibility
$\overline{\chi}_{\text{SG}}$ (Eq. 12) approaches the true susceptibility from
above and the replica susceptibility $\chi_{\text{SG}}$ (Eq. 14) from below.

After sufficiently long equilibration times, $\chi_{\text{SG}}$ and
$\overline{\chi}_{\text{SG}}$ agree. We define the system to be equilibrated
if

$\Delta\chi_{SG}=\frac{\displaystyle|\chi_{\text{SG}}-\overline{\chi}_{\text{SG}}|}{\displaystyle\frac{1}{2}\left(\chi_{\text{SG}}+\overline{\chi}_{\text{SG}}\right)}$
(15)

is less than 5% for three consecutive times in the $t_{0}$ sequence; then we
declare it equilibrated at the fourth time. For example, if the last three
equilibration times are $t_{1}=3000$, $t_{2}=10000$, $t_{3}=30000$, then the
equilibration time $t_{4}=100000$. At each temperature, the initial
equilibration time is $1\times 10^{5}$ MCS, and it is increased if the system
is not equilibrated.

###  II.6 Noise Power

The time series of each component of the total magnetic moment is given by
$m_{a}(t_{j})$, where $a=x,y,z$. The deviation from the average is $\delta
m_{a}(t_{j})=m_{a}(t_{j})-\overline{m_{a}(t_{j})}$. The noise power spectral
density can be determined from the Fourier transform of the autocorrelation
function of the time series. The autocorrelation is given by

$C_{a}(t_{k})=\frac{1}{N_{\tau}}\sum_{j=0}^{N_{\tau}-1}\delta
m_{a}(t_{j})\delta m_{a}^{*}(t_{j}-t_{k}),$ (16)

where $\delta m_{a}^{*}(t_{j}-t_{k})$ is the complex conjugate of $\delta
m_{a}(t_{j}-t_{k})$. For a discrete time series $m_{a}(t_{j})$ of length
$N_{\tau}$, the discrete inverse Fourier transform is given by

$\delta
m_{a}(t_{j})=\frac{1}{N_{\tau}}\sum_{k=0}^{N_{\tau}-1}\tilde{m}_{a}(f_{k})e^{2\pi
if_{k}t_{j}/N_{\tau}},$ (17)

where $\tilde{m}_{a}(f_{k})$ is the Fourier transform of the time series.
Using Eq. (17), the autocorrelation function becomes

$C_{a}(t_{k})=\frac{1}{N_{\tau}^{3}}\sum_{j=0}^{N_{\tau}-1}\left[\left(\sum_{l=0}^{N_{\tau}-1}\tilde{m}_{a}(f_{l})e^{2\pi
if_{l}t_{j}/N_{\tau}}\right)\right.\times\\\
\left.\left(\sum_{n=0}^{N_{\tau}-1}\tilde{m}^{*}_{a}(f_{n})e^{-2\pi
if_{n}(t_{j}-t_{k})/N_{\tau}}\right)\right].$ (18)

This can be simplified to

$C_{a}(t_{k})=\frac{1}{N_{\tau}^{2}}\sum_{n=0}^{N_{\tau}-1}|\tilde{m}_{a}(f_{n})|^{2}e^{2\pi
if_{n}t_{k}/N_{\tau}}.$ (19)

Taking the Fourier transform of Eq. (19) gives the periodogram estimate [20]
for computing the noise power $S_{a}(f_{k})$ for the axes $a=x,y,z$:

$S_{a}(f_{k})=\frac{1}{N_{\tau}}\frac{1}{N}|\tilde{m}_{a}(f_{k})|^{2},$ (20)

where the power spectra is divided by $N$ so that $S_{a}(f_{k})$ is the noise
power per site for spin component $a$. The power spectra are normalized so
that at a given temperature, the total noise power is equal to the variance of
the total magnetic moment time series divided by the number of sites.

The Fourier transform $\tilde{m}_{a}(f_{k})$ is computed using the C
subroutine library FFTW [25]. At a given temperature, the time series $\delta
m_{a}(t_{j})$ for a given run is split into either 10 or 100 segments (blocks)
of equal length. The power spectrum is found for each segment and is averaged
over these segments to give a smoother power spectrum. At each temperature,
the spectra are averaged over 200 independent runs. The power spectra used in
all plots shown in this paper are averaged over 100 segments.

##  III Results

###  III.1 Low-Frequency Knee

At lower frequencies of the power spectra, there are “knees” where the power
spectra transition from $A^{2}/f^{\alpha}$ to white noise at low frequency as
shown in Fig. 2. The frequency of the knee decreases as the temperature
decreases.

Figure 2: Noise power $S_{x}(f)$ of the x-component of the total magnetic
moment of a 2D isotropic ($A=0$) Heisenberg ferromagnet ($J_{ij}$=1) versus
frequency. The low-frequency knee of the power spectrum averaged over 100
segments is shown for $T=2$.

The low-frequency knee is due to finite size effects of the lattice [26]. The
correlation length increases as the temperature approaches the transition
temperature $T_{c}$ associated with the relevant order parameter. The
relaxation time increases as the correlation length increases, so the system
takes longer to equilibrate. In previous simulations of the ferromagnetic
Ising model and 5-state Potts model, it was found that the knee frequency is
proportional to the system’s minimum relaxation rate $\tau^{-1}$, i.e., the
system’s maximum relaxation time [26]. This means that as the temperature
decreases and approaches $T_{c}$, the knee frequency decreases. Since the
relaxation time increases with system size, we expect the knee frequency to
decrease with system size [26]. (So in real materials or devices where the
lattice size is huge, one does not see a knee in the noise power spectrum.)
One should look at the noise in the relevant order parameter to see the knee
frequency decrease. The order parameters are magnetization for the
ferromagnetic systems, staggered magnetization for the antiferromagnetic
systems, and the spin glass order parameter $q$ for spin glass systems. We
investigate this below.

According to the Mermin-Wagner-Hohenberg theorem [27, 28], two-dimensional
Heisenberg spin systems should not undergo a phase transition. The order
parameters of the ferromagnet, antiferromagnet, and spin glass systems show
that Heisenberg systems exhibit ordering at a positive $T_{c}$. This
discrepancy is because the Mermin-Wagner-Hohenberg theorem holds in the
thermodynamic limit, while these simulations are for finite lattices.
Simulations show that as the system size increases, the value of $T_{c}$
decreases toward zero. This system size dependence is shown in Fig. 3. We
define the magnetic susceptibility of the total magnetic moment per site as

$\chi_{\rm
total}=\frac{1}{k_{B}T}\left(\left\langle\left|\frac{m}{N}\right|^{2}\right\rangle-\left\langle\left|\frac{m}{N}\right|\right\rangle^{2}\right),$
(21)

where we set $k_{B}=1$, $T$ is temperature, $N$ is the number of lattice
sites, and the magnitude of the total magnetization
$\left|m\right|=\sqrt{m_{x}^{2}+m_{y}^{2}+m_{z}^{2}}$. As the system size
increases, the peaks in the magnetic susceptibility occur at lower
temperatures.

Figure 3: Magnetic susceptibility of the magnitude of the total magnetic
moment per site versus temperature for the 2D Heisenberg Ferromagnet as a
function of system size.

###  III.2 Aliasing

At higher frequencies near $0.5\text{ MCS}^{-1}$, there is an upturn of the
noise power due to aliasing that is shown in Fig. 4. The aliasing is due to
the periodicity of the factor $e^{-2\pi ift}$ used in the Fourier transform
[20]. When calculating the discrete Fourier transform, frequency components of
the power spectra that are greater than $0.5\text{ MCS}^{-1}$ are translated
into the range $0\text{ MCS}^{-1}<f<0.5\text{ MCS}^{-1}$ [20]. As we explain
later in this section, under certain conditions, this aliasing can cause
spectral pivoting where the power spectra at different temperatures cross
within a narrow range of frequency.

Figure 4: Noise power $S_{x}(f)$ of the x-component of the total magnetic
moment of a 2D isotropic ($A=0$) Heisenberg ferromagnet ($J_{ij}$=1) versus
frequency. Labeled regions of the power spectrum averaged over 100 segments at
$T=2$.

###  III.3 Noise Exponents

To determine the noise amplitude ($A^{2}$) and the noise exponent ($\alpha$),
the function $A^{2}/f^{\alpha}$ is fit to the region of the power spectra that
is linear on a log-log plot and that lies between the low-frequency knee and
high-frequency upturn due to aliasing as shown in Fig. 4. More details about
this fitting process can be found in Appendix A. In performing fits,
$A^{2}/f^{\alpha}$ is fit to the 10-averaged power spectra, so the fits are
linear on log-log plots. Each x, y, and z component of the spin results in a
power spectrum, and the amplitudes and exponents are determined for each
component. The noise amplitudes and exponents as a function of temperature are
shown in Figs. 5 and 6, respectively. In Fig. 6, the shaded region indicates
the experimental range of noise exponents where $0.5\leq\alpha\leq 1$ and
$1\leq T\leq 2$. (The lower bound of $T\geq 1$ is set by the energy scale of
the exchange constant $J_{ij}$.) From the plot of noise exponents, we can see
that the spin glass systems are the most consistent with experiment.

Figure 5: Noise amplitude averaged over spin components as a function of
temperature for ($0.5\leq T\leq 10$). Figure 6: Noise exponents averaged over
spin components as a function of temperature for ($0.5\leq T\leq 10$). The
shaded region indicates the experimental range of $0.5\leq\alpha\leq 1$ and
$1\leq T\leq 2.$

###  III.4 Pivoting

If the knee and aliasing upturn are close in frequency, then the two regions
around these features overlap, and the frequency range of the power law fit is
reduced. This effect is shown in Fig. 7 for the isotropic Heisenberg
ferromagnet at $T=10$.

Figure 7: Noise power $S_{x}(f)$ of the x-component of the total magnetic
moment of a 2D isotropic ($A=0$) Heisenberg ferromagnet ($J_{ij}=1$) versus
frequency. Black lines indicate the frequency fit range of the noise power
spectra averaged over 100 segments at $T=10$ and $T=1$.

Another effect of the aliasing upturn being close to the knee is that the
power spectra at different temperature pivot about a common frequency as seen
in Figs. 8 and 9. In Fig. 8, the noise power as a function of frequency is
shown for different temperatures for the isotropic ferromagnetic, isotropic
antiferromagnetic, isotropic spin glass, and noninteracting anisotropic
$(A=-1)$ Heisenberg spin models. Since the power spectra for the
antiferromagnetic spin system do not cross at a common frequency, they do not
pivot. The power spectra for the isotropic ferromagnetic, isotropic spin
glass, and noninteracting, anisotropic $(A=-1)$ systems pivot at high
temperature. To examine this further, in Fig. 9, the noise power as a function
of frequency is shown for the Heisenberg ferromagnet for $0.5\leq T\leq 10$.
Over this large range of temperature, the power spectra does not pivot. The
inset of Fig. 9 shows pivoting for the high-temperature range $4\leq T\leq
10$. Fig. 10 shows the pivoting of the power spectra for the 2D Ising spin
glass for $1.7\leq T\leq 2.3$. The anomalous high-temperature pivoting of
power spectra for $4\leq T\leq 10$ is presented in Appendix B for the 2D Ising
spin glass.

Figure 8: Noise power $S_{x}(f)$ of the x-component of the total magnetic
moment versus frequency. Spectral pivoting of power spectra averaged over 100
segments of the 2D (a) noninteracting ($J_{ij}=0$) anisotropic ($A=-1$), (b)
isotropic ($A=0$) spin glass ($\langle J_{ij}\rangle=0,\sigma_{J_{ij}}=1$),
(c) isotropic ($A=0$) antiferromagnetic ($J_{ij}=-1$), (d) isotropic ($A=0$)
Poisson ferromagnetic ($\langle J_{ij}\rangle=1$, $\sigma^{2}_{J_{ij}}=0.2$),
(e) isotropic ($A=0$) ferromagnetic ($J_{ij}=1$), and (f) isotropic ($A=0$)
spin glass ferromagnetic ($\langle J_{ij}\rangle=0.5,\sigma_{J_{ij}}=1$)
Heisenberg models for $1\leq T\leq 10$.

Figure 9: Noise power $S_{x}(f)$ of the x-component of the total magnetic
moment of a 2D isotropic $(A=0)$ Heisenberg ferromagnet ($J_{ij}=0$) versus
frequency. Spectral pivoting of the power spectra averaged over 100 segments
for $0.5\leq T\leq 10$. The inset shows the power spectra averaged over 100
segments for $4\leq T\leq 10$. Figure 10: Noise power $S(f)$ of the total
magnetic moment of a 2D Ising spin glass ($\langle
J_{ij}\rangle=0,\sigma_{J_{ij}}=1$) versus frequency. Spectral pivoting of
power spectra averaged over 100 segments for $1.7\leq T\leq 2.3$.

Spectral pivoting occurs for all systems tested except for the Heisenberg
antiferromagnets with $A=0$ and $A=-1$ as well as the XY antiferromagnet (see
Fig. 8). Spectral pivoting is most evident at high temperatures well above the
magnetic transition temperature $T_{c}$. This is consistent with experiments
on flux noise where there is no conclusive evidence for a magnetic phase
transition, indicating that the experimentally observed pivoting occurs at
temperatures above any magnetic phase transition.

For our simulated spin systems that exhibit pivoting, as their temperature
decreases, the spins change their orientations more slowly; the low-frequency
noise power increases, and the noise power at high frequencies decreases. Note
that the total noise power i.e. the integrated area under spectral density
curve, is the same for all temperatures for the noninteracting Heisenberg
model with anisotropy $A=-1$ and spin glass models with $T_{c}=0$ is zero. As
a result, increasing noise power at low frequencies means decreasing noise
power at high frequencies.

We find that that the crossing frequency of the power spectra has a weak
temperature dependence of the form $f_{c}=B\cdot T+f_{0}$. For noninteracting
Heisenberg spins with anisotropy $A=-1$, $B=0.023$ and $f_{0}=-0.16$ for
$0.5\leq T\leq 2.5$. For comparison, $B=0.11$ and $f_{0}=-0.12$ for the
isotropic Heisenberg ferromagnet for $1.25\leq T\leq 2.5$. Spectral pivoting
is an artifact that occurs at high temperature where the knee and aliasing
regions are close to each other as seen in the inset of Fig. 9. At low
temperatures where the knee and aliasing are not close, we do not see pivoting
which can be seen in Fig. 9.

By changing a few simulation parameters, the pivoting can be affected. The
pivoting of the power spectra of the Heisenberg ferromagnet when recording the
magnetic moment at every time step using the “every-site method” is shown in
Fig. 11(a). Using the method of randomly selecting sites for reorientation
outlined in Sec. II.3 results in a lower crossing frequency by lowering the
frequency where aliasing occurs as shown in Fig. 11(b). Recording the magnetic
moment time series every ten time steps as described in Sec. II.4 also lowers
the crossing frequency compared to recording the magnetic moment time series
at every time step as shown in Fig. 11(c).

Figure 11: Noise power $S_{x}(f)$ of the x-component of the total magnetic
moment of a 2D isotropic ($A=0$) Heisenberg ferromagnet ($J_{ij}=1$) versus
frequency. Spectral pivoting of x-component of the power spectra averaged over
100 segments resulting from recording the magnetic moment time series (a) at
every time step using the “every-site method,” (b) at every time step using
the “random-site method,” and (c) at every tenth time step using the “every-
site method.”

###  III.5 Mean-Square Flux Noise

As we mentioned earlier, the mean-square flux noise is given by

$\langle\Phi^{2}\rangle=\int_{f_{1}}^{f_{2}}S_{\Phi}(f)df.$ (22)

In experiments by Anton et al., the mean-square flux noise in SQUIDs was found
to increase with increasing temperature with $f_{1}=10^{-4}\text{ Hz}$ and
$f_{2}=10^{9}\text{ Hz}$ [3]. In our simulations, the mean-square flux noise
is equivalent to the total noise power with $f_{1}=0\text{ MCS}^{-1}$ and
$f_{2}=0.5\text{ MCS}^{-1}$. The total noise power is calculated for the
isotropic ferromagnetic, noninteracting anisotropic, and isotropic
antiferromagnetic Heisenberg systems for $0.5\leq T\leq 10$. The total noise
power as a function of temperature is shown in Fig. 12.

In our simulations, only the Heisenberg antiferromagnet shows a total noise
power that increases with temperature, though the curvature of the plot
differs from that of Anton et al. [3]. Furthermore, Anton et al. [3] find that
the mean square flux noise increases by two to three orders of magnitude as
the temperature increases from 0.1 K to 4 K while our antiferromagnetic
simulations find an increase of less than one order of magnitude.

Figure 12: Total noise power ($\int_{0\text{ MCS}^{-1}}^{0.5\text{
MCS}^{-1}}S_{x}(f)df$) of the total magnetic moment versus temperature of the
2D isotropic ($A=0$) ferromagnetic ($J_{ij}=1$), noninteracting ($J_{ij}=0$)
anisotropic ($A=-1$), and isotropic ($A=0$) antiferromagnetic ($J_{ij}=-1$)
Heisenberg models for $0.5\leq T\leq 10$.

##  IV Summary and Discussion

Monte Carlo simulations of various spin models with nearest-neighbor exchange
were carried out on 2D lattices to determine which interactions yield
$\alpha\sim 1$ at low temperatures as well as to find cases of spectral
pivoting. We find that the spin glass systems produce noise exponents that
best match experiment [2, 3, 4]. In simulations, pivoting at high frequencies
occurs as a result of the proximity of a low frequency knee and the aliasing
of the noise power spectra. In experiments, aliasing can be avoided by using
low-pass filters [29]. This paper does not explain pivoting seen in
experiments but it can explain the pivoting we see in some simulations. We
note that in our simulations pivoting is most evident at temperatures high
compared with the magnetic transition temperature. Presumably this is
consistent with experiments where there is no consistent evidence of a
magnetic transition.

This paper also does not explain why the mean-square flux noise increases with
temperature in the experiment by Anton et al. [3]. We find that this increase
is characteristic of antiferromagnetic interactions between spins, for which
there is no other experimental evidence. Furthermore, the increase in the mean
square flux noise is 2 to 3 orders of magnitude in experiment compared to less
than one order of magnitude in simulations. One way to interpret the
experimental results would be to say that it implies an increase by a couple
of orders of magnitude in the number of spins as the temperature increases;
this is an interpretation that would be difficult to explain in the context of
most spin models that have a fixed number of spins. Thus, the experimental
findings imply that additional sources of flux noise become more prominent
with increasing temperature. What these additional sources are is unknown.

##  V Acknowledgements

We thank Hervé Carruzzo for helpful discussions. RW acknowledges the support
by U.S. Department of Energy, Office of Science, Basic Energy Sciences, under
Award DE-FG02-05ER46237. This work was performed in part at the Aspen Center
for Physics, which is supported by National Science Foundation grant
PHY-1607611.

##  Appendix A Fitting

The process for fitting $A^{2}/f^{\alpha}$ to the noise power spectra to
determine the noise exponent starts with dividing the noise power spectra into
frequency segments as shown in Fig. 13(a). The segments from $f_{i}$ to
$f_{i+1}$ range from $i=1$ to $i=N_{\text{segments}}-1$. The segments follow
the condition $f_{i+1}/f_{i}=10^{0.1}$. The function $A^{2}/f^{\alpha}$ is fit
to each segment starting at the lowest frequency segment.

Figure 13: Fitting method using the power spectrum averaged over 100 segments
of the isotropic ($A=0$) Heisenberg ferromagnet ($J_{ij}$=1) at $T=0.9$. (a)
Segmented power spectra: noise power $S_{x}(f)$ of the x-component of the
total magnetic moment versus frequency. (b) Relative change in exponent versus
frequency for the segmented power spectra. (c) Percent error between the power
spectra fit and data versus frequency. (d) Noise power $S_{x}(f)$ of the
x-component of the total magnetic moment versus frequency. The data segment
used for fitting is shown in black.

The relative change in the noise exponent between the current segment and the
previous one is calculated. An example of the relative change is shown in Fig.
13(b). If the relative change is less than $1\%$ and the noise exponent is
greater than 0.2, then the frequency range of the segment is noted. The lower
and upper limits of the new frequency fit range are found by combining all of
these segments. In the figure, this is shown by two vertical dotted lines. A
fit is performed in this region.

The percent error of the fit is the percent difference between
$A^{2}/f^{\alpha}$ evaluated at a particular $f$ and the noise power at $f$
from the data. The percent error is evaluated for all frequency data points
within the fitted frequency range and is shown in Fig. 13(c). A new frequency
range is defined by the maximum and minimum frequencies corresponding to
percent differences of less than $0.03\%$ as shown by vertical dotted lines.

This frequency range is shortened by increasing the lower limit by $30\%$ and
decreasing the upper limit by $20\%$. By reducing the frequency range, the fit
region will not be within the knee and aliasing regions. These percentages
that were found by trial and error work well for all models presented. The
final fit of $A^{2}/f^{\alpha}$ is performed within this new frequency range.
The power spectra and final fit region of the power spectra are shown in Fig.
13(d).

##  Appendix B High-Temperature Pivot of the Ising Spin Glass

As seen in Figs. 8 and 9 in Sec. III.4, the noise power spectra tend to pivot
at high temperature. The noise power as a function of frequency at high
temperature for the 2D Ising spin glass is shown in Fig. 14. Although the
power spectra pivots, the noise exponent is negative. This is because Ising
spin flips are always $180^{\circ}$ rotations. At high temperature, Ising
spins have a high probability of flipping at every time step which increases
the noise at $f=0.5\text{ MCS}^{-1}$. In the infinite-temperature limit, the
power spectrum would be a delta function peaked at $f=0.5\text{ MCS}^{-1}$.

Figure 14: Noise power $S(f)$ of the total magnetic moment of a 2D Ising spin
glass ($\langle J_{ij}\rangle=0,\sigma_{J_{ij}}=1$) versus frequency. Spectral
pivoting of power spectra averaged over 100 segments for $4\leq T\leq 10$.

##  Appendix C Crossing Condition

The experimental results of Anton et al. indicate that at high temperatures,
the temperature dependence of the noise amplitude can be described by
$A^{2}(T)=A^{2}_{0}T^{\gamma}$ and the noise exponent can be described by
$\alpha(T)=\alpha_{0}\ln(T)+\alpha_{1}$ [3]. Monte Carlo simulations also find
that these relations hold at high temperature. In addition, in both theory and
experiment, the crossing frequency $f_{c}$ as a function of temperature does
not change significantly. This can be expressed as

$\frac{dS(f,T)}{dT}\bigg{|}_{f=f_{c}}=0.$ (23)

Using $S(f,T)=A^{2}(T)/f^{\alpha(T)}$, we can relate the noise amplitudes to
noise exponents:

$\displaystyle\frac{dS(f,T)}{dT}\bigg{|}_{f=f_{c}}$
$\displaystyle=\frac{d}{dT}\left(\frac{A^{2}(T)}{f^{\alpha(T)}}\right)\bigg{|}_{f=f_{c}}$
(24)
$\displaystyle=\left(A^{2}(T)\frac{d}{dT}\frac{1}{f^{\alpha(T)}}+\frac{1}{f^{\alpha(T)}}\frac{d}{dT}A^{2}(T)\right)\bigg{|}_{f=f_{c}}$
$\displaystyle=\left(A^{2}(T)\left(-\frac{\ln(f)}{f^{\alpha(T)}}\frac{d\alpha(T)}{dT}\right)+\frac{A^{2}_{0}}{f^{\alpha(T)}}T^{\gamma-1}\gamma\right)\bigg{|}_{f=f_{c}}$
$\displaystyle=-S(f_{c})\ln(f_{c})\frac{\alpha_{0}}{T}+\frac{S(f_{c})}{T}\gamma$
$\displaystyle=0.$

This gives the crossing condition

$\gamma=\alpha_{0}\ln(f_{c}).$ (25)

## References

  * Bialczak  _et al._ [2007] R. C. Bialczak, R. McDermott, M. Ansmann, M. Hofheinz, N. Katz, E. Lucero, M. Neeley, A. D. O’Connell, H. Wang, A. N. Cleland, and J. M. Martinis, $1/f$ flux noise in josephson phase qubits, Phys. Rev. Lett. 99, 187006 (2007).
  * Wellstood  _et al._ [1987] F. C. Wellstood, C. Urbina, and J. Clarke, Low‐frequency noise in dc superconducting quantum interference devices below 1 k, Applied Physics Letters 50, 772 (1987).
  * Anton  _et al._ [2013] S. M. Anton, J. S. Birenbaum, S. R. O’Kelley, V. Bolkhovsky, D. A. Braje, G. Fitch, M. Neeley, G. C. Hilton, H.-M. Cho, K. D. Irwin, F. C. Wellstood, W. D. Oliver, A. Shnirman, and J. Clarke, Magnetic flux noise in dc squids: Temperature and geometry dependence, Phys. Rev. Lett. 110, 147002 (2013).
  * Kempf  _et al._ [2016] S. Kempf, A. Ferring, and C. Enss, Towards noise engineering: Recent insights in low-frequency excess flux noise of superconducting quantum devices, Appl. Phys. Lett. 109, 162601 (2016).
  * Bluhm  _et al._ [2009] H. Bluhm, J. A. Bert, N. C. Koshnick, M. E. Huber, and K. A. Moler, Spinlike susceptibility of metallic and insulating thin films at low temperature, Phys. Rev. Lett. 103, 026805 (2009).
  * Sendelbach  _et al._ [2008] S. Sendelbach, D. Hover, A. Kittel, M. Mück, J. M. Martinis, and R. McDermott, Magnetism in squids at millikelvin temperatures, Phys. Rev. Lett. 100, 227006 (2008).
  * Sendelbach  _et al._ [2009] S. Sendelbach, D. Hover, M. Mück, and R. McDermott, Complex inductance, excess noise, and surface magnetism in dc squids, Phys. Rev. Lett. 103, 117001 (2009).
  * Wang  _et al._ [2015] H. Wang, C. Shi, J. Hu, S. Han, C. C. Yu, and R. Q. Wu, Candidate source of flux noise in squids: Adsorbed oxygen molecules, Phys. Rev. Lett. 115, 077002 (2015).
  * Kumar  _et al._ [2016] P. Kumar, S. Sendelbach, M. A. Beck, J. W. Freeland, Z. Wang, H. Wang, C. C. Yu, R. Q. Wu, D. P. Pappas, and R. McDermott, Origin and reduction of $1/f$ magnetic flux noise in superconducting devices, Phys. Rev. Applied 6, 041001(R) (2016).
  * Dutta and Horn [1981] P. Dutta and P. M. Horn, Low-frequency fluctuations in solids: $\frac{1}{f}$ noise, Rev. Mod. Phys. 53, 497 (1981).
  * Koch  _et al._ [2007] R. H. Koch, D. P. DiVincenzo, and J. Clarke, Model for $1/f$ flux noise in squids and qubits, Phys. Rev. Lett. 98, 267003 (2007).
  * Choi  _et al._ [2009] S. K. Choi, D.-H. Lee, S. G. Louie, and J. Clarke, Localization of metal-induced gap states at the metal-insulator interface: Origin of flux noise in squids and superconducting qubits, Phys. Rev. Lett. 103, 197001 (2009).
  * de Sousa [2007] R. de Sousa, Dangling-bond spin relaxation and magnetic $1/f$ noise from the amorphous-semiconductor/oxide interface: Theory, Phys. Rev. B 76, 245306 (2007).
  * Faoro and Ioffe [2008] L. Faoro and L. B. Ioffe, Microscopic origin of low-frequency flux noise in josephson circuits, Phys. Rev. Lett. 100, 227005 (2008).
  * Chen and Yu [2010] Z. Chen and C. C. Yu, Comparison of ising spin glass noise to flux and inductance noise in squids, Phys. Rev. Lett. 104, 247204 (2010).
  * De [2014] A. De, Ising-glauber spin cluster model for temperature-dependent magnetization noise in squids, Phys. Rev. Lett. 113, 217002 (2014).
  * De [2019] A. De, $1/f$ flux noise in low-${T}_{c}$ squids due to superparamagnetic phase transitions in defect clusters, Phys. Rev. B 99, 024305 (2019).
  * Lanting  _et al._ [2014] T. Lanting, M. H. Amin, A. J. Berkley, C. Rich, S.-F. Chen, S. LaForest, and R. de Sousa, Evidence for temperature-dependent spin diffusion as a mechanism of intrinsic flux noise in squids, Phys. Rev. B 89, 014503 (2014).
  * Davis and Chamberlin [2018] B. F. Davis and R. V. Chamberlin, 1/f noise from a finite entropy bath: comparison with flux noise in squids, Journal of Statistical Mechanics: Theory and Experiment 2018, 103206 (2018).
  * Press  _et al._ [2007] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,  _Numerical Recipes in C: The Art of Scientific Computing_ , 3rd ed. (Cambridge University Press, New York, 2007).
  * Harris  _et al._ [1973] R. Harris, M. Plischke, and M. J. Zuckermann, New model for amorphous magnetism, Phys. Rev. Lett. 31, 160 (1973).
  * Brown [1963] W. F. Brown, Thermal fluctuations of a single-domain particle, Phys. Rev. 130, 1677 (1963).
  * Metropolis  _et al._ [1953] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, Equation of state calculations by fast computing machines, J. Chem. Phys. 21, 1087 (1953).
  * Bhatt and Young [1988] R. N. Bhatt and A. P. Young, Numerical studies of ising spin glasses in two, three, and four dimensions, Phys. Rev. B 37, 5606 (1988).
  * Frigo and Johnson [2005] M. Frigo and S. G. Johnson, The design and implementation of FFTW3, Proceedings of the IEEE 93, 216 (2005), special issue on “Program Generation, Optimization, and Platform Adaptation”.
  * Chen and Yu [2007] Z. Chen and C. C. Yu, Measurement-noise maximum as a signature of a phase transition, Phys. Rev. Lett. 98, 057204 (2007).
  * Mermin and Wagner [1966] N. D. Mermin and H. Wagner, Absence of ferromagnetism or antiferromagnetism in one- or two-dimensional isotropic heisenberg models, Phys. Rev. Lett. 17, 1133 (1966).
  * Hohenberg [1967] P. C. Hohenberg, Existence of long-range order in one and two dimensions, Phys. Rev. 158, 383 (1967).
  * Clarke and Braginski [2004] J. Clarke and A. I. Braginski,  _The SQUID handbook_ , Vol. 1 (Wiley-VCH, 2004).


# Learning Dense Reward with Temporal Variant Self-Supervision

Yuning Wu1,2, Jieliang Luo2, Hui Li2 1Carnegie Mellon University, Pittsburgh.
2Autodesk Research, San Francisco. This research was conducted during Yuning
Wu’s internship at the Autodesk AI Lab and Autodesk Robotics Lab.

###### Abstract

Rewards play an essential role in reinforcement learning. In contrast to rule-
based game environments with well-defined reward functions, complex real-world
robotic applications, such as contact-rich manipulation, lack explicit and
informative descriptions that can directly be used as a reward. Previous
effort has shown that it is possible to algorithmically extract dense rewards
directly from multimodal observations. In this paper, we aim to extend this
effort by proposing a more efficient and robust way of sampling and learning.
In particular, our sampling approach utilizes temporal variance to simulate
the fluctuating state and action distribution of a manipulation task. We then
proposed a network architecture for self-supervised learning to better
incorporate temporal information in latent representations. We tested our
approach in two experimental setups, namely joint-assembly and door-opening.
Preliminary results show that our approach is effective and efficient in
learning dense rewards, and the learned rewards lead to faster convergence
than baselines.

##  I Introduction

Reinforcement learning (RL) is gaining momentum in solving complex real-world
robotics problems. One challenging category is contact-rich manipulation
tasks. The success of RL in these scenarios depends on a reliable reward
system. While this genre of problems is marked by rich, high-dimensional,
continuous observations, it is typically hard to come up with a dense reward
that can harness such richness to guide RL training. The conventional way of
using sparse, boolean rewards (e.g., 1 if the task is successfully completed
and 0 otherwise) is often challenging and inefficient. The difficulty has led
to the practice of reward engineering, where rewards are hand-crafted based on
domain knowledge and trial-and-error. However, such solutions often require
extensive robotics expertise and can be quite task-specific.

In this research, we propose an end-to-end learning framework that can extract
dense rewards from multimodal observations, inspired by [1]. The reward is
learned in a self-supervised manner by combining one or two human
demonstrations with a physics simulator, and can then be directly used in
training RL algorithms. We evaluate our framework in two contact-rich
manipulation tasks, joint-assembly and door-opening tasks.

There are two main contributions in this paper: 1) a temporal variant forward
sampling (TVFS) method that is more robust and cost-efficient in generating
samples from human demonstrations for contact-rich manipulation tasks, 2) a
self-supervised latent representation learning architecture that can utilize
sample pairs from TVFS.

##  II Problem Statement & Related Work

###  II-A Problem Statement

We focus on contact-rich tasks that can be suitably framed as discrete-time
Markov Decision Processes (MDPs) [2], which is described by a set of states
$S$, a set of actions $A$, a set of conditional probabilities
$p(s_{t+1}|s_{t},a_{t})$ for the state transition $s_{t}\rightarrow s_{t+1}$,
a reward function $R:S\times A\rightarrow\mathbb{R}$, and a discount factor
$\gamma\in[0,1]$. The MDPs can be solved by using RL algorithms to train an
optimal policy $\pi(s)\rightarrow a$ that maximizes the expected total reward.
Our goal is to learn a dense reward function $R$, which can be used by RL
algorithms to reach the optimal policy for the MDPs.

###  II-B Inverse Reinforcement Learning

To tackle the reward engineering problem, Inverse Reinforcement Learning (IRL)
has arisen as a prominent solution [3]. Rather than crafting the reward, IRL
methods learn reward functions [4], [5], [6] from expert demonstrations.
However, conventional IRL methods mostly deal with ideal scenarios where
states and representations are discrete and low-dimensional. Recent advances
in deep learning have extended classical IRL’s capability to continuous, high-
dimensional observation space [7], [8]. However, despite much improved
performance, learning is often conducted using generative adversarial learning
frameworks [9], [10], meaning that one must train the reward function
alongside a policy. Besides instability in training, this framework
essentially diverges from our goal, i.e. to learn a reward function
independently without concurrently learning a policy.

###  II-C Learning Dense Reward for Contact-Rich Manipulation

[11] and [12] explored the idea of training a dense reward function directly
from human feedback. The methods integrated human experts in the RL training
process and periodically ask their preference on a group of pairwise videos
clips of the agent’s behavior. A reward function is gradually trained that
eventually can best explain the human’s judgments. The methods showed
impressive results on training complex robotic locomotion tasks, but haven’t
been tested on contact-rich manipulation tasks where the behaviors are hard to
be observed by pure visual cues. [1] proposed a DREM framework that extracts
dense reward from multimodal observation through sampling and self-supervised
learning. The framework shows great potential in translating rich, continuous,
high-dimensional observations into a task progress variable that can be used
to guide RL training. Our work builds on top of this research, by proposing
improvements to the sampling method and self-supervised learning architecture.
[13] also proposed to learn a multimodal representation of sensor inputs and
use the representation for policy learning on contact-rich tasks, such as peg
insertion. However, it uses crafted reward functions for various sub-tasks.

##  III Approach: Learning Dense Reward with Temporal Variant Self-Supervision

Similar to ideas presented in [1], our method also aims to learn a task
progress variable $p\in[0,1]$ that captures the progress towards finishing a
task. The variable can then be used as a dense reward. With $p=0$ representing
the initial state, and $p=1$ representing the goal state, the variable is
structured as a similarity score in the latent space $\mathcal{H}$. The latent
representation $h_{\phi}:\mathcal{S}\rightarrow\mathcal{H}$ is learned in a
self-supervised manner with two major objectives. The first is to capture an
efficient, low-dimensional embedding of the multimodal observation space
$\mathcal{S}$. The second is to encode temporal information in the learned
representation. Adopted from [1], for contact-rich manipulation tasks with
relatively determined and repeatable goal state, the task progress can be
derived with distance measure $d$ in $\mathcal{H}$:

$\displaystyle p$
$\displaystyle=1-\frac{d(h_{\phi}(s),h_{\phi}(s_{g}))}{d(h_{\phi}(s_{0}),h_{\phi}(s_{g}))}$
(1)

Prior work has explored ways to learn the representation through explicitly
enforcing temporal order through a triplet loss function [1]. Such enforcement
by design involves tuning multiple hyperparameters. We propose a framework
where temporal information can be injected in a more natural, self-consistent
manner by utilizing dynamic relation among pairs of adjacent observations
$(s_{t},s_{t+1})$.

$h_{\phi}(s_{t})+\Delta h_{\psi}(s_{t})=h_{\phi}(s_{t+1})$ (2)

$h_{\phi}$ is the latent representation, whereas $\Delta h_{\psi}$ is the
change of latent representation between $t$ and $t+1$ resulted from dynamics.
$h_{\phi}$ and $\Delta h_{\psi}$ are learned using different modalities. The
insight behind such constraint is that latent representation of $s_{t+1}$
should be consistent with the representation of $s_{t}$, plus any dynamic
change happened within the time step.

Our proposed improvements are the following. The first is temporal variant
forward sampling, which generates a tree of data (i.e. observations with
temporal information) from a single human demonstration. The second is self-
supervised representation learning network architecture, which uses generated
observation pairs to learn representations through Eq.(2).

###  III-A Temporal Variant Forward Sampling

Figure 1: Sampling variance along different stages of a manipulation task. The
blue $V_{1}(t)$ and red $V_{2}(t)$ curves show potential temporal variance
control functions that can be used in sampling.

Collecting observations with temporal information is an essential step for
training, but can be challenging if only from human demonstrations. Therefore,
the idea of sampling has been broadly experimented in [1], [14]. By combining
one or two human demonstrations with a physics simulator, it is possible to
obtain a tree of data through sampling. [1] proposed a backward sampling
process based on the insight that variance of the goal state is smaller than
the initial state. However, although it is feasible to sample backward
positions and generate visual images from the positions in simulation, it is
typically hard to sample backward force and torque (F/T). Through experiments,
we have found that restoring a state with the exact F/T reading can be
computationally intensive and simulator-dependent. Also, when playing forward
a backward sampled action sequence, the F/T readings in the forward pass do
not necessarily match the F/T readings recorded in the backward pass.

Without loss of generality, we propose a new sampling process, named Temporal
Variant Forward Sampling (TVFS) that aims to tackle the aforementioned
challenges, while capturing the fluctuating variance of manipulation tasks.
The insight behind our method is to roughly control sampled actions with a
temporal variance $V(t)$, such that sampled actions do not diverge too much
from the potential distribution of an expert demonstration, and that the
actions are mostly progressing forward. For instance, an action that is the
opposite of an expert action may not appear at certain stable stages. As shown
in Fig.1, at the starting stage ($p=0$), the sampling variance is limited. At
the intermediate stage, the variance is high due to lack of constraints and
high moving flexibility. At the ending stage ($p=1$), the variance is low
because the goal state is relatively deterministic. $V(t)$ can be depicted
with a chosen kernel function. In our experiment, we have chosen the quadratic
function for simplicity. The general process of our sampling method is
illustrated in Fig.2, and is described as follows.

Figure 2: Temporal variant forward sampling (TVFS). The difference between
sampled actions (blue) and demo actions (black) are controlled through $V(t)$,
which is also a variance measure that changes along the task process.

  1. 1.

We record an expert demonstration in simulation, and choose the sampling seeds
(states) $\\{Q_{0},Q_{1},\cdots,Q_{M}\\}$ by certain sampling interval.

  2. 2.

At each seed (state) $Q_{i}$, randomly sample $N$ branches. Each branch may
contain multiple forward steps. At each step, control the variance between
sampled action $a^{\text{sampled}}_{t}$ and demo action $a^{\text{demo}}_{t}$.
The variance can be measured using any similarity score. In our case, we have
chosen the cosine similarity.

  3. 3.

Record the sampled observations in pairs $(s_{t},s_{t+1})$, such that they can
be used in learning Eq.(2).

###  III-B Multimodal Representation Learning

With the generated multimodal observation pairs from TVFS, we have designed a
network architecture and a loss function to incorporate temporal information
in representation learning. As mentioned above, we use different modalities to
learn different components of Eq.(2). $h_{\phi}$ is learned with static
modalities such as images, depth maps and poses, while $\Delta h_{\psi}$ is
learned using dynamic modalities such as F/T and velocities. The two
separately learned components should be consistent in the latent space. We
accentuate such consistency with a hybrid loss function, consisting of
temporal enforcement loss and reconstruction loss. The architecture and loss
functions are detailed as follows.

Figure 3: Self-supervised learning network architecture

  1. 1.

Static Modality Encoder $\mathbf{E^{s}}$. To learn $h_{\phi}$, we use a fixed
RGB-D camera as input for static modality encoders. The RGB image (256x256x3)
and depth map (256x256x1) are handled separately. Similar to [13], the network
is composed of a 6-layer Convolutionary Neural Network and a fsully-connected
layer. Depending on the experiment scenarios, one may switch or combine the
modalities. An extra Multi-Layer Perceptron may for output fusion. The final
embedding is a 64 dimensional hidden vector.

  2. 2.

Dynamic Modality Encoder $\mathbf{E^{d}}$. To learn $\Delta h_{\psi}$, we use
F/T reading and velocity as input for the dynamic modality encoder. Due to the
accumulative nature of F/T, we use a window size of 32 to better capture the
momentum. The 32x6 input is passed into a 4-layer Causal Convolution Network.
The output is concatenated and fused with velocity to produce another 64
dimensional hidden vector.

  3. 3.

Static Modality Decoder $\mathbf{D^{s}}$. Through experiments, we have found
that instead of enforcing self-supervised learning on both static and dynamic
modalities, it is better to focus on one side only. This choice will be
explained further in later descriptions. In our case, we are proposing an
auto-encoding architecture on the static modality side. The decoder takes
input of a 64 dimensional hidden vector and use transposed Convolutional
Neural Network to reconstruct the RGB image / depth map.

  4. 4.

Latent Representation Learning. As mentioned briefly in previous context, the
temporal order is injected through learning with pairs of adjacent
observations $(s_{t},s_{t+1})$. By enforcing Eq.(2) among each pair, the
temporal relation among latent representations are broadcasted. Fig.3 is an
illustration of the whole network architecture. We first encode and decode
$s_{t}$ with $\mathbf{E^{s}}$, $\mathbf{E^{d}}$ and $\mathbf{D^{s}}$, then use
the embedding to craft an latent representation for $s_{t+1}$.

$\displaystyle\hat{h}(s_{t+1})=\mathbf{E^{s}}(s_{t})+\mathbf{\bar{E}^{d}}(s_{t})$
(3)

The hybrid loss function is structured around $\hat{h}(s_{t+1})$. The first
component is temporal enforcement loss, which enforces Eq.(2) in the latent
space. To ensure effectiveness of $\Delta h_{\psi}$, we are applying $L2$
normalization to $\mathbf{\bar{E}^{d}}(s_{t})$.

$\displaystyle l^{\text{temporal}}$
$\displaystyle=\text{MSE}\left[\hat{h}(s_{t+1}),\mathbf{E^{s}}(s_{t+1})\right]$
(4)

The second component is reconstruction loss, which provides supervision signal
to the auto-encoding architecture. As apposed to directly decoding
$\mathbf{E^{s}}(s_{t+1})$, we are decoding $\hat{h}(s_{t+1})$ so that Eq.(2)
is also enforced in self-supervision.

$\displaystyle l^{\text{recon}}$
$\displaystyle=l^{\text{recon}}(s_{t})+\hat{l}^{\text{recon}}(s_{t+1})$ (5)
$\displaystyle l^{\text{recon}}(s_{t})$
$\displaystyle=\text{MSE}\left[\mathbf{D^{s}}(\mathbf{E^{s}}(s_{t})),s_{t}\right]$
(6) $\displaystyle\hat{l}^{\text{recon}}(s_{t+1})$
$\displaystyle=\text{MSE}\left[\mathbf{D^{s}}\left(\hat{h}(s_{t+1})\right),s_{t+1}\right]$
(7)

The two loss components are then combined through a hyperparameter $\lambda$.

$\displaystyle l=l^{\text{recon}}+\lambda\cdot l^{\text{temporal}}$ (8)

We set $\lambda=10$ in training to accentuate the temporal relation, so that
the representation learning does not converge suboptimally too early. The
learned embedding is then used in Eq.(1) for the dense reward.

##  IV Experimental Results

The experiments are conducted in simulation. In order to test that our
sampling method can be generalized to different simulators and robot
controllers, we tested lap-joint assembly in PyBullet [15] with an robot-
agnostic environment, and door-opening in Robosuite [16] with a Panda robot.
We conducted similar sampling process on both tasks, setting sampling interval
$I=50$, number of branches $N=5$, number of steps per branch $K=10$. The
temporal variance is controlled in $\theta\in[\frac{\pi}{12},\frac{\pi}{4}]$.
We trained the model with an NVIDIA 3060 GPU for around 5,000 iterations.
Compared to the training iterations mentioned in [1], our method is
potentially more efficient. We defer ablation study and further examination of
this comparison to future work.

###  IV-A Visualization of Learned Dense Reward

We visualized the learned dense reward in two cases. The results indicate that
the dense reward learned by our approach is effective. In the first case
(Fig.4), we compare the rewards between a successful trajectory and a failed
trajectory in the lap-joint task. The plots suggest that a successful
trajectory has rewards gradually increasing from 0 to 1, which matches the
definition of task progress. A failed trajectory have a decreasing reward
dropping below 0, which can happen when the agent get into unexpected
scenarios.

Figure 4: Comparison of a successful trajectory (left) and a failed trajectory
(right). Visualization produced in the lap-joint task.

In the second case (Fig.5), we examine rewards of an inexpert demonstration in
door-opening. The demonstrator experienced a plateau of trial-and-error when
rotating back and forth the door handle. While the hand-crafted reward mostly
gives analogous signals during this period, our rewards provides fluctuations
indicating more detailed feedback for learning.

Figure 5: Comparison of dense reward and hand-crafted rewards in door opening
task

###  IV-B Dense Rewards for Policy Training

To examine the performance of the learned reward in policy training, we have
chosen Soft Actor-Critic [17] as the RL algorithm for benchmarking. We
compared three types of rewards in the door-opening task, namely our dense
reward, a hand-crafted reward based on distance ($\gamma\|x_{t}-x_{g}\|_{2}$),
and the sparse boolean reward. In particular, we trained the policy for door-
opening task for 500 epochs. For each type of reward, we conducted three
training experiments with different random seeds. The results (Fig.6) indicate
that our dense reward leads to faster convergence, and more training
stability.

Figure 6: RL training comparison among three types of rewards: our dense
reward (blue), hand-crafted distance reward (orange), and sparse boolean
reward (red).

##  V Conclusions and Future Works

In this paper, we propose an improved framework for learning dense reward for
contact-rich manipulation tasks. The framework includes a more robust sampling
method, namely temporal variant forward sampling (TVFS), that can generate
samples from one or two human demonstrations with a physics simulator. A self-
supervised learning architecture is also designed to efficiently utilize the
generated sample pairs.

For future work, we intend to conduct more ablation studies regarding the
framework’s adaptability and modalities. For instance, during experiments we
observe that camera setup can have a substantial impact on the learning
result. Therefore one potential is to study whether we can mainly rely on pure
tactile sensors for reward inference. Another potential is to test whether the
reward can be transferred to manipulation tasks with nondeterministic goal
state.

## Acknowledgement

We thank Tonya Custis and Sachin Chitta for budgetary support of the project;
Yotto Koga for simulation support; our colleagues at Autodesk Research for the
valuable feedback, and Zheng Wu for the discussions.

## References

  * [1] Z. Wu, W. Lian, V. Unhelkar, M. Tomizuka, and S. Schaal, “Learning dense rewards for contact-rich manipulation tasks,” _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , pp. 6214–6221, 2021.
  * [2] M. L. Puterman, _Markov decision processes: discrete stochastic dynamic programming_. John Wiley & Sons, 2014\.
  * [3] A. Y. Ng, S. J. Russell, _et al._ , “Algorithms for inverse reinforcement learning.” in _Icml_ , vol. 1, 2000, p. 2.
  * [4] P. Abbeel and A. Y. Ng, “Apprenticeship learning via inverse reinforcement learning,” in _Proceedings of the twenty-first international conference on Machine learning_ , 2004, p. 1.
  * [5] B. D. Ziebart, A. L. Maas, J. A. Bagnell, A. K. Dey, _et al._ , “Maximum entropy inverse reinforcement learning.” in _Aaai_ , vol. 8. Chicago, IL, USA, 2008, pp. 1433–1438.
  * [6] D. Ramachandran and E. Amir, “Bayesian inverse reinforcement learning.” in _IJCAI_ , vol. 7, 2007, pp. 2586–2591.
  * [7] M. Wulfmeier, P. Ondruska, and I. Posner, “Maximum entropy deep inverse reinforcement learning,” _arXiv preprint arXiv:1507.04888_ , 2015.
  * [8] S. K. Seyed Ghasemipour, S. S. Gu, and R. Zemel, “Smile: Scalable meta inverse reinforcement learning through context-conditional policies,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.
  * [9] J. Ho and S. Ermon, “Generative adversarial imitation learning,” _Advances in neural information processing systems_ , vol. 29, 2016.
  * [10] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial nets,” _Advances in neural information processing systems_ , vol. 27, 2014.
  * [11] P. F. Christiano, J. Leike, T. Brown, M. Martic, S. Legg, and D. Amodei, “Deep reinforcement learning from human preferences,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [12] K. Lee, L. Smith, and P. Abbeel, “Pebble: Feedback-efficient interactive reinforcement learning via relabeling experience and unsupervised pre-training,” _arXiv preprint arXiv:2106.05091_ , 2021.
  * [13] M. A. Lee, Y. Zhu, K. Srinivasan, P. Shah, S. Savarese, L. Fei-Fei, A. Garg, and J. Bohg, “Making sense of vision and touch: Self-supervised learning of multimodal representations for contact-rich tasks,” in _2019 International Conference on Robotics and Automation (ICRA)_. IEEE, 2019, pp. 8943–8950.
  * [14] C. Florensa, D. Held, M. Wulfmeier, M. Zhang, and P. Abbeel, “Reverse curriculum generation for reinforcement learning,” in _Conference on robot learning_. PMLR, 2017, pp. 482–495.
  * [15] E. Coumans and Y. Bai, “Pybullet, a python module for physics simulation for games, robotics and machine learning,” http://pybullet.org, 2016–2021.
  * [16] Y. Zhu, J. Wong, A. Mandlekar, and R. Martín-Martín, “robosuite: A modular simulation framework and benchmark for robot learning,” in _arXiv preprint arXiv:2009.12293_ , 2020.
  * [17] T. Haarnoja, A. Zhou, P. Abbeel, and S. Levine, “Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor,” in _International conference on machine learning_. PMLR, 2018, pp. 1861–1870.