# Magnetic properties and field-driven dynamics of chiral domain walls
in epitaxial Pt/Co/AuxPt1-x trilayers

Kowsar Shahbazi  School of Physics and Astronomy, University of Leeds, Leeds
LS2 9JT, United Kingdom  Aleš Hrabec  Present Address: Mesoscopic Systems,
Paul Scherrer Institut, ODRA/113, 5232 Villigen, Switzerland.  School of
Physics and Astronomy, University of Leeds, Leeds LS2 9JT, United Kingdom
Laboratoire de Physique des Solides, CNRS, Université Paris-Sud, Université
Paris-Saclay, 91405 Orsay Cedex, France  Simone Moretti  Departamento Fisica
Aplicada, University of Salamanca, Plaza de los Caidos s/n E-37008, Salamanca,
Spain  Fachbereich Physik, Universität Konstanz, 78457 Konstanz, Germany
Michael B. Ward  Leeds Electron Microscopy and Spectroscopy Centre, School of
Chemical and Process Engineering, University of Leeds, Leeds LS2 9JT, United
Kingdom  Thomas A. Moore  School of Physics and Astronomy, University of
Leeds, Leeds LS2 9JT, United Kingdom  Vincent Jeudy  Laboratoire de Physique
des Solides, CNRS, Université Paris-Sud, Université Paris-Saclay, 91405 Orsay
Cedex, France  Eduardo Martinez  Departamento Fisica Aplicada, University of
Salamanca, Plaza de los Caidos s/n E-37008, Salamanca, Spain  Christopher H.
Marrows<EMAIL_ADDRESS>School of Physics and Astronomy, University
of Leeds, Leeds LS2 9JT, United Kingdom

###### Abstract

Chiral domain walls in ultrathin perpendicularly magnetised layers have a Néel
structure stabilised by a Dzyaloshinskii-Moriya interaction (DMI) that is
generated at the interface between the ferromagnet and a heavy metal.
Different interface materials or properties are required above and below a
ferromagnetic film in order to generate the structural inversion asymmetry
needed to ensure that the DMI arising at the two interfaces does not cancel.
Here we report on the magnetic properties of epitaxial Pt/Co/AuxPt1-x
trilayers grown by sputtering onto sapphire substrates with 0.6 nm thick Co.
As $x$ rises from 0 to 1 a structural inversion asymmetry is progressively
generated. We characterise the epilayer structure with x-ray diffraction and
cross-sectional transmission electron microscopy, revealing (111) stacking.
The saturation magnetization falls as the proximity magnetisation in Pt is
reduced, whilst the perpendicular magnetic anisotropy $K_{\mathrm{u}}$ rises.
The micromagnetic DMI strength $D$ was determined using the bubble expansion
technique and also rises from a negligible value when $x=0$ to $\sim 1$ mJ/m2
for $x=1$. The depinning field at which field-driven domain wall motion
crosses from the creep to the depinning regime rises from $\sim 40$ to $\sim
70$ mT, attributed to greater spatial fluctuations of the domain wall energy
with increasing Au concentration. Meanwhile, the increase in DMI causes the
Walker field to rise from $\sim 10$ to $\sim 280$ mT, meaning that only in the
$x=1$ sample is the steady flow regime accessible. The full dependence of
domain wall velocity on driving field bears little resemblance to the
prediction of a simple one-dimensional model, but can be described very well
using micromagnetic simulations with a realistic model of disorder. These
reveal a rise in Gilbert damping as $x$ increases.

###### pacs:

75.70.Tj; 75.76.+j; 75.50.-y

##  I Introduction

Research into magnetic thin films has made huge progress over the last few
decades, largely due to the additional functionality emerging at interfaces
between a magnetic film of few or sub-nanometer thickness and its
surroundings. In addition to a wide range of novel magnetotransport effects
Hoffmann and Bader (2015), the magnetic properties are profoundly affected by
the choice of materials with which those interfaces are made and the nature
and quality of the interfaces. Emergent properties span both the statics and
dynamics of the magnetization. One of the first to be discovered was the
perpendicular magnetic anisotropy (PMA) de Jonge  _et al._ (1994), which can
be strong enough to cause a spin-reorientation transition and lift the
magnetisation out of the film plane in the ground state. Forming an interface
with a metal close to satisfying the Stoner criterion can lead to proximity
magnetism Antel  _et al._ (1999), which can account for a significant
proportion of the total moment if the magnetic layer is very thin. Long-
predicted Fert (1991), but only recently observed Chen  _et al._ (2013), is
the interfacial Dzyaloshinskii-Moriya interaction (DMI), which arises due to
inversion symmetry breaking and favours chiral magnetic states Benitez  _et
al._ (2015); Tetienne  _et al._ (2015). Meanwhile, the damping of
magnetisation dynamics can be enhanced by spin-pumping into the surrounding
layers Tserkovnyak  _et al._ (2002a).

Magnetising a film perpendicular to its plane is a highly demagnetising
configuration and domain structures Hubert and Schäfer (1998), separated by
domain walls (DWs), commonly occur. Whilst magnetostatics alone favours a
Bloch wall configuration, the interfacial DMI takes a form that would prefer a
Néel structure with a fixed chirality Thiaville  _et al._ (2012). This is
important for the latest forms of racetrack memory that exploit PMA materials
Parkin and Yang (2015), since the DWs are driven by spin-orbit torques Manchon
(2014) that require a component of magnetization collinear to the electron
flow. This is not trivial to arrange by other means and it explains, for
example, why Bloch walls are insensitive to this type of torque Khvalkovskiy
_et al._ (2013), and why the wall chirality determines the direction of wall
motion under current Emori  _et al._ (2013); Ryu  _et al._ (2014).

Many of these effects, for instance the PMA and DMI, require strong spin-orbit
coupling (SOC) and so it is heavy metals that are commonly used to form the
interfaces with the ultrathin magnetic layer, usually in polycrystaline
multilayers prepared by sputtering. Nevertheless, properties such as the DMI
are extremely sensitive to the interface quality Hrabec  _et al._ (2014);
Wells  _et al._ (2017). To avoid such ambiguities and in order to explore this
systematically one has to grow such materials in a controlled way. Here we
experimentally investigate epitaxial layers of Pt/Co/AuxPt1-x prepared by
sputtering onto single crystal sapphire substrates Mihai  _et al._ (2013). Pt
is a common pairing with Co that has strong PMA and DMIHrabec  _et al._
(2014), and is close to satisfying the Stoner criterionStoner (1938), giving
strong proximity magnetism. The choice of Au, a heavy element with fully
filled $5d$ band, is motivated by a negligible proximity effect Wilhelm  _et
al._ (2004); Ryu  _et al._ (2014) and the expectation of a small DMI arising
at Co/Au interface Kashid  _et al._ (2014); Yang  _et al._ (2015). Gold
therefore serves as a textbook element giving the opportunity to study the
effect of controlled broken symmetry on either side of the ferromagnet and its
importance for the various effects arising at these interfaces.

##  II Epilayer growth and characterization

Figure 1: Structural characterization of Pt/Co/AuxPt1-x trilayers. (a) X-ray
diffraction spectra confirming the fcc (111) orientation of the samples. (b)
Dark field HRTEM cross-section images. (c) EDX spectra showing distinguishable
layers with possible intermixing at the interfaces.

The trilayers of Pt(3 nm)/Co(0.6 nm)/AuxPt1-x (3 nm), in which $0\leq x\leq
1$, were prepared by sputtering at high temperatures, as described in Ref.
Mihai  _et al._ , 2013. The seed Pt layer was sputtered directly on a C-plane
sapphire substrate at $500^{\circ}$C followed by Co at $100^{\circ}$C. The
AuxPt1-x layer was grown by co-sputtering from Pt and Au targets at
$100^{\circ}$C in a way that the sputtering powers were adjusted to keep the
rate $\sim 1$ Å/s. The multilayers were capped with $\sim$2 nm of Ta, which is
polycrystalline.

Figure 2: Magnetic characterization of Pt/Co/AuxPt1-x trilayers. (a) Polar
MOKE loops. (b) Anisotropy field, $H_{\mathrm{K}}$, and effective anisotropy
energy, $K_{\mathrm{eff}}$ both increase with Au concentration. (c) Saturation
magnetization $M_{\mathrm{s}}$ supposing that all the magnetic moment is
confined in the Co layer (black dots), extracted proximity-induced magnetism
in Pt layers (red) and corrected Co layer values $M_{\mathrm{Co}}$ (blue). (d)
Variation of relative saturation magnetization with temperature for Pt/Co/Pt.
The solid red line is a fit of the Bloch law.

The quality of these epilayers was confirmed by x-ray diffraction (XRD) and
transmission electron microscopy (TEM). The XRD spectra for $x=0$, 0.5, and 1,
shown in Fig. 1(a), are consistent with the previous study Mihai  _et al._
(2013). In each case, a broad principal peak (due to the extreme thinness of
the layers) is found close to $\theta=40^{\circ}$, surrounded by satellite
Pendellösung fringes. This peak is very close to bulk Pt (111) Bragg peak at
$\theta=39.74^{\circ}$, showing the matching of the multilayer structure with
the substrate. The appearance of Pendellösung fringes in the XRD of a
multilayer emphasise the uniformity of the strain in the layers and smoothness
of the interfaces Ploog  _et al._ (2013); Tanner (2013). The peak moves to
slightly lower angles as Au is introduced into the upper layer, as might be
expected from a Vegard law-type behaviour.

Cross-sections of the three samples are shown as dark field high-resolution
transmission electron microscopy (HRTEM) images in Fig. 1(b). Lattice fringes
are visible extending coherently across all three layers, emphasising the
epitaxial nature of the samples. Energy dispersive x-ray analyses of each of
the three samples are shown in Fig. 1(c), showing distinct layers. The lattice
spacing in the vertical direction, $d_{111}$, can be evaluated using position
of the main peak in the XRD pattern and Bragg’s law, or by averaging the
lattice spacings in the cross-sectional TEM images. Comparing evaluated values
from both measurements, there is just about 2% difference in $d_{111}$ between
them (within error bars) which is a good agreement.

The contrast in TEM images is related to the atomic number inducing a darker
appearance for the Co layer in comparison with Pt and Au. Investigating the
changes of greyness in vertical direction of images results in Co thicknesses
of 1.03 nm, 1.26 nm and 1.17 nm for $x=1$, 0.5 and 0, respectively. (The
thickness was measured for stripes of $\sim 0.4$ nm width across the entire
image, and then the average was calculated). This is much higher than the
nominal thickness of 0.6 nm. In fact, this is expected as TEM shows the
additive effect of all atomic layers in the depth of FIB-cut lamella (which is
usually $\sim 50$ nm). So, when the electron beam goes through the TEM
specimen, the thickness fluctuations of the thin layer will be added to the
actual thickness resulting in the overestimation of thickness. It has also
been reported previously than evaluating thickness from TEM cross-sections
might result in overestimated thicknesses, sometimes even two times higher
than the real value Gabureac  _et al._ (2014). Hence, x-ray reflectivity
patterns of samples were fitted using the GenX package Björck and Andersson
(2007), which yielded thicknesses very close to the nominal ones for the Co
layers. Considering all these facts, the nominal thickness of 0.6 nm was taken
to be a good estimate of the true thickness value for further analysis of
parameters describing the samples.

All the films showed square magneto-optical Kerr effect (MOKE) hysteresis
loops as a function of perpendicular magnetic field, confirming that the PMA
generated at the Co-heavy metal interfaces is strong enough to overcome the
shape anisotropy (Fig. 2(a)). Superconducting quantum interference device
vibrating sample magnetometry (SQUID-VSM) was used to measure the anisotropy
field $H_{\mathrm{K}}$, inferred from the in-plane saturation field, plotted
in Fig. 2(b), and the saturation magnetization $M_{\mathrm{s}}$ shown in Fig.
2(c). To calculate $M_{\mathrm{s}}$, it was assumed that all the measured
moment is confined in the Co layer of thickness $t$. However, due to the
proximity effects, induced magnetic moments also exist in the Pt layers. The
induced moment in Pt can be estimated by assuming a negligible proximity
effect in Au Wilhelm  _et al._ (2004), and so the difference in moment between
the $x=0$ and $x=1$ films corresponds to the induced moment in a single Pt
layer. This is summarized in Fig. 2(c) which reveals the contribution of the
magnetization in the Pt, to the constant magnetization of the Co layer
($1.0\pm 0.1$ MA/m). This is comparable to the previously published values for
Co/Pt multilayers with the same nominal thickness of cobalt Lin  _et al._
(1991); Bandiera  _et al._ (2012); Lavrijsen  _et al._ (2015); Metaxas  _et
al._ (2007); Hrabec  _et al._ (2014); Miron  _et al._ (2010). On the other
hand, if we consider thicknesses deduced from the TEM measurements, the
saturation magnetization of mulilayers would be about 0.6-0.7 MA/m which is
unusually low for supposedly 1 nm of Co Lin  _et al._ (1991); Vaňatka  _et
al._ (2015). This also emphasises the unreliability of thickness estimations
from TEM cross-sections. Following this analysis, the effective anisotropy can
be calculated by
$K_{\mathrm{eff}}=\frac{1}{2}\mu_{0}M_{\mathrm{Co}}H_{\mathrm{K}}$ that is
also plotted in Fig. 2(b), which rises markedly with Au content $x$. The
values, along with the other magnetic data for each sample, are given in Table
1. The temperature dependence of $M_{\mathrm{Co}}$ between 5 and 400 K for the
$x=0$ sample shown in Fig. 2(d) was fitted with the Bloch lawChikazumi (1997),

$\frac{M_{\mathrm{Co}(T)}}{M_{0}}=1-c\left(\frac{k_{\mathrm{B}}T}{J}\right)^{\frac{3}{2}},$
(1)

where $M_{0}$ is the magnetization at 0 K, $c=0.0294$ for an fcc lattice,
$k_{\mathrm{B}}$ is the Boltzmann constant and $J$ is the exchange integral.
Knowing $J$, one can obtain the exchange stiffness as $A=QJS^{2}/a$, with
$Q=4$ as the number of atoms per unit cell for fcc, $S$ the spin quantum
number, and $a=0.355$ nm the lattice constant. The value we obtain is $A=17\pm
1$ pJ/m. This reduction in $A$ with respect to the bulk value for Co of 30
pJ/m is in good agreement with our previous findings for Co layers of this
thickness Shepley  _et al._ (2018).

Table 1: Magnetic properties of epitaxial Pt/Co/AuxPt1-x for the three different values of Au concentration $x$ studied here. These are saturation magnetization $M_{\mathrm{CO}}$, uniaxial anisotropy constant $K_{\mathrm{u}}=K_{\mathrm{eff}}+\frac{1}{2}\mu_{0}M_{\mathrm{Co}}^{2}$, DW width $\Delta=\sqrt{A/K_{\mathrm{eff}}}$ , Gilbert damping constant from creep and depinning fits $\alpha_{\mathrm{exp}}$, Walker breakdown field $\mu_{0}H_{\mathrm{W}}$ calculated using DW mobility, Gilbert damping constant from micromagnetics $\alpha_{\mathrm{\upmu M}}$, depinning field $H_{\mathrm{d}}$ obtained from fitting with the creep law, thickness fluctuation $\delta$ used in the micromagnetics simulations, correlation length of the disorder $\xi$, strength of the pinning disorder $f_{\mathrm{pin}}/\xi$, and effective DMI constant $D_{\mathrm{eff}}$, as well as the inferred handedness of the chirality of the domain walls. $x$ | $M_{\mathrm{Co}}$ | $K_{\mathrm{u}}$ | $\Delta$ | $\alpha_{\mathrm{exp}}$ | $\mu_{0}H_{\mathrm{W}}$ | $\alpha_{\mathrm{\upmu M}}$ | $\mu_{0}H_{\mathrm{d}}$ | $\delta$ | $\xi$ | $f_{\mathrm{pin}}/\xi$ | $D_{\mathrm{eff}}$ | Chirality
---|---|---|---|---|---|---|---|---|---|---|---|---
| MA/m | MJ/m3 | nm |  | mT |  | mT |  | nm | meV | mJ/m2 |
1 | $1.0\pm 0.1$ | $1.2\pm 0.2$ | $5.5\pm 0.9$ | $0.9\pm 0.1$ | $280\pm 10$ | $0.35\pm 0.05$ | $72\pm 2$ | $9\%$ | $25\pm 3$ | $510\pm 50$ | $-1.0\pm 0.3$ | Left
0.5 | $1.1\pm 0.2$ | $1.1\pm 0.2$ | $7\pm 1$ | $0.4\pm 0.1$ | $49\pm 3$ | $0.35\pm 0.05$ | $56\pm 2$ | $8.5\%$ | $21\pm 2$ | $330\pm 30$ | $-0.35\pm 0.09$ | Left
0 | $1.0\pm 0.2$ | $1.0\pm 0.2$ | $8\pm 1$ | $0.2\pm 0.1$ | $11\pm 1$ | $0.17\pm 0.01$ | $39\pm 2$ | $7\%$ | $22\pm 2$ | $240\pm 20$ | $-0.07\pm 0.07$ | None

##  III Domain wall chirality

Ferromagnet-heavy metal interfaces also exhibit interfacial DMIFert (1991),
which enforces the Néel form of the DW profile with a fixed chiralityThiaville
_et al._ (2012); Benitez  _et al._ (2015). Here we have studied this property
using the asymmetric bubble expansion methodJe  _et al._ (2013); Hrabec  _et
al._ (2014). A bubble domain is nucleated and then expanded by an out-of-plane
applied field $H_{z}$, whilst $H_{\mathrm{DMI}}$, the effective field induced
at the DW by the DMIThiaville  _et al._ (2012) is enhanced or compensated by
an in-plane applied field $H_{x}$, resulting in different creep-regime
velocities for walls on either side of the bubble moving with or against the
in-plane field Kabanov  _et al._ (2010); Je  _et al._ (2013). The effective
micromagnetic DMI strength is given by
$D_{\mathrm{eff}}=\mu_{0}H_{\mathrm{DMI}}M_{\mathrm{Co}}\Delta$, where
$\Delta$ is the DW width. To measure $H_{\mathrm{DMI}}$ we have used Kerr
microscopy with a two-coil setup to apply $H_{z}$ and $H_{x}$ independently.
We measure the DW creep velocity as a function of $H_{x}$ by pulsing $H_{z}$
with a fixed amplitude. The DW motion at low magnetic fields follows the creep
law Metaxas  _et al._ (2007), and its velocity can be expressed as

$v=v_{0}\exp\left(\zeta H_{z}^{-\frac{1}{4}}\right)$ (2)

where $v_{0}$ is the velocity prefactor, $H_{z}$ is out-of-plane driving
field, and $\zeta$ is the scaling coefficient

$\zeta=\zeta_{0}\left[\frac{\sigma(H_{x})}{\sigma(0)}\right]^{\frac{1}{4}}$
(3)

where $\zeta_{0}$ can be considered as a scaling constant for pinning
properties independent of in-plane field. $\sigma$ is the Néel wall energy
density Je  _et al._ (2013); Hrabec  _et al._ (2014), which contains a term
proportional to $|H_{x}+H_{\mathrm{DMI}}|$. Hence the highest wall energy, and
(through Eqs 2 and 3) minimum DW creep velocity $v$ corresponds to the point
where $\mu_{0}H_{x}=-\mu_{0}H_{\mathrm{DMI}}$, i.e. the point at which
external field transforms the Néel wall to a Bloch wallHrabec  _et al._
(2014).

We performed bubble expansion measurements on all three epilayers, the results
are shown in Fig. 3. Measurements as a function of $H_{z}$ at given values of
$H_{x}$ (not shown) yielded a $\ln v\propto H_{z}^{-1/4}$ scaling, confirming
that measurements are in the creep regime. A typical Kerr difference
micrograph of an expanded DW bubble in Pt/Co/Au film is displayed in inset of
Fig. 3(a), showing the asymmetric expansion. Examples of datasets for
$v(H_{x})$ for Pt, Pt50Au50, and Au capping layers are shown in Fig. 3(b),
(c), and (d), respectively. One can immediately see a shift of velocity
minimum away from zero in-plane field, indicating the presence of a DMI. The
solid lines correspond to the fits of Eq. 2. The effective DMI field for each
sample, determined from these fits, is plotted in Fig. 3(a) as a function of
Au concentration $x$. It should be noted that the high uncertainty in the $D$
value for Au capped trilayer is due to not reaching the minimum velocity (i.e.
the inversion symmetry axis of the curve) cannot be seen, as a result of
experimental limitations, and so our value is based just on the fit to the
data that is shown. Negligible DMI in Pt/Co/Pt confirms that the system has an
high degree of inversion symmetry for $x=0$, which is not a trivial result
Hrabec  _et al._ (2014); Lavrijsen  _et al._ (2015), since small differences
in top and bottom interface structure can easily break inversion symmetry and
lead to an appreciable DMIWells  _et al._ (2017). The negligible DMI in this
case can thus be expected to lead to walls of of Bloch form with no fixed
chirality.

Also, Kim et al. stated that this assumption that the velocity minimum occurs
at $-\mu_{0}H_{\mathrm{DMI}}$ only holds if the minimum velocity falls on the
inversion symmetry axis of the velocity curve Kim  _et al._ (2018). In other
words, it holds if the asymmetric contribution to the velocity is negligible.
In such a case, any contribution from chiral damping is negligible and DMI can
be accounted for all the chiral behaviour in the sample. Examination of
symmetric and asymmetric contributions of our velocity curves, showed that the
latter is one order of magnitude smaller than the former, proving the
dominance of the symmetric contribution and ruling out any chiral damping
effect.

The effective DMI field and the strength of DMI constant $D_{\mathrm{eff}}$
increase with $x$. This is expected since the trilayers are becoming more
inversion asymmetric. From the symmetries of the applied fields we deduce that
for $x>0$, the DMI enforces left-handed chirality of the Néel walls, which is
consistent with ab initio calculations for Pt/Co interfaces Yang  _et al._
(2015).

Figure 3: Asymmetric bubble expansion measurements. (a) Effective DMI field
$\mu_{0}H_{\mathrm{DMI}}$ and effective DMI constant $D_{\mathrm{eff}}$ as a
function of Au concentration $x$ in Pt/Co(6 Å)/AuxPt1-x epilayers. An example
of a Kerr difference image of the bubble expansion in Pt/Co/Au is shown as an
inset. The values plotted were derived from the data in panels (b)-(d), which
show DW velocity as a function of in-plane magnetic field $\mu_{0}H_{x}$ in
case of Pt, Pt50Au50 and Au capping layers, respectively, for
$\mu_{0}H_{z}=7$, 8, and 10.5 mT.

##  IV Field-driven domain wall dynamics

In the previous section, we considered some aspects of the creep regime, where
the applied field is small with respect to the depinning field and DW motion
is thermally activated. For applied fields exceeding the depinning field, DW
dynamics first goes through a so-called depinning transition and then move in
the viscous flow regime where temperature plays no role: the DW velocity is
controlled by damping, parameterised by the Gilbert damping coefficient
$\alpha$. We measured the DW velocity across both regimes by investigating
bubble expansion driven by pulsed out-of-plane applied field using a fast
micro-coil of diameter 1 mm that produced pulses with a 200 ns rise time. The
DW velocity is then calculated by dividing the distance which the DW covered
during the pulse by the pulse duration. The velocities as a function of out-
of-plane field $\mu_{0}H_{z}$ are shown in Fig. 4(a). The DW velocity for a
specific value of $H_{z}$ is rising with the broken inversion symmetry of the
trilayer. The measured velocities cover the thermally activated creep regime,
the depinning transition, and the beginning of the appearance of the expected
linear change of velocity in the flow regime. The data were fitted
simultaneously with the universal creep and depinning regime functionsJeudy
_et al._ (2018):

$v(H)=\begin{cases}v(H_{\mathrm{d}})\exp\left(-\frac{\Delta
E}{k_{\mathrm{B}}T}\right)&\textrm{Creep}\\\
\frac{v_{\mathrm{T}}(H_{\mathrm{d}})}{x_{0}}\left(\frac{H_{z}-H_{\mathrm{d}}}{H_{\mathrm{d}}}\right)^{\beta}&\textrm{Depinning}\end{cases}$
(4)

where $v(H_{\mathrm{d}})$ is the creep prefactor corresponding to the velocity
at depinning, $T$ is the temperature, $k_{\mathrm{B}}$ is the Boltzmann
constant, $H_{\mathrm{d}}$ is the depinning field, $\Delta
E=k_{\mathrm{B}}T_{\mathrm{d}}[(H_{z}/H_{\mathrm{d}})^{-1/4}-1]$ is the
pinning potential, $v_{T}=v(H_{\mathrm{d}})(T_{\mathrm{d}}/T)^{\psi}$ is the
DW velocity in a disorder-free case, $x_{0}=0.65\pm 0.04$ is a dimensionless
constantPardo  _et al._ (2017), and $\psi=0.15$ and $\beta=0.25$ are depinning
exponentsJeudy  _et al._ (2018). Here we define the creep regime as $H\leq
H_{\mathrm{d}}$ and the depinning regime as $H\gtrsim
H_{\mathrm{d}}\left[1+(0.8(T_{\mathrm{d}}/T)^{-\psi})^{1/\beta}\right]$. The
fitting results, which yield a good agreement with the acquired experimental
data, are shown in Fig. 4(a) and Fig. 4(b). There is good agreement with both
the creep law and with the depinning law for $H>H_{\mathrm{d}}$.

The depinning field $H_{\mathrm{d}}$ that separates the creep and depinning
regimes of DW motion is marked on each curve in Fig. 4 by solid stars. The
only-non-universal (i.e. material and temperature dependent) parameters of the
fits are the depinning field $H_{\mathrm{d}}$, depinning velocity
$v_{\mathrm{T}}(H_{\mathrm{d}})$, and depinning temperature $T_{\mathrm{d}}$.
All these three parameters are increasing with the Au concentration in the
system (Fig. 4(c-e)). Following Ref. Jeudy  _et al._ , 2018, the material
dependent parameters $H_{\mathrm{d}}$ and $T_{\mathrm{d}}$ can be used to
estimate the variation with gold concentration of the coherence length
($\xi=\left[(k_{\mathrm{B}}T_{\mathrm{d}})^{2}/(2M_{\mathrm{s}}H_{\mathrm{d}}\sigma
t^{2})\right]^{1/3}$ ), and the pinning strength
($f_{\mathrm{pin}}\xi=\left[(k_{\mathrm{B}}T_{\mathrm{d}})^{3}/(\sigma\xi
t)\right]^{1/2}$). The results are reported in Table 1. It can be seen there
that the correlation length of the disorder is almost constant($\xi=23\pm 2$
nm), which suggests that the typical length scale of the variations of the
pinning potential (i.e. of the fluctuation of DW energy) is independent of the
gold concentration. Moreover, the increase of the pinning strength suggests
that fluctuations of DW energy are larger at the CoAu interface than at the
CoPt interface. This is compatible with a smearing of DW energy fluctuations
at the CoPt interface due to intermixing, which do not exist at the CoAu
interface since Co and Au do not form any alloyMassalski  _et al._ (1986).

Figure 4: OoP-only field-driven DW motion for Pt/Co/AuxPt1-x epilayers. (a)
Fits to the data by universal models for creep (dashed-dotted curve) and
depinning regimes (solid curve) on a linear-linear plot. The dotted lines are
linear fits to the viscous flow regime. (b) Semi-logarithmic plot of the same
velocity data as a function of $(\mu_{0}H)^{-1/4}$ to show fitting to the
creep regime more clearly. Data are shifted to the right respectively for
clarity of presentation. In each case, the solid stars show the depinning
field $H_{\mathrm{d}}$. (c)-(f) Changes of derived parameters with
concentration of gold $x$: depinning field $H_{\mathrm{d}}$, depinning
temperature $T_{\mathrm{d}}$, DW velocity in a film with pinning $v_{T}$, and
DW mobility $m$, respectively.

In each case the DW enters the flow regime, where $v\propto H_{z}$, at high
fields. As stated before, $v_{T}$ is assumed to be the velocity of DW in the
absence of pinning, hence DW mobility will be
$m=\mu_{0}^{-1}v_{T}/H_{\mathrm{d}}$. Using this calculated $m$ [Fig. 4(f)], a
line can be drawn with this slope, going through origin, shown in Fig. 4(a) as
a dotted line. Non-trivially, the line falls onto the experimental data-points
at high fields where $v\propto H_{z}$. This is interesting, as the calculated
DW mobility (that is a flow regime related parameter) is derived from creep
and depinning fits to the experimental data which are measured with magnetic
fields much lower than what viscous flow motion needs. In other words, we
could learn about the flow regime and its parameters without a need to reach
high, flow-related fields.

The standard one-dimensional (1D) model description of DW dynamics Thiaville
_et al._ (2012), which describes the DW motion in terms of its position $q$
and internal spin configuration $\phi$, predicts that the DW mobility $m$ in
the flow regime depends on Gilbert damping, $\alpha$, as follows:

$m=\begin{cases}\gamma\Delta/\alpha\qquad\qquad&\textrm{Steady flow}\\\
\gamma\Delta/\left(\alpha+\alpha^{-1}\right)&\textrm{Precessional
flow}\end{cases}$ (5)

where $\gamma$ is the gyromagnetic ratioMetaxas  _et al._ (2007). The Gilbert
damping for the three layers, $\alpha_{\mathrm{exp}}$, determined from these
dashed lines are given in Table 1. Eq. 5, which relates the DW mobility to
$\alpha_{\mathrm{exp}}$, only gives a solution in the case of precesssional
flow the trilayers capped with Pt and AuPt, whilst yielding a only steady flow
solution for the epilayer where $x=1$.

The Walker field, which separates the steady and precessional flow regimes,
can be estimated by using relation Thiaville  _et al._ (2012)

$H_{\mathrm{W}}=\alpha_{\mathrm{exp}}\sin\Phi_{\mathrm{W}}(H_{\mathrm{DMI}}-H_{\mathrm{demag}}\cos\Phi_{\mathrm{W}}),$
(6)

where $H_{\mathrm{demag}}$ is the DW demagnetizing field, and
$\cos\Phi_{\mathrm{W}}=([H_{\mathrm{DMI}}/H_{\mathrm{demag}}]-\sqrt{[H_{\mathrm{DMI}}/H_{\mathrm{demag}}]^{2}+8})/4.$
The values we obtain for $H_{\mathrm{W}}$ are also given in Table 1. The
Walker field increases markedly with $x$, due to the stronger DMI as the
inversion asymmetry is progressively broken. It can be seen that knowledge of
the DMI field and Gilbert damping are needed to determine the regime of the DW
dynamics with respect to the Walker field. Nevertheless, comparing these
values with $H_{\mathrm{d}}$, we can see that our results are self-consistent,
confirming which flow regime each sample enters on exceeding the depinning
field.

##  V Micromagnetic simulation

The Gilbert damping values $\alpha_{\mathrm{exp}}$ obtained from Eq. 5 is
unexpectedly large in the Pt/Co/Au sample Metaxas  _et al._ (2007), suggesting
that such a description, based on a 1D model, might not always be applicable
in these systems. This could be due to a strong effect of disorder, since the
1D model usually assumes DW motion in perfectly homogeneous conditions.

Figure 5: Field-induced DW motion in Pt/Co/AuxPt1-x epilayers. (a) and (b)
Visualisations of the micromagnetic simulations for Pt/Co/Au with
$\delta=8.5\%$, $\alpha=0.35$ at $t=0$ns. (c) Field-dependent velocity curves
along with the simulation from micromagnetic simulations and 1D model
calculations. It is shown that 1D model calculations cannot reproduce the data
while the micromagnetic simulations can follow change in velocity precisely.

In order to test this conjecture, we performed numerical micromagnetic
simulations using the mumax3 package Vansteenkiste  _et al._ (2014). We used
the experimentally measured magnetic parameters ($A$, $K_{\mathrm{u}}$,
$M_{\mathrm{Co}}$, and $D_{\mathrm{eff}}$) and implemented disorder by means
of Voronoi tessellation (see Fig. 5(a)) of the sort usually used to represent
a grain structure Kim and Yoo (2017); Zeissler  _et al._ (2017). Since our
trilayers are epitaxial they do not have a grain structure per se, but there
are magnetic layer thickness fluctuations, as revealed by the TEM cross-
sections. In each region the micromagnetic parameters were changed as expected
from a random thickness fluctuations that affect $K_{\mathrm{u}}$,
$D_{\mathrm{eff}}$ and $M_{\mathrm{Co}}$Hrabec  _et al._ (2017); Gross  _et
al._ (2018); Moretti  _et al._ (2017). The standard deviation of this
fluctuation is given by the parameter $\delta$. The simulated sample is
discretized in rectangular cells of $(2\times 2\times t)~{}\mathrm{nm}^{3}$,
where $t$ is the thickness of the sample. We consider a region of
$1.0~{}\upmu\mathrm{m}\times 0.5~{}\upmu\mathrm{m}$, with periodic boundary
condition along the $y$ direction. The DW is relaxed at the center as shown in
Fig. 5(b). The ‘grain’ size is set to 15 nm, which is a scale consistent with
our TEM analysis here, and close to the values of $\xi$ given in Table 1. We
simulate the DW dynamics for $50$ ns, in order to reach a stationary state in
which the calculated DW velocity is not influenced by the simulation time. We
performed different simulations combining different values of $\alpha$ and
$\delta$. As all the experimental measurements were done at room temperature,
the simulation temperature was set to $T=300$ K, accounted for by a thermal
field of the form given by Brown Vansteenkiste  _et al._ (2014); Brown (1963).
The best fits are shown in Fig. 5(c) for the three samples. The description of
the data with this approach is excellent. The corresponding best fitting
parameters ($\alpha_{\mathrm{\upmu M}}$ and $\delta$) are reported in Table 1.
The larger $\delta$ we had to introduce for the samples doped with Au reflects
the higher values of $f_{\mathrm{pin}}$ due to the more abrupt interface
between Co and Au, as previously commented on.

In every case, $\alpha_{\mathrm{\upmu M}}<\alpha_{\mathrm{exp}}$, although the
largest difference is for $x=1$. That damping parameter for Pt/Co/Au of is in
reasonable agreement with the values reported for typical Pt/Co systems
$\alpha\approx 0.3$Metaxas  _et al._ (2007). For Pt/Co/PtAu and Pt/Co/Pt are
slightly lower but just within an error bar of $\alpha_{\mathrm{exp}}$. We
calculated the prediction of the 1D model for $\alpha_{\mathrm{\upmu M}}$,
also shown Fig. 5(c) as solid lines, which do not coincide with the data at
all. The failure of the 1D model can be easily understood if we look at the
depinning fields obtained from the creep fits. When $H_{z}\approx
H_{\mathrm{d}}$ for $H_{\mathrm{d}}<H_{\mathrm{W}}$, as is the case for
Pt/Co/Au, disorder has a strong effect on the DW dynamics, which can be
reproduced successfully only with full micromagnetic simulations. On the other
hand, if $H_{z}\gg H_{\mathrm{d}}>H_{\mathrm{W}}$ disorder has a smaller
effect on the DW dynamics and the predictions of the 1D model are in agreement
with micromagnetic calculations in that field regime, as is the case for
Pt/Co/Pt and Pt/Co/PtAu. Not taking into account the DW’s internal dynamics
(such as the presence and collision of Bloch lines) can only lead to an upper
limit estimation for damping Moretti  _et al._ (2017), due to the fact that
all dissipation mechanisms are–implicitly–attributed to damping in that case
Yoshimura  _et al._ (2015).

The larger damping in the more strongly Au-doped samples could be due to a
larger spin pumping Tserkovnyak  _et al._ (2002b) at the Co/Au interface. An
increased spin Hall effect was indeed recently observed in Pt/Au/Co stacks
Mann and Beach (2017) and attributed to a larger spin transmission at the
Au/Co interface, which could also affect the spin pumping mechanism. However,
it is unlikely that it could explain alone the large difference with the
Pt/Co/Pt sample. Moreover, similar spin mixing conductances were observed in
Pt/Co and Au/Co by FMR Singh  _et al._ (2017). Non-local damping (dependent on
the spin configuration) Hankiewicz  _et al._ (2008); Zhang and Zhang (2009);
Weindler  _et al._ (2014) could also play an important role for narrow DWs and
could lead to an increased effective Gilbert damping 111The inclusion of the
non-local damping would imply a modification of the LLG equation Hankiewicz
_et al._ (2008); Zhang and Zhang (2009). Thus the value fitted with the
standard LLG must be seen as an effective damping, which averages out these
inhomogeneous contributions.. Hankiewicz et al., in particular, showed that
non-local damping increases with the disorder scattering rate and is strongly
enhanced by electron-electron interactions for weak spin polarization
Hankiewicz  _et al._ (2008), as could occur at the Co/Au interface or in case
of intermixing. Non-local effects will also be stronger for the larger
magnetization gradients around the narrower walls found for trilayers
containing more Au. Also, it was reported that for Au/Co interfaces the very
interfacial atomic layers are playing the most important role in damping of
the system, while for Pt/Co bilayers the effect can be extended to a few
atomic layers Barati  _et al._ (2014). So the observed higher thickness
fluctuations (i.e. interface roughness) and probable intermixing of layers for
Pt/Co/Au trilayers of this study can affect the effective damping.

##  VI Conclusions

We have prepared and studied a series of Pt/Co/Pt1-xAux epitaxial trilayers
with a well-defined crystallographic structure and controlled levels of
proximity magnetism and inversion symmetry breaking. We have determined their
static magnetic properties and the field-driven domain wall dynamics that they
support. Our main findings are as follows. Comparing the measured moments of
our samples allows us to estimate the degree of proximity magnetization in Pt,
which is about one-third of the total moment in a Pt/Co(0.6 nm)/Pt trilayer.
Adding Au increases the PMA, and also breaks the inversion symmetry and
increases the overall DMI. This DMI leads to homochiral left-handed Néel
domain walls, as predicted by theory for Pt/CoYang  _et al._ (2015),
suggesting that the DMI arising from the Au interface is weak, mirroring
findings for the spin-orbit torque Ryu  _et al._ (2014). Adding Au also
increases the depinning field and temperature for these walls, which is
consistent with a more abrupt magnetic interface in the absence of proximity
magnetism. The field-driven domain wall dynamics, including the crossover from
creep to viscous flow, can be described well by a micromagnetic model that
takes appropriate account of the disorder in the sample whereas a simple 1-D
model fails to do so. This micromagnetic modelling reveals lower values of
Gilbert damping than the 1-D case, due to the proper treatment of pinning,
although still rising up to $\alpha\approx 0.3$ as Au is added. The rise is
likely to be due to higher interface transparency leading to increased spin
pumping or to higher levels of non-local damping. The overall approach to the
study of these trilayers used here could be also extended to studying the
physics of skyrmionsFert  _et al._ (2013); Nagaosa and Tokura (2013), which
contain the same physics as the chiral walls although expanded into two
dimensions.

###### Acknowledgements.

This work was supported by the EU via WALL network (grant number
FP7-PEOPLE-2013-ITN 608031), the UK EPSRC (grant number EP/I011668/1), and the
Agence Nationale de la Recherche (France) under Contract No. ANR-14-CE26-0012
(Ultrasky). The work by E. M. was supported by Projects MAT2014-52477-C5-4-P
and MAT2017-87072- C4-1-P from the Spanish government, and Project No.
SA090U16 from the Junta de Castilla y Leon.

## References

  * Hoffmann and Bader (2015) A. Hoffmann and S. D. Bader, “Opportunities at the frontiers of spintronics,” Phys. Rev. Appl. 4, 047001 (2015).
  * de Jonge  _et al._ (1994) W. J. M. de Jonge, P. J. H. Bloemen, and F. J. A. den Broeder, “Experimental investigations of magnetic anisotropy,” (Springer, Berlin, 1994) Chap. 2.3.
  * Antel  _et al._ (1999) W. J. Antel, M. M. Schwickert, T. Lin, W. L. O’Brien, and G. R. Harp, “Induced ferromagnetism and anisotropy of Pt layers in Fe/Pt(001) multilayers,” Phys. Rev. B 60, 12933 (1999).
  * Fert (1991) A. Fert, “Magnetic and transport properties of metallic multilayers,” in  _Metallic Multilayers_ , Materials Science Forum, Vol. 59, edited by A. Chamberod and J. Hillairet (1991) pp. 439–480.
  * Chen  _et al._ (2013) G. Chen, J. Zhu, A. Quesada, J. Li, A. T. N’Diaye, Y. Huo, T. P. Ma, Y. Chen, H. Y. Kwon, C. Won, Z. Q. Qiu, A. K. Schmid Schmid, and Y. Z. Wu, “Novel chiral magnetic domain wall structure in Fe/Ni/Cu $(001)$ films,” Phys. Rev. Lett. 110, 177204 (2013).
  * Benitez  _et al._ (2015) M.-J. Benitez, A. Hrabec, A. P. Mihai, T. A. Moore, G. Burnell, D. McGrouther, C. H. Marrows, and S. McVitie, “Magnetic microscopy of topologically protected homochiral domain walls in an ultrathin perpendicularly magnetized Co film,” Nature Commun. 6, 8957 (2015).
  * Tetienne  _et al._ (2015) J.-P. Tetienne, T. Hingant, L. J. Martínez, S. Rohart, A. Thiaville, L. Herrera Diez, K. Garcia, J.-P. Adam, J.-V. Kim, J.-F. Roch, I. M. Miron, G. Gaudin, L. Vila, B. Ocker, D. Ravelosona, and V. Jacques, “The nature of domain walls in ultrathin ferromagnets revealed by scanning nanomagnetometry,” Nature Commun. 6, 6733 (2015).
  * Tserkovnyak  _et al._ (2002a) Y. Tserkovnyak, A. Brataas, and G. E. W. Bauer, “Enhanced Gilbert damping in thin ferromagnetic films,” Phys. Rev. Lett. 88, 117601 (2002a).
  * Hubert and Schäfer (1998) A. Hubert and R. Schäfer,  _Magnetic Domains - The Analysis of Magnetic Microstructures_ (Springer, Berlin, 1998).
  * Thiaville  _et al._ (2012) A. Thiaville, S. Rohart, É. Jué, V. Cros, and A. Fert, “Dynamics of Dzyaloshinskii domain walls in ultrathin magnetic films,” Europhys. Lett. 100, 57002 (2012).
  * Parkin and Yang (2015) S. S. P. Parkin and S.-H. Yang, “Memory on the racetrack,” Nature Nanotech. 10, 195 (2015).
  * Manchon (2014) A. Manchon, “Spin-orbitronics: A new moment for Berry,” Nature Phys. 10, 340 (2014).
  * Khvalkovskiy  _et al._ (2013) A. V. Khvalkovskiy, V. Cros, D. Apalkov, V. Nikitin, M. Krounbi, K. A. Zvezdin, A. Anane, J. Grollier, and A. Fert, “Matching domain-wall configuration and spin-orbit torques for efficient domain-wall motion,” Phys. Rev. B 87, 020402 (2013).
  * Emori  _et al._ (2013) S. Emori, U. Bauer, S.-M. Ahn, E. Martinez, and G.S.D. Beach, “Current-driven dynamics of chiral ferromagnetic domain walls,” Nature Mater. 12, 611–616 (2013).
  * Ryu  _et al._ (2014) K.-S. Ryu, S.-H. Yang, L. Thomas, and S. S. P. Parkin, “Chiral spin torque arising from proximity-induced magnetization,” Nature Commun. 5, 3910 (2014).
  * Hrabec  _et al._ (2014) A. Hrabec, N. A. Porter, A. Wells, M. J. Benitez, G. Burnell, S. McVitie, D. McGrouther, T. A. Moore, and C. H. Marrows, “Measuring and tailoring the Dzyaloshinskii-Moriya interaction in perpendicularly magnetized thin films,” Phys. Rev. B 90, 020402 (2014).
  * Wells  _et al._ (2017) A. W. J. Wells, P. M. Shepley, C. H. Marrows, and T. A. Moore, “Effect of interfacial intermixing on the Dzyaloshinskii-Moriya interaction in Pt/Co/Pt,” Phys. Rev. B 95, 054428 (2017).
  * Mihai  _et al._ (2013) A. P. Mihai, A. L. Whiteside, E. J. Canwell, C. H. Marrows, M.-J. Benitez, D. McGrouther, S. McVitie, S. McFadzean, and T. A. Moore, “Effect of substrate temperature on the magnetic properties of epitaxial sputter-grown Co/Pt,” Appl. Phys. Lett. 103, 262401 (2013).
  * Stoner (1938) E. C. Stoner, “Collective electron ferromagnetism,” Proc. R. Soc. London, Ser. A 165, 372 (1938).
  * Wilhelm  _et al._ (2004) F. Wilhelm, M. Angelakeris, N. Jaouen, P. Poulopoulos, E. Th. Papaioannou, Ch. Mueller, P. Fumagalli, A. Rogalev, and N. K. Flevaris, “Magnetic moment of Au at Au/Co interfaces: A direct experimental determination,” Phys. Rev. B 69, 220404 (2004).
  * Kashid  _et al._ (2014) V. Kashid, T. Schena, B. Zimmermann, Y. Mokrousov, S. Blügel, V. Shah, and H. G. Salunke, “Dzyaloshinskii-Moriya interaction and chiral magnetism in $3d$-$5d$ zigzag chains: Tight-binding model and ab initio calculations,” Phys. Rev. B 90, 054412 (2014).
  * Yang  _et al._ (2015) H. Yang, A. Thiaville, S. Rohart, A. Fert, and M. Chshiev, “Anatomy of Dzyaloshinskii-Moriya interaction at Co/Pt interfaces,” Phys. Rev. Lett. 115, 267210 (2015).
  * Ploog  _et al._ (2013) K. Ploog, W. Stolz, and L. Tapfer, “Effect of barrier configuration and interface quality on structural and electronic properties and MBE-grown AlxGa1-xAs/GaAs, AlxGa1-xSb/GaSb and AlxIn1-xAs/GaxIn1-xAs superlattices,” in  _Thin Film Growth Techniques for Low-Dimensional Structures_ , Nato Science Series B, Vol. 163, edited by R. F. C. Farrow, S. S. P. Parkin, P. J. Dobson, J. H. Neave, and A. S. Arrott (Springer Science & Business Media, 2013).
  * Tanner (2013) B. K. Tanner,  _Characterization of Crystal Growth Defects by X-Ray Methods_, Nato Science Series B: (Springer US, 2013).
  * Gabureac  _et al._ (2014) M. S. Gabureac, D. A. MacLaren, H. Courtois, and C. H. Marrows, “Long-ranged magnetic proximity effects in noble metal-doped cobalt probed with spin-dependent tunnelling,” New J. Phys. 16, 043008 (2014).
  * Björck and Andersson (2007) M. Björck and G. Andersson, “GenX: an extensible X-ray reflectivity refinement program utilizing differential evolution,” J. Appl. Crystallogr. 40, 1174 (2007).
  * Lin  _et al._ (1991) C.-J. Lin, G. L. Gorman, C. H. Lee, R. F. C. Farrow, E. E. Marinero, H. V. Do, H. Notarys, and C. J. Chien, “Magnetic and structural properties of Co/Pt multilayers,” J. Magn. Magn. Mater. 93, 194 (1991).
  * Bandiera  _et al._ (2012) S. Bandiera, R. C. Sousa, B. Rodmacq, and B. Dieny, “Enhancement of perpendicular magnetic anisotropy through reduction of Co-Pt interdiffusion in (Co/Pt) multilayers,” Appl. Phys. Lett. 100, 142410 (2012).
  * Lavrijsen  _et al._ (2015) R. Lavrijsen, D. M. F. Hartmann, A. van den Brink, Y. Yin, B. Barcones, R. A. Duine, M. A. Verheijen, H. J. M. Swagten, and B. Koopmans, “Asymmetric magnetic bubble expansion under in-plane field in Pt/Co/Pt: effect of interface engineering,” Phys. Rev. B 91, 104414 (2015).
  * Metaxas  _et al._ (2007) P. J. Metaxas, J. P. Jamet, A. Mougin, M. Cormier, J. Ferré, V. Baltz, B. Rodmacq, B. Dieny, and R. L. Stamps, “Creep and flow regimes of magnetic domain-wall motion in ultrathin Pt/Co/Pt films with perpendicular anisotropy,” Phys. Rev. Lett. 99, 217208 (2007).
  * Miron  _et al._ (2010) I. M. Miron, G. Gaudin, S. Auffret, B. Rodmacq, A. Schuhl, S. Pizzini, J. Vogel, and P. Gambardella, “Current-driven spin torque induced by the Rashba effect in a ferromagnetic metal layer,” Nature Mater. 9, 230 (2010).
  * Vaňatka  _et al._ (2015) M. Vaňatka, J.-C. Rojas-Sánchez, J. Vogel, M. Bonfim, M. Belmeguenai, Y. Roussigné, A. Stashkevich, A. Thiaville, and S. Pizzini, “Velocity asymmetry of Dzyaloshinskii domain walls in the creep and flow regimes,” J. Phys.: Condens. Matter 27, 326002 (2015).
  * Chikazumi (1997) S. Chikazumi,  _Physics of Ferromagnetism_ , 2nd ed. (Oxford University Press, Oxford, 1997).
  * Shepley  _et al._ (2018) P. M. Shepley, H. Tunnicliffe, K. Shahbazi, G. Burnell, and T. A. Moore, “Magnetic properties, domain-wall creep motion, and the Dzyaloshinskii-Moriya interaction in Pt/Co/Ir thin films,” Phys. Rev. B 97, 134417 (2018).
  * Je  _et al._ (2013) S.-G. Je, D.-H. Kim, S.-C. Yoo, B.-C. Min, K.-J. Lee, and S.-B. Choe, “Asymmetric magnetic domain-wall motion by the Dzyaloshinskii-Moriya interaction,” Phys. Rev. B 88, 214401 (2013).
  * Kabanov  _et al._ (2010) Y. P. Kabanov, Y. L. Iunin, V. I. Nikitenko, A. J. Shapiro, R. D. Shull, L. Y. Zhu, and C. L. Chien, “In-plane field effects on the dynamics of domain walls in ultrathin Co films with perpendicular anisotropy,” IEEE Trans. Magn. 46, 2220 (2010).
  * Kim  _et al._ (2018) D.-Y. Kim, M.-H. Park, Y.-K. Park, J.-S. Kim, Y.-S. Nam, D.-H. Kim, S.-G. Je, H.-C. Choi, B.-C. Min, and S.-B. Choe, “Chirality-induced antisymmetry in magnetic domain wall speed,” NPG Asia Mater. 10, e464 (2018).
  * Jeudy  _et al._ (2018) V. Jeudy, R. Dìaz Pardo, W. Savero Torres, S. Bustingorry, and A. Kolton, “Pinning of domain walls in thin ferromagnetic films,” Phys. Rev. B 98, 054406 (2018).
  * Pardo  _et al._ (2017) R. Diaz Pardo, W. Savero Torres, A. B. Kolton, S. Bustingorry, and V. Jeudy, “Universal depinning transition of domain walls in ultrathin ferromagnets,” Phys. Rev. B 95, 184434 (2017).
  * Massalski  _et al._ (1986) T. B. Massalski, H. Okamoto, P. R. Subramanian, L. Kacprzak, and W. W. Scott,  _Binary alloy phase diagrams_ , Vol. 1 (American Society for Metals, Metals Park, OH, 1986).
  * Vansteenkiste  _et al._ (2014) A. Vansteenkiste, J. Leliaert, M. Dvornik, M. Helsen, F. Garcia-Sanchez, and B. Van Waeyenberge, “The design and verification of MuMax3,” AIP Advances 4, 107133 (2014).
  * Kim and Yoo (2017) J.-V. Kim and M.-W. Yoo, “Current-driven skyrmion dynamics in disordered films,” Appl. Phys. Lett. 110, 132404 (2017).
  * Zeissler  _et al._ (2017) K. Zeissler, M. Mruczkiewicz, S. Finizio, J. Raabe, P. M. Shepley, A. V. Sadovnikov, S. A. Nikitov, K. Fallon, S. McFadzean, S. McVitie, T. A. Moore, G. Burnell, and C. H. Marrows, “Pinning and hysteresis in the field dependent diameter evolution of skyrmions in Pt/Co/Ir superlattice stacks,” Sci. Rep. 7, 15125 (2017).
  * Hrabec  _et al._ (2017) A. Hrabec, J. Sampaio, M. Belmeguenai, I. Gross, R. Weil, S. M. Chérif, A. Stashkevich, V. Jacques, A. Thiaville, and S. Rohart, “Current-induced skyrmion generation and dynamics in symmetric bilayers,” Nature Commun. 8 (2017).
  * Gross  _et al._ (2018) I. Gross, W. Akhtar, A. Hrabec, J. Sampaio, L. J. Martínez, S. Chouaieb, B. J. Shields, P. Maletinsky, A. Thiaville, S. Rohart, and V. Jacques, “Skyrmion morphology in ultrathin magnetic films,” Phys. Rev. Mater. 2, 024406 (2018).
  * Moretti  _et al._ (2017) S. Moretti, M. Voto, and E. Martinez, “Dynamical depinning of chiral domain walls,” Phys. Rev. B 96, 054433 (2017).
  * Brown (1963) W. F. Brown, “Thermal fluctuations of a single-domain particle,” J. Appl. Phys. 34, 1319 (1963).
  * Yoshimura  _et al._ (2015) Y. Yoshimura, K.-J. Kim, T. Taniguchi, T. Tono, K. Ueda, R. Hiramatsu, T. Moriyama, K. Yamada, Y. Nakatani, and T. Ono, “Soliton-like magnetic domain wall motion induced by the interfacial Dzyaloshinskii-Moriya interaction,” Nature Phys. 12, 157 (2015).
  * Tserkovnyak  _et al._ (2002b) Y. Tserkovnyak, A. Brataas, and G. E. W. Bauer, “Spin pumping and magnetization dynamics in metallic multilayers,” Phys. Rev. B 66, 224403 (2002b).
  * Mann and Beach (2017) M. Mann and G. S. D. Beach, “Reduction of in-plane field required for spin-orbit torque magnetization reversal by insertion of Au spacer in Pt/Au/Co/Ni/Co/Ta,” APL Materials 5, 106104 (2017).
  * Singh  _et al._ (2017) B. B. Singh, S. K. Jena, and S. Bedanta, “Study of spin pumping in co thin film vis-à-vis seed and capping layers using ferromagnetic resonance spectroscopy,” J. Phys. D: Appl. Phys. 50, 345001 (2017).
  * Hankiewicz  _et al._ (2008) E. M. Hankiewicz, G. Vignale, and Y. Tserkovnyak, “Inhomogeneous Gilbert damping from impurities and electron-electron interactions,” Phys. Rev. B 78, 020404 (2008).
  * Zhang and Zhang (2009) S. Zhang and S. S.-L. Zhang, “Generalization of the Landau-Lifshitz-Gilbert equation for conducting ferromagnets,” Phys. Rev. Lett. 102, 086601 (2009).
  * Weindler  _et al._ (2014) T. Weindler, H. G. Bauer, R. Islinger, B. Boehm, J.-Y. Chauleau, and C. H. Back, “Magnetic damping: Domain wall dynamics versus local ferromagnetic resonance,” Phys. Rev. Lett. 113, 237204 (2014).
  * Note (1) The inclusion of the non-local damping would imply a modification of the LLG equation Hankiewicz  _et al._ (2008); Zhang and Zhang (2009). Thus the value fitted with the standard LLG must be seen as an effective damping, which averages out these inhomogeneous contributions.
  * Barati  _et al._ (2014) E. Barati, M. Cinal, D. M. Edwards, and A. Umerski, “Gilbert damping in magnetic layered systems,” Phys. Rev. B 90, 014420 (2014).
  * Fert  _et al._ (2013) A. Fert, V. Cros, and J. Sampaio, “Skyrmions on the track,” Nature Nanotech. 8, 152 (2013).
  * Nagaosa and Tokura (2013) N. Nagaosa and Y. Tokura, “Topological properties and dynamics of magnetic skyrmions,” Nature Nanotech. 8, 899 (2013).


# 1D Anisotropic Surface Wave Tomography with Bayesian Inference

John Keith. Magali  Université de Lyon, UCBL, CNRS, LGL-TPE, 69622
Villeurbanne, France
Email<EMAIL_ADDRESS>
###### Abstract

Classically, anisotropic surface wave tomography is treated as an optimisation
problem where it proceeds through a linearised two-step approach. It involves
the construction of 2D group or phase velocity maps for each considered
period, followed by the inversion of local dispersion curves inferred from
these maps for 1D depth-functions of the elastic parameters. Here, we cast the
second step into a fully Bayesian probability framework. Solutions to the
inverse problem are thus an ensemble of model parameters (i.e. 1D elastic
structures) distributed according to a posterior probability density function
and their corresponding uncertainty limits. The method is applied to
azimuthally-varying synthetic surface wave dispersion curves generated by a
3D-deforming upper mantle. We show that such a procedure captures essential
features of the upper mantle structure. The robustness of these features
however strongly depend on the wavelength of the wavefield considered and the
choice of the model parameterisation. Additional information should therefore
be incorporated to regularise the problem such as the imposition of
petrological constraints to match the geodynamic predictions.

##  1 Probabilistic approaches to surface wave tomography: An introduction

Conventional surface wave tomography is usually implemented using a two-step
approach (e.g Nakanishi and Anderson, 1983; Nataf et al., 1984; Trampert and
Woodhouse, 1995; Romanowicz, 2002; Ritzwoller et al., 2002). The first step
involves the inversion of the arrival times of each period considered in the
measured source-receiver dispersion data to infer 2D group or phase velocity
maps at a given period. Using the 2D velocity maps, the second step proceeds
by inverting a dispersion curve at a given geographical location to estimate
the 1D velocity structure beneath this location. One may then build a smooth
3D velocity model by the juxtaposition of the inferred 1D models followed by
interpolating them.

The tomography problem is often solved by applying first-order corrections of
the forward function $g$ around a reference model $\textbf{m}_{0}$.
Mathematically. this translates to:

$\mathbf{d}=g(\mathbf{m_{0}})+\frac{\partial g}{\partial m}\Delta\mathbf{m}.$
(1)

Doing so allows it to be treated as a linearised inverse problem. For
instance, one may estimate the 2D phase velocity maps by minimizing an
objective function containing a data residual term and more than one
regularization terms such as:

$S=||\mathbf{Gm}-\mathbf{d}||+R_{1}||\mathbf{m}||+R_{2}||\mathbf{Dm}||$ (2)

where $G=\frac{\partial g}{\partial m}$ is now a mathematical forward operator
which now refers to the linearised physics between the model parameters (in
this case surface wave velocities), and the data (dispersion data), $D$ is a
second-order smoothing operator, and $R_{1}$ and $R_{2}$ are the damping and
smoothing parameters, respectively. The last two parameters will often dictate
the trade-off between the data misfit (i.e., how well the model parameters
predict the data), the proximity of the estimated model from its reference
state, and the degree of smoothing in the inverted model. Moreover, the
regularization parameters are chosen ad hoc. Hence, the inverted model may be
susceptible to non-data driven constraints and that these constraints shroud
essential information provided by the data. Lastly, optimization techniques
such as this lack the capability to estimate model uncertainties.

The problems associated with non-uniqueness and quantification of
uncertainties, coupled by the ever-growing computational capacity of modern
supercomputers led to the advancement of probabilistic approaches to
geophysical inverse problems, as first exemplified by Mosegaard and Tarantola
(1995). Following this pioneering study, a volume of studies that involve
probabilistic approaches to geophysical inverse problems have been published
in seismology (e.g. Lomax et al., 2000; Shapiro and Ritzwoller, 2002; Husen et
al., 2003; Bodin and Sambridge, 2009; Debski, 2010; Bodin et al., 2016), and
rapidly growing in the field of geodynamics (e.g. Baumann et al., 2014;
Baumann and Kaus, 2015; Morishige and Kuwatani, 2020; Ortega-Gelabert et al.,
2020).

Casting the inverse problem in a probabilistic framework allows one to utilise
the original non-linear mapping $g$ between the data and the model within its
forward procedure. However due to the use of sampling-based methods, one is
then forced to solve the forward model numerous times depending on the number
of candidate models to be sampled. Fortunately, various techniques are
available to address such complications that are not solely based on
heuristics. Here, we restrict ourselves with Bayesian inference, that is, a
form of statistical inference that formulates our solution as an a posteriori
probability initially based on the information we have prior to evaluating the
inverse problem. The goal is therefore not to create an ensemble of solutions
that follow a certain probability distribution ex nihilo, but to update a
prior probability based on valuable information provided by the data.

In such schemes, the parameter space has to be explored for best possible
model candidates that could match the predictions observed at the surface.
Grid-search algorithms however are time consuming and thus uniform sampling
may not be performed efficiently. As such, direct-search algorithms have been
introduced in geophysical inverse problems that sample candidate models within
a subset instead of the entire parameter space. A specific class of ergodic
algorithms is called Markov chain Monte Carlo (McMC) methods, which has been
initially used to solve problems in physics (Metropolis and Ulam, 1949;
Metropolis et al., 1953), but is now widely used in geophysical inverse
problems. In this algorithm, the parameter space is randomly sampled based on
our current state of knowledge of a given model candidate. The sequence of
searching new model candidates depend on an acceptance probability that is
often determined by satisfying a detailed balance condition to ensure
stationarity of the desired solution, the most common being the Metropolis-
Hastings algorithm (Hastings, 1970).

In this manuscript, we cast the second step of the surface wave tomography
problem in a Bayesian framework. That is, we assume that we have completely
inferred 2D phase velocity maps of the regions considered, and invert for 1D
velocity structures that explain the dispersion curves built from these maps.
The problem is applied to synthetic data for each geographical location
considered whose anisotropic signatures are solely influenced by convective
flow in the mantle beneath it. The solution is a marginal posterior
distribution of 1D velocity models that best explain the data accompanied by
their uncertainty limits. We show that even with strong a priori constraints,
conventional surface tomography falls short to capture the complete picture
related to mantle deformation. Still, some of the notable features are
resolved more or less. This work builds upon the hypothesis that adding
geodynamical and petrological constraints would allows us to reduce the number
of acceptable tomographic models that are consistent with the geodynamical
predictions.

##  2 Methods and experiments design

In this section, we discuss the full implementation of a 1D anisotropic
surface wave tomography in a full Bayesian parameter search approach. We
highlight in full detail the (1) model parameterisation, (2) the forward
problem, (3) the data, and (4) the inverse approach. The method is applied to
synthetic surface wave dispersion curves produced by simple setups of an
intrinsically anisotropic upper mantle.

###  2.1 Model parameterisation of a 1D Earth structure

####  2.1.1 Radial anisotropy component

Surface waves are sensitive to 13 depth parameters which are just a linear
combination of the full elastic tensor $S_{ij}$ (Montagner and Nataf, 1986).
In particular, Love and Rayleigh phase velocities are sensitive to five depth
depth parameters which make up an azimuthally-averaged vertically transverse
isotropic (VTI) medium. These parameters are also known as the Love
parameters, and by convention, are designated as $A$, $C$, $F$, $L$, and $N$.
Note that these functions are constrained by the isotropic phase velocities,
and are independent of the azimuth of surface wave anisotropy. As a
supplementary, the seismic wave velocities propagating either parallel or
perpendicular to the symmetry axis can be written as:

$\displaystyle V_{PH}=\sqrt{\frac{A}{\rho}}$ (3) $\displaystyle
V_{PV}=\sqrt{\frac{C}{\rho}}$ (4) $\displaystyle V_{SH}=\sqrt{\frac{N}{\rho}}$
(5) $\displaystyle V_{SV}=\sqrt{\frac{L}{\rho}},$ (6)

where $\rho$ is the density of the medium. Here, it is important to understand
that the vertical symmetry axis need not be the fast axis of anisotropy.
Hence, the relative magnitude between $N$ and $L$, and $C$ and $A$ are
interchangeable. Most anisotropic tomography studies interpret $L>N$ as
vertical flow, since to first-order, the direction of shear is presumed to be
vertical and thus aligned with vertically propagating $S-$waves (e.g.
Montagner, 1994). Likewise, $L<N$ is often interpreted in terms of horizontal
flow. When the flow has both horizontal and vertical components, then the
resulting anisotropy will be ambiguous. This requires resolving the tilt of
anisotropy (Montagner and Nataf, 1988; Montagner and Jobert, 1988).

We constrain radial anisotropy by using a more compact form related to the
Love parameters. For $P-$waves, the strength of radial anisotropy can be
expressed as $\phi=C/A$, whereas for $S-$waves, it is given by $\xi=N/L$.
Finally, there exists another anisotropic parameter which relates to the
ellipticity $\eta=F/A-2L$. The phenomenology of $\eta$ can be understood
through the parameter $F$ which controls the velocity along the direction
between the fast and slow velocities.

####  2.1.2 Azimuthal anisotropy component

Assuming a slightly anisotropic medium, azimuthal anisotropy in surface waves
can be decomposed into two terms that depend on its azimuth $\theta$. These
two terms are usually called the 2$\theta$ and 4$\theta$ components, and are
small perturbations around the isotropic phase velocities. The 2$\theta$ and
4$\theta$ components are sensitive to eight depth functions 2$\theta$:
$G_{s},G_{c},B_{s},B_{c},H_{s},H_{c}$, and 4$\theta$: $C_{s}$, $C_{c}$
(Montagner and Nataf, 1986).

We will only work with azimuthal anisotropy in Rayleigh waves. Rayleigh waves
are much more sensitive to the 2$\theta$ terms than the 4$\theta$ terms
(Maupin and Park, 2015). Thus, to first-order, we could eliminate the
parameters $C_{s}$ and $C_{c}$ from the inversions. Surface waves also poorly
resolve the parameters $H_{s}$ and $H_{c}$ (Bodin et al., 2016). Thus, we
could reduce the model dimensionality associated with azimuthal anisotropy
down to four parameters $G_{s},G_{c},B_{s}$, and $B_{c}$.

####  2.1.3 Pseudo-regularisation

By using compact notations and first-order approximations, we are able to
reduce the possible number of parameters to be inverted for. Four of which:
$A$, $L$, $\xi$, $\phi$, and $\eta$ are sensitive to Love and Rayleigh phase
velocities, and the remaining four: $G_{s},G_{c},B_{s}$, and $B_{c}$ are
sensitive to the azimuthal variations in Rayleigh phase velocities. Working
with synthetic data allows us to access the correct values of the parameters
defining our Earth model. As we wish to compare conventional tomography with
geodynamic tomography, here, we will be placing ourselves in the best case
scenario. We impose strong a priori constraints to our solution in the
tomographic problem. These constraints can be regarded as regularisation
parameters which limit the regions in the parameter space to search through.
Here, we assume that we have the correct values relating to the $P-$wave
structure $A$, $\phi$, $\eta$, $B_{s}$, and $B_{c}$; and thus, only invert for
$S-$wave-related structures $L$, $\xi$, $G_{c}$, and $G_{s}$. The list of
parameters to/what not to invert for is summarized in Table 1.

Table 1: Depth functions constrained by surface waves and their azimuthal variations. Parameter | Fixed | Inverted for
---|---|---
$A$ | $\checkmark$ |
$L$ |  | $\checkmark$
$\xi$ |  | $\checkmark$
$\phi$ | $\checkmark$ |
$\eta$ | $\checkmark$ |
$G_{c}$ |  | $\checkmark$
$G_{s}$ |  | $\checkmark$
$B_{c}$ | $\checkmark$ |
$B_{s}$ | $\checkmark$ |

The 1D Earth structure spans from the surface down to a depth of 400 km. It is
subdivided into 62 layers of equal thicknesses where each layer has a given
value of the nine depth functions (see Table 1). Thus, the total number of
unknowns to be inverted for is 62 $\times$ 4 = 248. Bayesian inversion is
therefore a suitable method to treat such an under-determined inverse problem.
To reduce the model dimensions further, the mantle is parameterised using a
piecewise cubic Hermite polynomials at fixed control points. In the
inversions, we vary the unknown parameters $L$, $\xi$, $G_{c}$, and $G_{s}$ at
the control points and in-between these points, the parameters are
interpolated. Below the 400 km, we impose isotropic PREM (Dziewonski and
Anderson, 1981). The 1D structure is illustrated in Fig. 1

Figure 1: Schematic diagram of the 1D parameterisation. The entire region is
parameterised with a piecewise cubic Hermite polynomials based on a number of
control points. Here, the control points are fixed in depth. The model
parameters are then varied at these points using an McMC sampling algorithm.
The layers in between are interpolated, and anything below the 400 km is
isotropic PREM.

The model vector is thus defined as:

$\mathbf{m}=[\mathbf{L},\bm{\xi},\mathbf{G_{s}},\mathbf{G_{c}}],$ (7)

where bold faces indicate that each parameter is also a vector of size 62.

###  2.2 The forward problem

For each step in the McMC algorithm, the forward problem is evaluated using
the proposed model (see Fig. 1) as input. The predicted data from this model
is then compared with the observed synthetics.

Isotropic Rayleigh $c_{R}(T)$ and Love $c_{L}(T)$ dispersion curves are
computed using normal mode summation in a spherical Earth (Smith and Dahlen,
1973). Here, the computations are carried out in a fully non-linear approach
following the method developed by Saito (1967, 1988). The software package
`DISPER80` (Saito, 1988) takes 1D depth profiles of $V_{p}$, $V_{s}$, $\rho$,
$\xi$, $\phi$, and $\eta$ as inputs to compute for $c_{R}(T)$ and $c_{L}(T)$
and their associated sensitivity kernels using a Runge-Kutta matrix
integration scheme.

Following the pioneering work of Montagner and Nataf (1986), the azimuthal
variations in surface wave phase velocities $c_{1}$ and $c_{2}$ can be
evaluated using the following expressions:

$\displaystyle c_{1}(T)=\int_{z=0}^{\infty}\left(B_{c}(z)\frac{\partial
c_{r}(T)}{\partial A}+G_{c}(z)\frac{\partial c_{r}(T)}{\partial L}\right)dz$
(8) $\displaystyle c_{2}(T)=\int_{z=0}^{\infty}\left(B_{s}(z)\frac{\partial
c_{r}(T)}{\partial A}+G_{s}(z)\frac{\partial c_{r}(T)}{\partial L}\right)dz.$
(9)

Eqs. (8) and (9) imply that the azimuthal variations in Rayleigh waves are
linearized around the reference VTI model after averaging azimuthally. Such
approximations are valid assuming the medium is quasi-anisotropic (Montagner
and Nataf, 1986; Maupin and Park, 2015).

###  2.3 The data

For each geographical location, we can express the local dispersion curve as
the sum of the isotropic dispersion curves and their azimuthal variations
giving:

$c(T,\theta)=c_{0}(T)+c_{1}(T)\cos(2\theta)+c_{2}(T)\sin(2\theta),$ (10)

where $T$ is the period, and $\theta$ is the azimuth of the propagating
surface wave.

For Rayleigh waves, we invert $c_{0}$, $c_{1}$, and $c_{2}$ whereas for Love
waves, we only invert $c_{0}$. For simplicity, we neglect the higher-order
terms associated with the elastic parameter $N$. such assumptions are valid
due to sparse sampling, low sensitivity, or high noise levels.

###  2.4 Quantification of anisotropy

Anisotropic surface wave tomography is capable of constraining azimuthal and
radial anisotropy. The level of radial anisotropy can be quantified through
the parameter $\xi$. Conversely, there are a variety of ways to quantify the
strength of azimuthal anisotropy, and its fast azimuth $\Psi$. Here, we
quantify it in terms of the peak-to-peak anisotropy:

$\mathrm{azi}=\frac{2G}{L},$ (11)

where $G$ = $\sqrt{G_{c}^{2}+G_{s}^{2}}$. The azimuth of fast propagation is
given by:

$\Psi=0.5\arctan\left(\frac{G_{s}}{G_{c}}\right).$ (12)

###  2.5 1D surface wave tomography with Bayesian inference

The inverse problem is cast in a full Bayesian procedure where the solution is
an ensemble of models distributed according to the posterior pdf $p$(m—dobs),
accompanied by their uncertainty bounds. In this framework, Bayes’ theorem
holds:

$p(\mathbf{m}|\mathbf{d_{obs}})\propto
p(\mathbf{m})p(\mathbf{d_{obs}}|\mathbf{m}).$ (13)

The parameter space is searched using a Markov chain Monte Carlo (McMC)
algorithm. To produce reasonable acceptance rates, we employed the adaptive
perturbation scheme.

####  2.5.1 The likelihood function

The likelihood function $p$(dobs—m) quantifies how well the model parameters
fit the observed data. In the context of our problem, it is loosely based on
the $L^{2}$-norm cost function in that it measures the level of misfit between
the predictions and the observations. Here, it is essential to distinguish the
data residuals related from errors in measurement $\bm{\epsilon_{d}}$, and
from the modeling error due to the use of an incorrect forward model
$\bm{\epsilon_{g}}$. Assuming that the errors are independent and describe a
random process, the forward problem can be written as:

$\mathbf{d_{obs}}=g(\mathbf{m})+\bm{\epsilon_{d}}+\bm{\epsilon_{g}}.$ (14)

Thus, $\mathbf{d_{obs}}$ can also be seen as a random process, and the
likelihood distribution can be formulated in terms of the pdf of the data
errors:

$p(\mathbf{d_{obs}|{m}})=p(\bm{\epsilon_{d}}+\bm{\epsilon_{g}}).$ (15)

If we then assume that the errors are uncorrelated and follow a univariate
Gaussian distribution with zero mean, and variance $\sigma^{2}$ where
$\sigma^{2}=\sigma_{d}^{2}+\sigma_{g}^{2}$, we can write the likelihood
function as an exponential giving:

$p(\mathbf{d_{obs}|{m}})=\frac{1}{(2\pi\sigma^{2})^{N/2}}\exp\bigg{[}\frac{-||\mathbf{d_{obs}}-g(\mathbf{m})||^{2}}{2\sigma^{2}}\bigg{]}.$
(16)

where $N$ is the size of the data vector. Since the goal of any optimization
problem is to minimise the $L^{2}$ cost function, minimising it tantamounts to
maximising the probability of the Gaussian likelihood function given by eq.
(16).

The above formulation can be utilised to construct the likelihood function for
surface wave data. For instance, the likelihood function corresponding to a
single dispersion measurement can be written as:

$p(\mathbf{c_{obs}|m})=\frac{1}{(2\pi\sigma_{c}^{2})^{N/2}}\exp\bigg{[}\frac{-||\mathbf{c_{obs}}-c||^{2}}{2\sigma_{c}^{2}}\bigg{]},$
(17)

where $\mathbf{m}$ is the 1D velocity model described in Fig. 1, $N$ is the
number of discreet periods, and $\sigma_{c}^{2}$ is the estimated variance of
the data noise. The likelihood functions of the 2$\theta$ terms can be cast in
the same manner.

####  2.5.2 The prior distribution

One of the flexibilities of the Bayesian framework is that it enables one to
account for prior information provided that it can be formulated as a
probability distribution. The choice of the prior depends on the type of
geophysical process we aim to tackle. In the context of seismic tomography,
the prior information depends on a range velocity structures that are
reasonable for the Earth in general. In practice, the prior information is
constrained by existing studies.

At this point forward, we assume minimal prior knowledge and hence, make use
of uniform prior distributions with wide bounds. Although we acknowledge that
using uniform distributions may be a naive way to setup the prior, working
with such a simple distribution would already suffice when demonstrating
proofs of concept. Here, we know the exact values of the model parameters, and
by imposing wide uniform priors, we are able to assess the efficiency of the
method by placing ourselves in the worst case scenario. Indeed, a subject of
future work is to consider other forms of the prior distribution, for example
non-informative priors or even hiearachical Bayes (Malinverno and Briggs,
2004).

Let us now consider a given model parameter $m_{i}$. The prior $p(m_{i})$ is
prescribed a constant value over a given range of values defined by
$[m_{\mathrm{min}},m_{\mathrm{max}}]$. The prior distribution is thus given
by:

$p(m_{i})=\begin{cases}0&m_{i}>m_{\mathrm{max}},m_{i}<m_{\mathrm{min}}\\\
\frac{1}{\Delta m}&m_{\mathrm{min}}\leq m_{i}\leq
m_{\mathrm{max}}.\end{cases}$ (18)

Eq. (18) is interpreted as follows. Suppose that we draw a sample for the
specific model parameter ${}_{m}i$ from the proposal distribution $q$. If
$m_{i}$ is out of bounds, then the proposal is automatically rejected because
the value is not specified by the prior. If $m_{i}$ is within the prior
bounds, then the proposal is accepted with condition based on the the
acceptance probability $A$. Choosing narrow bounds therefore imposes hard
constraints to the model parameters giving less emphasis to the information
provided by the data.

Assuming prior independence, the prior $p(\mathbf{m})$ can be written as a
product of 1D priors on each unknown parameter considered giving us:

$p(\mathbf{m})=\prod_{i=1}^{\mathrm{Nlyrs}}\bigg{[}p(L_{i})p(\xi_{i})p(Gs_{i})p(Gc_{i})\bigg{]},$
(19)

where $\mathrm{Nlyrs}$ is the number of layers of the 1D Earth model. Eq. (19)
implies that the probability of accepting the transition is automatically zero
should one of the parameters be outside their respective bounds.

####  2.5.3 The Sampling algorithm

We use a Markov chain Monte Carlo (McMC) algorithm to search the parameter
space for 1D Earth model candidates that could explain the surface wave
dispersion measurements. The sampler initiates by randomly drawing a reference
model from the prior followed by the evaluation of the likelihood function.
Within the Markov chain, the current 1D Earth model is perturbed at their
control points (see the red arrows in Fig. 1 to transition into a new state.
This is performed by randomly selecting one of the following set of moves:

  1. 1.

Change the Love parameter $L$ values of all control points according to a
Gaussian distribution centered at the current value of $L$.

  2. 2.

Change $\xi$ values of all control points according to a Gaussian distribution
centered at the current value of $\xi$.

  3. 3.

Change $G_{s}$ values of all control points according to a Gaussian
distribution centered at the current value of $G_{s}$.

  4. 4.

Change $G_{c}$ values of all control points according to a Gaussian
distribution centered at the current value of $G_{c}$.

If the proposed 1D Earth model is within their respective prior bounds, we
then solve the forward problem completely. The computed dispersion curves are
then compared with the observed synthetics following the evaluation of the
likelihood function. The resulting probability is then used to evaluate the
acceptance probability via the Metropolis-Hastings algorithm. The outcome of
the algorithm determines whether the proposed model is added to the posterior
distribution. Should it be rejected, the current model is counted successively
in the next iteration.

####  2.5.4 An adaptive perturbation scheme

We opted for a more dynamic perturbation, that is, the perturbations vary
depending on a given situation. Here, we employ the McMC sampler with an
adaptive perturbation scheme based on the acceptance rate. This requires
keeping track of the acceptance rate for a given model parameter on the fly.
Let us now denote $N_{\mathrm{pop}}$ to be the population size, that is, total
number of samples in the inversion, and $M$ to be the total number of accepted
models within $N_{\mathrm{pop}}$. The acceptance rate corresponding to the
entire population is just:

$\mathrm{acceptance}=\frac{M}{N_{\mathrm{pop}}}\times 100\%$ (20)

Next, we need a sizeable amount of samples $n$ within $N_{\mathrm{pop}}$ to
allow for the relaxation of the acceptance rates, that is, the period at which
the acceptance rates are in stable conditions. The population
$N_{\mathrm{pop}}$ are then separated into different cycles (i.e., sample
window) that are multiples of $n$. For instance, choosing $n=N_{\mathrm{pop}}$
just pertains to the acceptance rates of the entire population, whereas
choosing $n=1$ would constantly reset the counters for the proposal and the
accepted models at every iteration. At the end of every cycle, the current
values of the acceptance rates are then used to determine whether to increase
or decrease the perturbation. If the acceptance is less than 20$\%$, then the
perturbations will be reduced by 25$\%$ of its current value. Likewise, if the
acceptance is more than 50$\%$, then the perturbations will be increased by
25$\%$ instead. The counters for the proposal and the acceptance are then
reset, and the new acceptance rate is recalculated in the subsequent cycle. As
an example, suppose that we have $N_{\mathrm{pop}}=50000$ samples, and a
subset of $n=5000$, we thus having ten cycles. After the first 5000 iterations
$l$ (first cycle), the acceptance rate at $l=5001$ will determine whether to
increase or decrease the perturbations by 25$\%$. Resetting of the counters
will then be ensued regardless. The appraisal of the acceptance rates are then
marked after every $n$ intervals,that is, at $l=10000,15000,20000$ and so on.

##  3 Synthetic example: Application to a 3D deforming upper mantle

We perform 1D anisotropic surface wave tomography at 32 different geographical
locations. To efficiently implement the inversions for each location, the
algorithm is parallellised in such a way that multiple Markov chains can
simultaneously search the parameter space independently from one another.
Here, one Markov chain is allocated to one processor. At the end of the
inversion procedure, the ensemble of models from each chain are then gathered
to construct the posterior probability distribution. Since we consider 20
independent chains for each geographical location, and we have 32 locations in
total, we need a hefty 32 $\times$ 20 = 640 cores to implement the parallel
scheme.

###  3.1 Upper mantle flow induced by a sinking anomaly

We first test the method to a deforming viscous upper mantle induced by a
negatively buoyant spherical anomaly, akin to a Stoke’s sinker. As illustrated
in Fig. 2, the surrounding material responds to the sinking anomaly by
producing a return flow. The return flow, together with the downward motion of
the anomaly, generates local convection cells whose scales are consistent with
that of the upper mantle (which we set at $L_{s}$ = 400 km). We setup the
geographical locations of the local surface wave dispersion curves in such a
way that they provide a good coverage of the anomaly. Here, we use 32
locations that slice the region evenly into two (refer to the dashed lines).
The observed synthetics are generated entirely by a 3D deforming upper mantle.
We refer the reader to Magali et al. (2020) for further details.

Figure 2: Snapshot of a 3D deforming mantle induced by a sinking spherical
anomaly. The model domain is of the size 400 km $\times$ 400 km $\times$ 400
km. The isovolumetric gradients correspond to the 3D temperature profile of
the region, and the superimposed vector field is the flow induced by the
spherical body. The dashed lines at the surface represent the 32 geographical
locations of the local surface wave dispersion curves.

(a) Figure 3: Reconstructed 2D surface wave maps with added noise. (A)
Rayleigh. (B) Love.

(b)

To mimic real-Earth observations, we tarnish the observed synthetics with
noise. We added random uncorrelated noise with standard deviation
$\sigma_{c_{R,L}}$ = 0.001 km/s for the isotropic dispersion curves
$\mathbf{c_{R}}$ and $\mathbf{c_{L}}$, and $\sigma_{c_{1,2}}$ = 0.0005 km/s
for the azimuthal components $\mathbf{c_{1}}$ and $\mathbf{c_{2}}$. Fig. 4
shows the resulting dispersion curves at one given location. Solid blue lines
are the synthetic dispersion curves without noise and the red dots represent
the synthetics with added random uncorrelated noise. By stacking the
dispersion curves laterally (i.e., placing them side-by-side), we are able to
visualize the true structure in the data space. Fig. 3(a) shows the effect of
the density anomaly onto the resulting surface wave dispersion maps. The
negative subsidence observed in the Rayleigh map corresponds to an increases
in its speed as it traverses the cold anomaly.

(a)

(b)

(c)

(d)

Figure 4: Surface wave phase velocity dispersion curves and their azimuthal
variations at a specific geographical location. Solid blue lines are the
correct values and the scatter plots are the ones added with noise and are to
be inverted.

For each geographical location, the inversion consists of 20 independent
Markov chains each containing 1.0 $\times$ $10^{6}$ samples initiated at a
random 1D Earth model drawn from the prior distribution. The prior bounds of
each model parameter is summarized in Table 2.

Table 2: Prior ranges of the unknown model parameters. |  $L$(GPa) | $\xi$ |  $G_{c}$(GPa) |  $G_{s}$(GPa)
---|---|---|---|---
min | 20 | 0.7 | -5 | -5
max | 150 | 1.3 | 5 | 5

The models are then collected after a burn-in period of 900,000 samples to
ensure the chains are in steady-state and are properly sampling the posterior
distribution. Note that even though the models are randomly initiated, the
positions of the control points are fixed, and thus may have strong
implications on the recovered structures.

Fig. 5 shows the mean velocity structure and the mean radial anisotropy
recovered from Bayesian inversion, and the correct structures. The 2D models
are constructed by placing the recovered 1D mean structures side-by-side.
Surface waves were able to successfully map the most prominent feature, that
is, the seismic anomaly associated with the denser sphere. Radial anisotropy
was also successfully recovered. However, one of the major drawbacks of
surface waves is that they are more sensitive at shallower depths. Hence, the
top layers are better resolved than the bottom layers. We thus expect the
upper portion to exhibit less model uncertainties that the bottom half.
Additionally, some essential features are smeared vertically, which we
attribute to the inherent long period nature of surface waves. The structure
also appears to be smooth with depth which in part is due to the choice of
parameterisation (i.e., cubic splines are smooth functions). The lateral
resolution however does not exhibit complete smoothness; instead, appears to
be degraded at some regions. This is a result of using randomly uncorrelated
data. In this case, the random data noise manifest as small-scale artifacts in
the model space.

One of the main advantages of a Bayesian framework is we can express the
solution as a marginal posterior pdf where the width of the distribution
quantifies the model uncertainties. These depth profiles were obtained by
merging the ensemble of models from the 20 Markov chains. The entire number of
models $\mathbf{m}$ used to build the posterior pdf is thus 2.0 $\times$
$10^{6}$. Results show that the profiles successfully capture the true
structure although sharp gradients fail to be resolved. The resulting
azimuthal anisotropy shows some peculiarities. Since we only make use of five
control points which are fixed, the saddle points of the true azimuthal
anisotropy fail to be captured. The increase in model uncertainty at depths is
again a result of surface waves being concentrated at shallower depths. Fig. 6
shows 1D marginal distributions versus depth of $L$, $\xi$, peak-to-peak
azimuthal anisotropy, and its fast azimuth $\Psi$ at a specific geographical
location.

(a)

(b)

(c)

(d)

Figure 5: True models (left), Mean models recovered from the inversion
(right). (A) and (B) $L-$ structure, (C) and (D) radial anisotropy.  Figure 6:
1D marginal posterior distributions of $L$, $\xi$, peak-to-peak azimuthal
anisotropy, and its fast azimuth $\Psi$ at a specific geographical location,
inferred from the Bayesian inversion of surface wave dispersion curves. The
true structures are plotted in solid red.

###  3.2 Upper mantle flow induced by subduction

We up the notch further by considering instantaneous flow in the upper mantle
induced by subduction as shown in Fig. 7. Similar to the previous example, we
added random uncorrelated noise with $\sigma_{c_{R,L}}$ = 0.001 km/s and
$\sigma_{c_{1,2}}$ = 0.0005 km/s. Fig. 8(a) shows the surface wave maps across
a subduction zone. As expected, surface waves (especially Rayleigh in this
case) undeniably map the subducting slab.

Figure 7: Snapshot of a 3D upper mantle flow induced by subduction. The model
dimensions are also of the size 400 km $\times$ 400 km $\times$ 400 km. The
dashed lines at the surface represent the 32 geographical locations of the
local surface wave dispersion measurements.

(a) Figure 8: Reconstructed 2D surface wave maps of a subduction zone with
added noise. (A) Rayleigh. (B) Love.

(b)

The inversion properties are preserved (e.g., number of independent chains,
same prior bounds, etc.). Fig. 9 shows the true structures and those recovered
by Bayesian inversion. As expected, we are still able to successfully recover
the subducting slab based from the inversions for $L$. Unlike the previous
example however, the added complexity brought about by such a type of tectonic
setting makes it more difficult for surface wave tomography. One of the
reasons for this is the choice of parameterisation as other means of
parameterising the model might increase the quality of the recovered images.
Other glaring features such as the presence of blobs in $\xi$ instead of sharp
edges ascribe to the use of long period observations. Finally, the lateral
resolution is very much degraded by random uncorrelated noise being mapped in
the parameter space.

Fig. 10 displays the 1D marginal posterior depth profiles of the model
parameters at a specific geographical location. As predicted, the sharp
gradients such as those shown in $L$ at 75 km are not recovered by surface
waves but are spatially-averaged instead. The complex azimuthal anisotropy
structure also fails to be rendered by surface wave inversion in part due to
the type of parameterisation. It is also quite evident that surface waves lose
resolution with depth as initially observed in the sphere case.

(a)

(b)

(c)

(d)

Figure 9: Subduction case. True models (left), Mean models recovered from the
inversion (right). (A) and (B) $L-$ structure, (C) and (D) radial anisotropy.
Figure 10: Subduction case. 1D marginal posterior distributions of $L$, $\xi$,
peak-to-peak azimuthal anisotropy, and its fast azimuth $\Psi$ at a specific
geographical location, inferred from the Bayesian inversion of surface wave
dispersion curves. The true structures are plotted in solid red.

##  4 Discussion and conclusion

We have demonstrated conventional anisotropic surface wave tomography. Here,
the inverse problem is fully Bayesian in that the solution is an ensemble of
models distributed according to a posterior probability. Here, the parameter
space is searched using a Markov chain Monte Carlo algorithm with importance
sampling. This method was applied to azimuthally-varying surface wave
dispersion curves computed from a 3D deforming upper mantle.

In both synthetic examples, conventional surface wave tomography was able to
recover robust features. However, the resolving power of surface waves limits
the vertical resolution of the recovered structures. As exemplified in both
cases, surface are long period observations and hence cannot resolve small-
scale features. These features are instead spatially-averaged and are smooth
as a result. It is also important to emphasize that the choice of
parameterisation regularises the inverse problem. As initially stated, the
positions of the control points where the values are perturbed are fixed. As a
consequence, the results from all the chains inflict strong dependency of the
final result on the parameterisation. Such is apparent in the recovered
azimuthal anisotropy structures. Finally, since the energy associated with
surface waves tend to be more concentrated near the surface, they thus tend to
decrease resolution with depth. This is evident in both inversions where below
the 250 km mark, the width of the posteriors begin to increase.

To handle such complications, possible future avenues include the
incorporation of higher-modes. Throughout the inversions, we only considered
fundamental-mode surface waves. Using higher modes would increase the
sensitivity of surface waves to deeper structures thus providing adequate
resolution with depth (e.g. Simons et al., 1999). Another alternative route we
could take is to consider trans-dimensional approaches where the number of
model parameters to be inverted for are treated as an unknown. In such cases,
the model adapts to the data itself, thus providing a state of balance between
model complexity and resolution (Bodin et al., 2012).

Across the horizontal, the recovered structures appear to be less resolved
since unlike teleseismic body waves, surface waves exhibit poor lateral
resolution. Still, we should expect the horizontal structures to be smooth.
This is not the case however as this was a result of using randomly
uncorrelated data noise. This random noise maps as small-scale artifacts in
the model, which clearly explains the lateral discontinuities. Since real data
noise are inherently spatially- and periodically-correlated, it is therefore
necessary to build the full data covariance matrix and account for them in the
likelihood function.

Finally, even though we placed ourselves in the best case scenario, that is,
imposing hard a priori constraints in the inversions by setting the correct
values for $A$, $\phi$, $\eta$, $B_{c}$, and $B_{s}$, we are still hindered by
the problems associated with conventional surface wave tomography. In
practice, these values are treated as unknowns in the inversion or are
determined based on empirical relations. Implementing the former complicates
the inversion procedure in a way that it increases the model complexity, and
that using only one type of data may not resolve these parameters at every
location. One is thus forced to go with the latter, where they utilise simple
relations such as the $V_{p}V_{s}$ ratio to constrain $V_{p}$-related
parameters. Velocity models inferred from such simplistic formulations however
may not be representative of the would-be recovered structure especially in
situations where complex underlying mechanisms dictate the lithological
integrity of the region. The tomographic problem should therefore be
approached from a different perspective, where in lieu of prescribing tensor
symmetries at the outset it is instead driven by key geophysical processes,
such as the work of Magali et al. (2020).

## Acknowledgements

The provision of the computational resources are made possible through the ERC
grant $-$ TRANSCALE: Reconciling Scales in Global Seismology project. The
manuscript is part of the PhD dissertation of JK Magali. All computations were
carried out through the in-house Transcale cluster situated in Lyon.

## References

  * Baumann et al. [2014] Tobias S Baumann, Boris JP Kaus, and Anton A Popov.  Constraining effective rheology through parallel joint geodynamic inversion.  _Tectonophysics_ , 631:197–211, 2014.
  * Baumann and Kaus [2015] TS Baumann and Boris JP Kaus.  Geodynamic inversion to constrain the non-linear rheology of the lithosphere.  _Geophysical Journal International_ , 202(2):1289–1316, 2015.
  * Bodin et al. [2016] T Bodin, J Leiva, B Romanowicz, V Maupin, and Huaiyu Yuan.  Imaging anisotropic layering with bayesian inversion of multiple data types.  _Geophysical Journal International_ , 206(1):605–629, 2016.
  * Bodin and Sambridge [2009] Thomas Bodin and Malcolm Sambridge.  Seismic tomography with the reversible jump algorithm.  _Geophysical Journal International_ , 178(3):1411–1436, 2009.
  * Bodin et al. [2012] Thomas Bodin, Malcolm Sambridge, H Tkalčić, Pierre Arroucau, Kerry Gallagher, and Nicholas Rawlinson.  Transdimensional inversion of receiver functions and surface wave dispersion.  _Journal of Geophysical Research: Solid Earth_ , 117(B2), 2012.
  * Debski [2010] Wojciech Debski.  Seismic tomography by monte carlo sampling.  _Pure and applied geophysics_ , 167(1-2):131–152, 2010.
  * Dziewonski and Anderson [1981] Adam M Dziewonski and Don L Anderson.  Preliminary reference earth model.  _Physics of the earth and planetary interiors_ , 25(4):297–356, 1981.
  * Hastings [1970] W Keith Hastings.  Monte carlo sampling methods using markov chains and their applications.  1970\.
  * Husen et al. [2003] S Husen, E Kissling, N Deichmann, S Wiemer, D Giardini, and M Baer.  Probabilistic earthquake location in complex three-dimensional velocity models: Application to switzerland.  _Journal of Geophysical Research: Solid Earth_ , 108(B2), 2003.
  * Lomax et al. [2000] Anthony Lomax, Jean Virieux, Philippe Volant, and Catherine Berge-Thierry.  Probabilistic earthquake location in 3d and layered models.  In _Advances in seismic event location_ , pages 101–134. Springer, 2000.
  * Magali et al. [2020] J K Magali, T Bodin, N Hedjazian, H Samuel, and S Atkins.  Geodynamic Tomography: Constraining Upper Mantle Deformation Patterns from Bayesian Inversion of Surface Waves.  _Geophysical Journal International_ , 12 2020.  ISSN 0956-540X.  doi: 10.1093/gji/ggaa577.
  * Malinverno and Briggs [2004] Alberto Malinverno and Victoria A Briggs.  Expanded uncertainty quantification in inverse problems: Hierarchical bayes and empirical bayes.  _Geophysics_ , 69(4):1005–1016, 2004.
  * Maupin and Park [2015] V Maupin and J Park.  1.09—theory and observations—seismic anisotropy.  _Treatise on Geophysics_ , pages 277–305, 2015.
  * Metropolis and Ulam [1949] Nicholas Metropolis and Stanislaw Ulam.  The monte carlo method.  _Journal of the American statistical association_ , 44(247):335–341, 1949.
  * Metropolis et al. [1953] Nicholas Metropolis, Arianna W Rosenbluth, Marshall N Rosenbluth, Augusta H Teller, and Edward Teller.  Equation of state calculations by fast computing machines.  _The journal of chemical physics_ , 21(6):1087–1092, 1953.
  * Montagner [1994] Jean-Paul Montagner.  Can seismology tell us anything about convection in the mantle?  _Reviews of Geophysics_ , 32(2):115–137, 1994\.
  * Montagner and Jobert [1988] Jean-Paul Montagner and Nelly Jobert.  Vectorial tomography—ii. application to the indian ocean.  _Geophysical Journal International_ , 94(2):309–344, 1988.
  * Montagner and Nataf [1986] Jean-Paul Montagner and Henri-Claude Nataf.  A simple method for inverting the azimuthal anisotropy of surface waves.  _Journal of Geophysical Research: Solid Earth_ , 91(B1):511–520, 1986.
  * Montagner and Nataf [1988] JP Montagner and HC Nataf.  Vectorial tomography. part i: Theory.  _Geophysical Journal International_ , 94:295–307, 1988\.
  * Morishige and Kuwatani [2020] M Morishige and T Kuwatani.  Bayesian inversion of surface heat flow in subduction zones: a framework to refine geodynamic models based on observational constraints.  _Geophysical Journal International_ , 222(1):103–109, 2020.
  * Mosegaard and Tarantola [1995] Klaus Mosegaard and Albert Tarantola.  Monte carlo sampling of solutions to inverse problems.  _Journal of Geophysical Research: Solid Earth_ , 100(B7):12431–12447, 1995.
  * Nakanishi and Anderson [1983] Ichiro Nakanishi and Don L Anderson.  Measurement of mantle wave velocities and inversion for lateral heterogeneity and anisotropy: 1. analysis of great circle phase velocities.  _Journal of Geophysical Research: Solid Earth_ , 88(B12):10267–10283, 1983.
  * Nataf et al. [1984] H-C Nataf, I Nakanishi, and Don L Anderson.  Anisotropy and shear-velocity heterogeneities in the upper mantle.  _Geophysical Research Letters_ , 11(2):109–112, 1984.
  * Ortega-Gelabert et al. [2020] O Ortega-Gelabert, S Zlotnik, JC Afonso, and P Díez.  Fast stokes flow simulations for geophysical-geodynamic inverse problems and sensitivity analyses based on reduced order modeling.  _Journal of Geophysical Research: Solid Earth_ , 125(3):e2019JB018314, 2020.
  * Ritzwoller et al. [2002] Michael H Ritzwoller, Nikolai M Shapiro, Mikhail P Barmin, and Anatoli L Levshin.  Global surface wave diffraction tomography.  _Journal of Geophysical Research: Solid Earth_ , 107(B12):ESE–4, 2002.
  * Romanowicz [2002] Barbara Romanowicz.  Inversion of surface waves: a review.  _International Geophysics Series_ , 81(A):149–174, 2002.
  * Saito [1967] Masanori Saito.  Excitation of free oscillations and surface waves by a point source in a vertically heterogeneous earth.  _Journal of Geophysical Research_ , 72(14):3689–3699, 1967.
  * Saito [1988] MASANORI Saito.  Disper80: A subroutine package for the calculation of seismic normal-mode solutions.  _Seismological algorithms_ , pages 293–319, 1988.
  * Shapiro and Ritzwoller [2002] NM Shapiro and MH Ritzwoller.  Monte-carlo inversion for a global shear-velocity model of the crust and upper mantle.  _Geophysical Journal International_ , 151(1):88–105, 2002.
  * Simons et al. [1999] Frederik J Simons, Alet Zielhuis, and Rob D Van Der Hilst.  The deep structure of the australian continent from surface wave tomography.  In _Developments in Geotectonics_ , volume 24, pages 17–43. Elsevier, 1999.
  * Smith and Dahlen [1973] Martin L Smith and FA Dahlen.  The azimuthal dependence of love and rayleigh wave propagation in a slightly anisotropic medium.  _Journal of Geophysical Research_ , 78(17):3321–3333, 1973.
  * Trampert and Woodhouse [1995] Jeannot Trampert and John H Woodhouse.  Global phase velocity maps of love and rayleigh waves between 40 and 150 seconds.  _Geophysical Journal International_ , 122(2):675–690, 1995.