AABI 20203rd Symposium on Advances in Approximate Bayesian Inference, 2020

# Preconditioned training of normalizing flows for variational inference in
inverse problems

Ali Siahkoohi<EMAIL_ADDRESS>
Georgia Institute of Technology  Gabrio Rizzuti<EMAIL_ADDRESS>
Georgia Institute of Technology Now at Utrecht University. Mathias Louboutin
<EMAIL_ADDRESS>
Georgia Institute of Technology  Philipp A. Witte<EMAIL_ADDRESS>
Microsoft  Felix J. Herrmann<EMAIL_ADDRESS>
Georgia Institute of Technology

###### Abstract

Obtaining samples from the posterior distribution of inverse problems with
expensive forward operators is challenging especially when the unknowns
involve the strongly heterogeneous Earth. To meet these challenges, we propose
a preconditioning scheme involving a conditional normalizing flow (NF) capable
of sampling from a low-fidelity posterior distribution directly. This
conditional NF is used to speed up the training of the high-fidelity objective
involving minimization of the Kullback-Leibler divergence between the
predicted and the desired high-fidelity posterior density for indirect
measurements at hand. To minimize costs associated with the forward operator,
we initialize the high-fidelity NF with the weights of the pretrained low-
fidelity NF, which is trained beforehand on available model and data pairs.
Our numerical experiments, including a 2D toy and a seismic compressed sensing
example, demonstrate that thanks to the preconditioning considerable speed-ups
are achievable compared to training NFs from scratch.

##  1 Introduction

Our aim is to perform approximate Bayesian inference for inverse problems
characterized by computationally expensive forward operators,
$F:\mathcal{X}\rightarrow\mathcal{Y}$, with a data likelihood,
$\pi_{\text{like}}(\bm{y}\mid\bm{x})$:

$\bm{y}=F(\bm{x})+\bm{\epsilon},$ (1)

where $\bm{x}\in\mathcal{X}$ is the unknown model, $\bm{y}\in\mathcal{Y}$ the
observed data, and $\bm{\epsilon}\sim\mathrm{N}(\bm{0},\sigma^{2}\bm{I})$ the
measurement noise. Given a prior density, $\pi_{\text{prior}}(\bm{x})$,
variational inference (VI, Jordan et al., 1999) based on normalizing flows
(NFs, Rezende and Mohamed, 2015) can be used where the Kullback-Leibler (KL)
divergence is minimized between the predicted and the target—i.e., _high-
fidelity_ , posterior density $\pi_{\text{post}}(\bm{x}\mid\bm{y})$ (Liu and
Wang, 2016; Kruse et al., 2019; Rizzuti et al., 2020; Siahkoohi et al., 2020;
Sun and Bouman, 2020):

$\min_{\theta}\,\mathbb{E}\,_{\bm{z}\sim\pi_{z}(\bm{z})}\bigg{[}\frac{1}{2\sigma^{2}}\left\|F\big{(}T_{\theta}(\bm{z})\big{)}-\bm{y}\right\|_{2}^{2}-\log\pi_{\text{prior}}\big{(}T_{\theta}(\bm{z})\big{)}-\log\Big{|}\det\nabla_{z}T_{\theta}(\bm{z})\Big{|}\bigg{]}.$
(2)

In the above expression, $T_{\theta}:\mathcal{Z}_{x}\to\mathcal{X}$ denotes a
NF with parameters $\bm{\theta}$ and a Gaussian latent variable
$\bm{z}\in\mathcal{Z}_{x}$. The above objective consists of the data
likelihood term, regularization on the output of the NF, and a log-determinant
term that is related to the entropy of the NF output. The last term is
necessarily to prevent the output of the NF from collapsing on the maximum a
posteriori estimate. For details regarding the derivation of the objective in
Equation (2), we refer to Appendix A. During training, we replace the
expectation by Monte-Carlo averaging using mini-batches of $\bm{z}$. After
training, samples from the approximated posterior,
$\pi_{\theta}(\bm{x}\mid\bm{y})\approx\pi_{\text{post}}(\bm{x}\mid\bm{y})$,
can be drawn by evaluating $T_{\theta}(\bm{z})$ for
$\bm{z}\sim\pi_{z}(\bm{z})$ (Kruse et al., 2019). It is important to note that
Equation (2) trains a NF specific to the observed data $\bm{y}$. While the
above VI formulation in principle allows us to train a NF to generate samples
from the posterior given a single observation $\bm{y}$, this variational
estimate requires access to a prior density, and the training calls for
repeated evaluations of the forward operator, $F$, as well as the adjoint of
its Jacobian, $\nabla{F}^{\top}$. As in multi-fidelity Markov chain Monte
Carlo (MCMC) sampling (Peherstorfer and Marzouk, 2019), the costs associated
with the forward operator may become prohibitive even though VI-based methods
are known to have computational advantages over MCMC (Blei et al., 2017).

Aside from the above computational considerations, reliance on having access
to a prior may be problematic especially when dealing with images of the
Earth’s subsurface, which are the result of complex geological processes that
do not lend themselves to be easily captured by hand-crafted priors. Under
these circumstances, data-driven priors—or even better data-driven posteriors
obtained by training over model and data pairs sampled from the joint
distribution, $\widehat{\pi}_{y,x}(\bm{y},\bm{x})$—are preferable. More
specifically, we follow Kruse et al. (2019), Kovachki et al. (2020), and
Baptista et al. (2020), and formulate the objective function in terms of a
block-triangular conditional NF,
$G_{\phi}:\mathcal{Y}\times\mathcal{X}\to\mathcal{Z}_{y}\times\mathcal{Z}_{x}$,
with latent space $\mathcal{Z}_{y}\times\mathcal{Z}_{x}$:

$\displaystyle\min_{\phi}\,\mathbb{E}\,_{\bm{y},\bm{x}\sim\widehat{\pi}_{y,x}(\bm{y},\bm{x})}\,\left[\frac{1}{2}\left\|G_{\phi}(\bm{y},\bm{x})\right\|^{2}-\log\Big{|}\det\nabla_{y,x}\,G_{\phi}(\bm{y},\bm{x})\Big{|}\right],$
(3) $\displaystyle\text{where}\quad
G_{\phi}(\bm{y},\bm{x})=\begin{bmatrix}G_{\phi_{y}}(\bm{y})\\\
G_{\phi_{x}}(\bm{y},\bm{x})\end{bmatrix},\
\bm{\phi}=\left\\{\bm{\phi}_{y},\bm{\phi}_{x}\right\\}.$

Thanks to the block-triangular structure of $G_{\phi}$, samples of the
approximated posterior,
$\pi_{\phi}(\bm{x}\mid\bm{y})\approx\pi_{\text{post}}(\bm{x}\mid\bm{y})$ can
be drawn by evaluating $G_{\phi_{x}}^{-1}(G_{\phi_{y}}(\bm{y}),\bm{z})$ for
$\bm{z}\sim\pi_{z}(\bm{z})$ (Marzouk et al., 2016). Unlike the objective in
Equation (2), training $G_{\phi}$ does not involve multiple evaluations of $F$
and $\nabla F^{\top}$, nor does it require specifying a prior density.
However, its success during inference heavily relies on having access to
training pairs from the joint distribution,
$\bm{y},\bm{x}\sim\widehat{\pi}_{y,x}(\bm{y},\bm{x})$. Unfortunately, unlike
medical imaging, where data is abundant and variability among patients is
relatively limited, samples from the joint distribution are unavailable in
geophysical applications. Attempts have been made to address this lack of
training pairs including the generation of simplified artificial geological
models (Mosser et al., 2019), but these approaches cannot capture the true
heterogeneity exhibited by the Earth’s subsurface. This is illustrated in
Figure 1, which shows several true seismic image patches drawn from the
Parihaka dataset. Even though samples are drawn from a single data set, they
illustrate significant differences between shallow (Figures 1a and 1b) and
deeper (Figures 1c and 1d) sections.

(a)

(b)

(c)

(d)

Figure 1: Subsurface images estimated from a real seismic survey, indicating
strong heterogeneity between (a), (b) shallow and (c), (d) deep parts of the
same survey area.

To meet the challenges of computational cost, heterogeneity and lack of access
to training pairs, we propose a preconditioning scheme where the two described
VI methods are combined to:

  1. 1.

take maximum advantage of available samples from the joint distribution
$\widehat{\pi}_{y,x}(\bm{y},\bm{x})$, to pretrain $G_{\phi}$ by minimizing
Equation (3). We only incur these costs once, by training this NF beforehand.
As these samples typically come from a different (neighboring) region, they
are considered as low-fidelity;

  2. 2.

exploit the invertibility of $G_{\phi_{x}}(\bm{y},\,\cdot\,)$, which gives us
access to a low-fidelity posterior density, $\pi_{\phi}(\bm{x}\mid\bm{y})$.
For a given $\bm{y}$, this trained (conditional) prior can be used in Equation
(2);

  3. 3.

initialize $T_{\theta}$ with weights from the pretrained $G_{\phi_{x}}^{-1}$.
This initialization can be considered as an instance of transfer learning
(Yosinski et al., 2014), and we expect a considerable speed-up when solving
Equation (2). This is important since it involves inverting $F$, which is
computationally expensive.

##  2 Related work

In the context of variational inference for inverse problems with expensive
forward operators, Herrmann et al. (2019) train a generative model to sample
from the posterior distribution, given indirect measurements of the unknown
model. This approach is based on an Expectation Maximization technique, which
infers the latent representation directly instead of using an inference
encoding model. While that approach allows for inclusion of hand-crafted
priors, capturing the posterior is not fully developed. Like Kovachki et al.
(2020), we also use a block-triangular map between the joint model and data
distribution and their respective latent spaces to train a network to generate
samples from the conditional posterior. By imposing an additional monotonicity
constraint, these authors train a generative adversarial network (GAN,
Goodfellow et al., 2014) to directly sample from the posterior distribution.
To allow for scalability to large scale problems, we work with NFs instead,
because they allow for more memory efficient training (Leemput et al., 2019;
Putzky and Welling, 2019; Peters et al., 2020; Peters and Haber, 2020). Our
contribution essentially corresponds to a reformulation of Parno and Marzouk
(2018) and Peherstorfer and Marzouk (2019). In that work, transport-based maps
are used as non-Gaussian proposal distributions during MCMC sampling. As part
of the MCMC, this transport map is then tuned to match the target density,
which improves the efficiency of the sampling. Peherstorfer and Marzouk (2019)
extend this approach by proposing a preconditioned MCMC sampling technique
where a transport-map trained to sample from a low-fidelity posterior
distribution is used as a preconditioner. This idea of multi-fidelity
preconditioned MCMC inspired our work where we setup a VI objective instead.
We argue that this formulation can be faster and may be easier to scale to
large-scale Bayesian inference problems (Blei et al., 2017).

Finally, there is a conceptional connection between our work and previous
contributions on amortized variational inference (Gershman and Goodman, 2014),
including an iterative refinement step (Hjelm et al., 2016; Krishnan et al.,
2018; Kim et al., 2018; Marino et al., 2018). Although similar in spirit, our
approach is different from these attempts because we adapt the weights of our
conditional generative model to account for the inference errors instead of
correcting the inaccurate latent representation of the out-of-distribution
data.

##  3 Multi-fidelity preconditioning scheme

For an observation $\bm{y}$, we define a NF
$T_{\phi_{x}}:\mathcal{Z}_{x}\to\mathcal{X}$ as

$T_{\phi_{x}}(\bm{z}):=G_{\phi_{x}}^{-1}(G_{\phi_{y}}(\bm{y}),\bm{z}),$ (4)

where we obtain $\bm{\phi}=\left\\{\bm{\phi}_{y},\bm{\phi}_{x}\right\\}$ by
training $G_{\phi}$ through minimizing the objective function in Equation (3).
To train $G_{\phi}$, we use available low-fidelity training pairs
$\bm{y},\bm{x}\sim\widehat{\pi}_{y,x}(\bm{y},\bm{x})$. We perform this
training phase beforehand, similar to the pretraining phase during transfer
learning (Yosinski et al., 2014). Thanks to the invertibility of $G_{\phi}$,
it provides an expression for the posterior. We refer to this posterior as
low-fidelity because the network is trained with often scarce and out-of-
distribution training pairs. Because the Earth’s heterogeneity does not lend
itself to be easily captured by hand-crafted priors, we argue that this NF can
still serve as a (conditional) prior in Equation (2):

$\pi_{\text{prior}}(\bm{x}):=\pi_{\bm{z}}\big{(}G_{\phi_{x}}(\bm{y},\bm{x})\big{)}\,\Big{|}\det\nabla_{x}G_{\phi_{x}}(\bm{y},\bm{x})\Big{|}.$
(5)

To train the high-fidelity NF given observed data $\bm{y}$, we minimize the KL
divergence between the predicted and the high-fidelity posterior density,
$\pi_{\text{post}}(\bm{x}\mid\bm{y})$ (Liu and Wang, 2016; Kruse et al., 2019)

$\min_{\phi_{x}}\,\mathbb{E}_{\bm{z}\sim\pi_{z}(\bm{z})}\bigg{[}\frac{1}{2\sigma^{2}}\left\|F\big{(}T_{\phi_{x}}(\bm{z})\big{)}-\bm{y}\right\|_{2}^{2}-\log\pi_{\text{prior}}\big{(}T_{\phi_{x}}(\bm{z})\big{)}-\log\Big{|}\det\nabla_{z}T_{\phi_{x}}(\bm{z})\Big{|}\bigg{]},$
(6)

where the prior density of Equation (5) is used. Notice that this minimization
problem differs from the one stated in Equation (2). Here, the optimization
involves “fine-tuning” the low-fidelity network parameters $\phi_{x}$
introduced in Equation (3). Moreover, this low-fidelity network is also used
as a prior. While other choices exist for the latter—e.g., it could be
replaced or combined with a hand-crafted prior in the form of constraints
(Peters et al., 2019) or by a separately trained data-driven prior (Mosser et
al., 2019), using the low-fidelity posterior as a prior (cf. Equation (5)) has
certain distinct advantages. First, it removes the need for training a
separate data-driven prior model. Second, use of the low-fidelity posterior
may be more informative (Yang and Soatto, 2018) than its unconditional
counterpart because it is conditioned by the observed data $\bm{y}$. In
addition, our multi-fidelity approach has strong connections with online
variational Bayes (Zeno et al., 2018) where data arrives sequentially and
previous posterior approximates are used as priors for subsequent
approximations.

In summary, the problem in Equation (6) can be interpreted as an instance of
transfer learning (Yosinski et al., 2014) for conditional NFs. This
formulation is particularly useful for inverse problems with expensive forward
operators, where access to high fidelity training samples, i.e. samples from
the target distribution, is limited. In the next section, we present two
numerical experiments designed to show the speed-up and accuracy gained with
our proposed multi-fidelity formulation.

##  4 Numerical experiments

We present two synthetic examples aimed at verifying the anticipated speed-up
and increase in accuracy of the predicted posterior density via our multi-
fidelity preconditioning scheme. The first example is a two-dimensional
problem where the posterior density can be accurately and cheaply sampled via
MCMC. The second example demonstrates the effect of the preconditioning scheme
in a seismic compressed sensing (Candes et al., 2006; Donoho, 2006) problem.
Details regarding training hyperparameters and the NF architectures are
included in Appendix B. Code to reproduce our results are made available on
GitHub. Our implementation relies on InvertibleNetworks.jl (Witte et al.,
2020), a recently-developed memory-efficient framework for training invertible
networks in the Julia programming language.

###  4.1 2D toy example

To illustrate, the advantages of working with our multi-fidelity scheme, we
consider the 2D Rosenbrock distribution,
$\pi_{\text{prior}}(\bm{x})\propto\exp\left(-\frac{1}{2}x_{1}^{2}-\left(x_{2}-x_{1}^{2}\right)^{2}\right)$,
plotted in Figure 2a. High-fidelity data $\bm{y}\in\mathbb{R}^{2}$ are
generated via $\bm{y}=\bm{A}\,\bm{x}+\bm{\epsilon}$, where
$\bm{\epsilon}\sim\mathrm{N}(\bm{0},0.4^{2}\bm{I})$ and
$\bm{A}\in\mathbb{R}^{2\times 2}$ is a forward operator. To control the
discrepancy between the low- and high-fidelity samples, we set $\bm{A}$ equal
to $\bar{\bm{A}}\mathbin{/}\rho(\bar{\bm{A}})$, where $\rho(\cdot)$ is the
spectral radius of $\bar{\bm{A}}=\bm{\Gamma}+\gamma\bm{I}$,
$\bm{\Gamma}\in\mathbb{R}^{2\times 2}$ has independent and normally
distributed entries, and $\gamma=3$. By choosing smaller values for $\gamma$,
we make $\bm{A}$ more dissimilar to the identity matrix, therefore increasing
the discrepancy between the low- and high-fidelity posterior.

Figure 2b depicts the low- (purple) and high-fidelity (red) data densities.
The dark star represents the unknown model. Low-fidelity data samples are
generated with the identity as the forward operator. During the pretraining
phase conducted beforehand, we minimize the objective function in Equation (3)
for $25$ epochs.

(a)

(b)

(c)

(d)

(e)

Figure 2: (a) Prior, (b) low- and high-fidelity data, (c) low- (blue) and
high-fidelity (orange) approximated posterior densities, and (d) approximated
posterior densities via MCMC (dark circles), and objectives in Equations (2)
in green and (6) in orange. (e) Objective value during training via Equations
(2) in green and (6) in orange.

The pretrained low-fidelity posterior is subsequently used to precondition the
minimization of (6) given observed data $\bm{y}$. The resulting low- and high-
fidelity estimates for the posterior as plotted in Figure 2c differ
significantly. In Figure (2d), the accuracy of the proposed method is verified
by comparing the approximated high-fidelity posterior density (orange
contours) with the approximation (in green) obtained by minimizing the
objective of Equation (2). The overlap between the orange contours and the
green shaded area confirms consistency between the two methods. To assess the
accuracy of the estimated densities themselves, we also include samples from
the posterior (dark circles) obtained via stochastic gradient Langevin
dynamics (Welling and Teh, 2011), an MCMC sampling technique. As expected, the
estimated posterior densities with and without the preconditioning scheme are
in agreement with the MCMC samples.

Finally, to illustrate the performance our multi-fidelity scheme, we consider
the convergence plot in Figure 2e where the objective values of Equations (2)
and (6) are compared. As explained in Appendix A, the values of the objective
functions correspond to the KL divergence (plus a constant) between the
posterior given by Equation (2) and the posterior distribution obtained by our
multi-fidelity approach (Equation (6)). As expected, the multi-fidelity
objective converges much faster because of the “warm start”. In addition, the
updates of $T_{\phi_{x}}$ via Equation (6) succeed in bringing down the KL
divergence within only five epochs (see orange curve), whereas it takes $25$
epochs via the objective in Equation (2) to reach approximately the same KL
divergence. This pattern holds for smaller values of $\gamma$ too as indicated
in Table 1. According to Table 1, the improvements by our multi-fidelity
method become more pronounced if we decrease the $\gamma$. This behavior is to
be expected since the samples used for pretraining are more and more out of
distribution in that case. We refer to Appendix C for additional figures for
different values of $\gamma$.

$\gamma$ | Low-fidelity | Without preconditioning | With preconditioning
---|---|---|---
$3$ | $6.13$ | $4.88$ | $4.66$
$2$ | $8.43$ | $5.60$ | $5.26$
$1$ | $8.84$ | $6.26$ | $6.51$
$0$ | $14.73$ | $8.41$ | $8.45$
Table 1: The KL divergence (plus some constant, see Appendix A) between
different estimated posterior densities and high-fidelity posterior
distribution for different values of $\gamma$. Second column corresponds to
the low-fidelity posterior estimate obtained via Equation (3), third column
relates to the posterior estimate via Equation (2), and finally, the last
column shows the same quantity for the posterior estimate via the multi-
fidelity preconditioned scheme.

###  4.2 Seismic compressed sensing example

This experiment is designed to show challenges with geophysical inverse
problems due to the Earth’s strong heterogeneity. We consider the inversion of
noisy indirect measurements of image patches $\bm{x}\in\mathbb{R}^{256\times
256}$ sampled from deeper parts of the Parihaka seismic dataset. The observed
measurements are given by $\bm{y}=\bm{A}\bm{x}+\bm{\epsilon}$ where
$\bm{\epsilon}\sim\mathrm{N}(\bm{0},0.2^{2}\bm{I})$. For simplicity, we chose
$\bm{A}={\bm{M}}^{T}{\bm{M}}$ with ${\bm{M}}$ a compressing sensing matrix
with $66.66\,$% subsampling. The measurement vector $\bm{y}$ corresponds to a
pseudo-recovered model contaminated with noise.

To mimic a realistic situation in practice, we change the likelihood
distribution by reducing the standard deviation of the noise to $0.01$ in
combination with using image patches sampled from the shallow part of the
Parihaka dataset. As we have seen in Figure 1, these patches differ in
texture. Given pairs $\bm{y},\bm{x}\sim\widehat{\pi}_{y,x}(\bm{y},\bm{x})$, we
pretrain our network by minimizing Equation (3). Figures 3a and 3b contain a
pair not used during pretraining. Estimates for the conditional mean and
standard deviation obtained by drawing $1000$ samples from the pretrained
conditional NF for the noisy indirect measurement (Figure 3b) are included in
Figures 3c and 3d. Both estimates exhibit the expected behavior because the
examples in Figure 3a and 3b are within the distribution. As anticipated, this
observation no longer holds if we apply this pretrained network to indirect
data depicted in Figure 3f, which is sampled from the deeper part. However,
these results are significantly improved when the pretrained network is fine-
tuned by minimizing Equation (6). After fine tuning, the fine details in the
image are recovered (compare Figures 3g and 3h). This improvement is confirmed
by the relative errors plotted in Figures 3i and 3j, as well as by the reduced
standard deviation (compare Figures 3k and 3l).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 3: Seismic compressed sensing. First row indicates the performance of
the pretrained network on data not used during pretraining. Second row is
compares the pretraining-based recovery with the accelerated scheme result.
Last row compares the recovery errors and pointwise STDs for the two
recoveries.

##  5 Conclusions

Inverse problems in fields such as seismology are challenging for several
reasons. The forward operators are complex and expensive to evaluate
numerically while the Earth is highly heterogeneous. To handle this situation
and to quantify uncertainty, we propose a preconditioned scheme for training
normalizing flows for Bayesian inference. The proposed scheme is designed to
take full advantage of having access to training pairs drawn from a joint
distribution, which for the reasons stated above is close but not equal to the
actual joint distribution. We use these samples to train a normalizing flow
via likelihood maximization leveraging the normalizing property. We use this
pretrained low-fidelity estimate for the posterior as a prior and
preconditioner for the actual variational inference on the observed data,
which minimizes the Kullback-Leibler divergence between the predicted and the
desired posterior density. By means of a series of examples, we demonstrate
that our preconditioned scheme leads to considerable speed-ups compared to
training a normalizing flow from scratch.

## References

  * Baptista et al. (2020) Ricardo Baptista, Olivier Zahm, and Youssef Marzouk.  An adaptive transport framework for joint and conditional density estimation.  _arXiv preprint arXiv:2009.10303_ , 2020.  URL https://arxiv.org/abs/2009.10303.
  * Blei et al. (2017) David M Blei, Alp Kucukelbir, and Jon D McAuliffe.  Variational inference: A review for statisticians.  _Journal of the American statistical Association_ , 112(518):859–877, 2017.
  * Bogachev, V.I. (2006) Bogachev, V.I.  _Measure Theory_.  Springer-Verlag Berlin Heidelberg, 2006.  ISBN 9783540345138.  10.1007/978-3-540-34514-5.  URL https://books.google.com/books?id=4vmOjgEACAAJ.
  * Candes et al. (2006) Emmanuel J Candes, Justin K Romberg, and Terence Tao.  Stable signal recovery from incomplete and inaccurate measurements.  _Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences_ , 59(8):1207–1223, 2006.
  * Dinh et al. (2016) Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio.  Density estimation using Real NVP.  _arXiv preprint arXiv:1605.08803_ , 2016.
  * Donoho (2006) David L Donoho.  Compressed sensing.  _IEEE Transactions on Information Theory_ , 52(4):1289–1306, 2006.  10.1109/TIT.2006.871582.
  * Gershman and Goodman (2014) Samuel Gershman and Noah Goodman.  Amortized Inference in Probabilistic Reasoning.  In _Proceedings of the Annual Meeting of the Cognitive Science Society_ , volume 36, pages 517–522, 2014.
  * Goodfellow et al. (2014) Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative Adversarial Nets.  In _Proceedings of the 27th International Conference on Neural Information Processing Systems_ , NIPS’14, pages 2672–2680, 2014.  URL http://papers.nips.cc/paper/5423-generative-adversarial-nets.pdf.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep Residual Learning for Image Recognition.  In _The IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 770–778, June 2016.  10.1109/CVPR.2016.90.  URL https://ieeexplore.ieee.org/document/7780459.
  * Herrmann et al. (2019) Felix J. Herrmann, Ali Siahkoohi, and Gabrio Rizzuti.  Learned imaging with constraints and uncertainty quantification.  In _Neural Information Processing Systems (NeurIPS) 2019 Deep Inverse Workshop_ , 12 2019.  URL https://arxiv.org/pdf/1909.06473.pdf.
  * Hjelm et al. (2016) Devon Hjelm, Russ R Salakhutdinov, Kyunghyun Cho, Nebojsa Jojic, Vince Calhoun, and Junyoung Chung.  Iterative Refinement of the Approximate Posterior for Directed Belief Networks.  In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 29, pages 4691–4699. Curran Associates, Inc., 2016.  URL https://proceedings.neurips.cc/paper/2016/file/20c9f5700da1088260df60fcc5df2b53-Paper.pdf.
  * Jordan et al. (1999) Michael I Jordan, Zoubin Ghahramani, Tommi S Jaakkola, and Lawrence K Saul.  An Introduction to Variational Methods for Graphical Models.  _Machine Learning_ , 37(2):183–233, 1999.  10.1023/A:1007665907178.
  * Kim et al. (2018) Yoon Kim, Sam Wiseman, Andrew Miller, David Sontag, and Alexander Rush.  Semi-amortized variational autoencoders.  In Jennifer Dy and Andreas Krause, editors, _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , pages 2678–2687, Stockholmsmässan, Stockholm Sweden, 10–15 Jul 2018. PMLR.  URL http://proceedings.mlr.press/v80/kim18e.html.
  * Kingma and Ba (2014) Diederik P. Kingma and Jimmy Ba.  Adam: A Method for Stochastic Optimization.  _CoRR_ , abs/1412.6980, 2014.
  * Kovachki et al. (2020) Nikola Kovachki, Ricardo Baptista, Bamdad Hosseini, and Youssef Marzouk.  Conditional Sampling With Monotone GANs.  _arXiv preprint arXiv:2006.06755_ , 2020.
  * Krishnan et al. (2018) Rahul Krishnan, Dawen Liang, and Matthew Hoffman.  On the challenges of learning with inference networks on sparse, high-dimensional data.  In Amos Storkey and Fernando Perez-Cruz, editors, _Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics_ , volume 84 of _Proceedings of Machine Learning Research_ , pages 143–151. PMLR, 09–11 Apr 2018.
  * Kruse et al. (2019) Jakob Kruse, Gianluca Detommaso, Robert Scheichl, and Ullrich Köthe.  HINT: Hierarchical Invertible Neural Transport for Density Estimation and Bayesian Inference.  _arXiv preprint arXiv:1905.10687_ , 2019.
  * Leemput et al. (2019) Sil C. van de Leemput, Jonas Teuwen, Bram van Ginneken, and Rashindra Manniesing.  MemCNN: A Python/PyTorch package for creating memory-efficient invertible neural networks.  _Journal of Open Source Software_ , 4(39):1576, 7 2019.  ISSN 2475-9066.  10.21105/joss.01576.  URL http://dx.doi.org/10.21105/joss.01576.
  * Liu and Wang (2016) Qiang Liu and Dilin Wang.  Stein Variational Gradient Descent: A General Purpose Bayesian Inference Algorithm.  In _Advances in Neural Information Processing Systems_ , volume 29, pages 2378–2386. Curran Associates, Inc., 2016.  URL https://proceedings.neurips.cc/paper/2016/file/b3ba8f1bee1238a2f37603d90b58898d-Paper.pdf.
  * Marino et al. (2018) Joseph Marino, Yisong Yue, and Stephan Mandt.  Iterative amortized inference.  _arXiv preprint arXiv:1807.09356_ , 2018.
  * Marzouk et al. (2016) Youssef Marzouk, Tarek Moselhy, Matthew Parno, and Alessio Spantini.  Sampling via measure transport: An introduction.  _Handbook of uncertainty quantification_ , pages 1–41, 2016.
  * Mosser et al. (2019) L. Mosser, O. Dubrule, and M. Blunt.  Stochastic Seismic Waveform Inversion Using Generative Adversarial Networks as a Geological Prior.  _Mathematical Geosciences_ , 84(1):53–79, 2019\.  10.1007/s11004-019-09832-6.
  * Parno and Marzouk (2018) Matthew D Parno and Youssef M Marzouk.  Transport Map Accelerated Markov Chain Monte Carlo.  _SIAM/ASA Journal on Uncertainty Quantification_ , 6(2):645–682, 2018.  10.1137/17M1134640.
  * Peherstorfer and Marzouk (2019) Benjamin Peherstorfer and Youssef Marzouk.  A transport-based multifidelity preconditioner for Markov chain Monte Carlo.  _Advances in Computational Mathematics_ , 45(5-6):2321–2348, 2019.
  * Peters and Haber (2020) Bas Peters and Eldad Haber.  Fully reversible neural networks for large-scale 3d seismic horizon tracking.  In _82nd EAGE Annual Conference & Exhibition_, pages 1–5. European Association of Geoscientists & Engineers, 2020.
  * Peters et al. (2019) Bas Peters, Brendan R Smithyman, and Felix J Herrmann.  Projection methods and applications for seismic nonlinear inverse problems with multiple constraints.  _Geophysics_ , 84(2):R251–R269, 2019.  10.1190/geo2018-0192.1.
  * Peters et al. (2020) Bas Peters, Eldad Haber, and Keegan Lensink.  Fully reversible neural networks for large-scale surface and sub-surface characterization via remote sensing.  _arXiv preprint arXiv:2003.07474_ , 2020.
  * Putzky and Welling (2019) Patrick Putzky and Max Welling.  Invert to Learn to Invert.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 32, pages 446–456. Curran Associates, Inc., 2019.  URL https://proceedings.neurips.cc/paper/2019/file/ac1dd209cbcc5e5d1c6e28598e8cbbe8-Paper.pdf.
  * Rezende and Mohamed (2015) Danilo Rezende and Shakir Mohamed.  Variational inference with normalizing flows.  volume 37 of _Proceedings of Machine Learning Research_ , pages 1530–1538. PMLR, 07–09 Jul 2015.  URL http://proceedings.mlr.press/v37/rezende15.html.
  * Rizzuti et al. (2020) Gabrio Rizzuti, Ali Siahkoohi, Philipp A. Witte, and Felix J. Herrmann.  Parameterizing uncertainty by deep invertible networks, an application to reservoir characterization.  _arXiv preprint arXiv:2004.07871_ , 4 2020.  URL https://arxiv.org/pdf/2004.07871.pdf.
  * Santambrogio (2015) Filippo Santambrogio.  Optimal Transport for Applied Mathematicians.  _Birkäuser, NY_ , 87, 2015.  10.1007/978-3-319-20828-2.
  * Siahkoohi et al. (2020) Ali Siahkoohi, Gabrio Rizzuti, Philipp A. Witte, and Felix J. Herrmann.  Faster uncertainty quantification for inverse problems with conditional normalizing flows.  Technical Report TR-CSE-2020-2, Georgia Institute of Technology, 07 2020\.  URL https://arxiv.org/abs/2007.07985.
  * Sun and Bouman (2020) He Sun and Katherine L Bouman.  Deep probabilistic imaging: Uncertainty quantification and multi-modal solution characterization for computational imaging.  _arXiv preprint arXiv:2010.14462_ , 2020.
  * Villani (2009) Cédric Villani.  _Optimal Transport: Old and New_ , volume 338.  Springer-Verlag Berlin Heidelberg, 2009.  10.1007/978-3-540-71050-9.
  * Welling and Teh (2011) Max Welling and Yee Whye Teh.  Bayesian Learning via Stochastic Gradient Langevin Dynamics.  In _Proceedings of the 28th International Conference on Machine Learning_ , ICML’11, pages 681–688, Madison, WI, USA, 2011. Omnipress.  ISBN 9781450306195.  10.5555/3104482.3104568.  URL https://dl.acm.org/doi/abs/10.5555/3104482.3104568.
  * Witte et al. (2020) Philipp Witte, Gabrio Rizzuti, Mathias Louboutin, Ali Siahkoohi, and Felix Herrmann.  InvertibleNetworks.jl: A Julia framework for invertible neural networks, November 2020.  URL https://doi.org/10.5281/zenodo.4298853.
  * Yang and Soatto (2018) Yanchao Yang and Stefano Soatto.  Conditional Prior Networks for Optical Flow.  In _Proceedings of the European Conference on Computer Vision (ECCV)_ , pages 271–287, 2018.
  * Yosinski et al. (2014) Jason Yosinski, Jeff Clune, Yoshua Bengio, and Hod Lipson.  How transferable are features in deep neural networks?  In _Proceedings of the 27th International Conference on Neural Information Processing Systems_ , NIPS’14, pages 3320–3328, 2014.  URL http://dl.acm.org/citation.cfm?id=2969033.2969197.
  * Zeno et al. (2018) Chen Zeno, Itay Golan, Elad Hoffer, and Daniel Soudry.  Task Agnostic Continual Learning Using Online Variational Bayes.  _arXiv preprint arXiv:1803.10123_ , 2018.

##  Appendix A Mathematical derivations

Let $f:\mathcal{Z}\to\mathcal{X}$ be a bijective transformation that maps a
random variable $\bm{z}\sim\pi_{z}(\bm{z})$ to $\bm{x}\sim\pi_{x}(\bm{x})$. We
can write the change of variable formula (Villani, 2009) that relates
probability density functions $\pi_{z}$ and $\pi_{x}$ in the following manner:

$\pi_{x}(\bm{x})=\pi_{z}(\bm{z})\,\Big{|}\det\nabla_{x}f^{-1}(\bm{x})\Big{|},\quad
f(\bm{z})=\bm{x},\quad\bm{x}\in\mathcal{X}.$ (7)

This relation serves as the basis for the objective functions used throughout
this paper.

###  A.1 Derivation of Equation (2)

In Equation 2, we train a bijective transformation, denoted by
$T_{\theta}:\mathcal{Z}_{x}\to\mathcal{X}$, that maps the latent distribution
$\pi_{z_{x}}(\bm{z}_{x})$ to the high-fidelity posterior density
$\pi_{\text{post}}(\bm{x}\mid\bm{y})$. We optimize the parameters of
$T_{\theta}$ by minimizing the KL divergence between the push-forward density
(Bogachev, V.I., 2006), denoted by
$\pi_{\theta}(\,\cdot\mid\bm{y}):=T_{\sharp}\,\pi_{z_{x}}$, and the posterior
density:

$\displaystyle\mathop{\rm
arg\,min}_{\theta}\,\mathbb{D}_{\text{KL}}\left(\pi_{\theta}(\,\cdot\mid\bm{y})\mid\mid\pi_{\text{post}}(\,\cdot\mid\bm{y})\right)$
(8) $\displaystyle=$ $\displaystyle\mathop{\rm
arg\,min}_{\theta}\,\mathbb{E}_{\bm{x}\sim\pi_{\theta}(\bm{x}\mid\bm{y})}\Big{[}-\log\pi_{\text{post}}(\bm{x}\mid\bm{y})+\log\pi_{\theta}(\bm{x}\mid\bm{y})\Big{]}.$

In the above expression, we can rewrite the expectation with respect to
$\pi_{\theta}(\bm{x}\mid\bm{y})$ as the expectation with respect to the latent
distribution, followed by a mapping via $T_{\theta}$—i.e.,

$\mathop{\rm
arg\,min}_{\theta}\,\mathbb{E}_{\bm{z}_{x}\sim\pi_{z_{x}}(\bm{z})}\Big{[}-\log\pi_{\text{post}}(T_{\theta}(\bm{z}_{x})\mid\bm{y})+\log\pi_{\theta}(T_{\theta}(\bm{z}_{x})\mid\bm{y})\Big{]}.$
(9)

The last term in the expectation above can be further simplified via the
change of variable formula in Equation (7). If
$\bm{x}=T_{\theta}(\bm{z}_{x})$, then:

$\pi_{\theta}(\bm{x}\mid\bm{y})=\pi_{z_{x}}(\bm{z}_{x})\,\Big{|}\det\nabla_{x}T_{\theta}^{-1}(\bm{x})\Big{|}=\pi_{z_{x}}(\bm{z}_{x})\,\Big{|}\det\nabla_{z_{x}}T_{\theta}(\bm{z}_{x})\Big{|}^{-1}.$
(10)

The last equality in Equation (10) holds due to the invertibility of
$T_{\theta}$ and the differentiability of its inverse (inverse function
theorem). By combining Equations (9) and (10), we arrive at the following
objective function for training $T_{\theta}$:

$\mathop{\rm
arg\,min}_{\theta}\,\mathbb{E}_{\bm{z}_{x}\sim\pi_{z_{x}}(\bm{z})}\Big{[}-\log\pi_{\text{post}}(T_{\theta}(\bm{z}_{x})\mid\bm{y})+\log\pi_{z_{x}}(\bm{z}_{x})-\log\Big{|}\det\nabla_{z_{x}}T_{\theta}(\bm{z}_{x})\Big{|}\Big{]}.$
(11)

Finally, by ignoring the $\log\pi_{z_{x}}(\bm{z}_{x})$ term, which is constant
with respect to $\bm{\theta}$, using Bayes’ rule, and inserting our data
likelihood model from Equation (1), we derive Equation (2):

$\min_{\theta}\,\mathbb{E}\,_{\bm{z}_{x}\sim\pi_{z_{x}}(\bm{z}_{x})}\bigg{[}\frac{1}{2\sigma^{2}}\left\|F\big{(}T_{\theta}(\bm{z}_{x})\big{)}-\bm{y}\right\|_{2}^{2}-\log\pi_{\text{prior}}\big{(}T_{\theta}(\bm{z}_{x})\big{)}-\log\Big{|}\det\nabla_{z_{x}}T_{\theta}(\bm{z}_{x})\Big{|}\bigg{]}.$
(12)

Next, based on this equation, we derive the objective function used in the
pretraining phase.

###  A.2 Derivation of Equation (3)

The derivation of objective in Equation (3) follows directly from the change
of variable formula in Equation 7, applied to a bijective map,
$G^{-1}_{\phi}:\mathcal{Z}_{y}\times\mathcal{Z}_{x}\to\mathcal{Y}\times\mathcal{X}$,
where $\mathcal{Z}_{y}$ and $\mathcal{Z}_{y}$ are Gaussian latent spaces. That
is to say:

$\begin{split}\widehat{\pi}_{y,x}(\bm{y},\bm{x})=\pi_{z_{y},z_{x}}(\bm{z}_{y},\bm{z}_{x})\,\Big{|}\det\nabla_{y,x}G_{\phi}(\bm{y},\bm{x})\Big{|},\quad
G_{\phi}(\bm{y},\bm{x})=\begin{bmatrix}\bm{z}_{y}\\\
\bm{z}_{x}\end{bmatrix}.\end{split}$ (13)

Given (low-fidelity) training pairs,
$\bm{y},\bm{x}\sim\widehat{\pi}_{y,x}(\bm{y},\bm{x})$, the maximum likelihood
estimate of $\bm{\phi}$ is obtained via the following objective:

$\displaystyle\mathop{\rm
arg\,max}_{\phi}\,\mathbb{E}\,_{\bm{y},\bm{x}\sim\widehat{\pi}_{y,x}(\bm{y},\bm{x})}\,\left[\log\widehat{\pi}_{y,x}(\bm{y},\bm{x})\right]$
(14) $\displaystyle=\mathop{\rm
arg\,min}_{\phi}\,\mathbb{E}\,_{\bm{y},\bm{x}\sim\widehat{\pi}_{y,x}(\bm{y},\bm{x})}\,\left[-\log\pi_{z_{y},z_{x}}(\bm{z}_{y},\bm{z}_{x})-\log\Big{|}\det\nabla_{y,x}\,G_{\phi}(\bm{y},\bm{x})\Big{|}\right]$
$\displaystyle=\mathop{\rm
arg\,min}_{\phi}\,\mathbb{E}\,_{\bm{y},\bm{x}\sim\widehat{\pi}_{y,x}(\bm{y},\bm{x})}\,\left[\frac{1}{2}\left\|G_{\phi}(\bm{y},\bm{x})\right\|^{2}-\log\Big{|}\det\nabla_{y,x}\,G_{\phi}(\bm{y},\bm{x})\Big{|}\right],$

that is the objective function in Equation (3). The NF trained via the
objective function, given samples from the latent distribution, draws samples
from the low-fidelity joint distribution, $\widehat{\pi}_{y,x}$.

By construction, $G_{\phi}$ is a block-triangular map—i.e.,

$\begin{split}\quad
G_{\phi}(\bm{y},\bm{x})=\begin{bmatrix}G_{\phi_{y}}(\bm{y})\\\
G_{\phi_{x}}(\bm{y},\bm{x})\end{bmatrix},\
\bm{\phi}=\left\\{\bm{\phi}_{y},\bm{\phi}_{x}\right\\}.\end{split}$ (15)

Kruse et al. (2019) showed that after solving the optimization problem in
Equation (3), $G_{\phi}$ approximates the well-known triangular Knothe-
Rosenblat map (Santambrogio, 2015). As shown in Marzouk et al. (2016), the
triangular structure and $G_{\phi}$’s invertibility yields the following
property,

$\left(G_{\phi_{x}}^{-1}(G_{\phi_{y}}(\bm{y}),\cdot\,)\right)_{\sharp}\,\pi_{z_{x}}=\widehat{\pi}_{\text{post}}(\,\cdot\mid\bm{y}),$
(16)

where $\widehat{\pi}_{\text{post}}$ denotes the low-fidelity posterior
probability density function. The expression above means we can get access to
low-fidelity posterior distribution samples by simply evaluating
$G_{\phi_{x}}^{-1}(G_{\phi_{y}}(\bm{y}),\bm{z}_{x})$ for
$\bm{z}_{x}\sim\pi_{z_{x}}(\bm{z}_{x})$ for a given observed data $\bm{y}$.

##  Appendix B Training details and network architectures

For our network architecture, we adapt the recursive coupling blocks proposed
by Kruse et al. (2019), which use invertible coupling layers from Dinh et al.
(2016) in a hierarchical way. In other words, we recursively divide the
incoming state variables and apply an affine coupling layer. The final
architecture is obtained by composing several of these hierarchical coupling
blocks. The hierarchical structure leads to dense triangular Jacobian, which
is essential in representation power of NFs (Kruse et al., 2019).

For all examples in this paper, we use the hierarchical coupling blocks as
described in Kruse et al. (2019). The affine coupling layers within each
hierarchal block contain a residual block as described in He et al. (2016).
Each residual block has the following dimensions: $64$ input, $128$ hidden,
and $64$ output channels, except for the first and last coupling layer where
we have $4$ input and output channels, respectively. We use the Wavelet
transform and its transpose before feeding seismic images into the network and
after the last layer of the NFs.

Below, we describe the network architectures and training details regarding
the two numerical experiments described in the paper. Throughout the
experiments, we use the Adam optimization algorithm (Kingma and Ba, 2014).

###  B.1 2D toy example

We use 8 hierarchal coupling blocks, as described above for both
$G_{\phi_{x}}$ and $G_{\phi_{y}}$ (Equation (3)). As a result, due to our
proposed method in Equation (4), we choose the same architecture for
$T_{\phi_{x}}$ (Equation (2)).

For pretraining $G_{\theta}$ according to Equation (3), we use $5000$ low-
fidelity joint training pairs,
$\bm{y},\bm{x}\sim\widehat{\pi}_{y,x}(\bm{y},\bm{x})$. We minimize Equation
(3) for $25$ epochs with a batch size of $64$ and a (starting) learning rate
of $0.001$. We decrease the learning rate each epoch by a factor of $0.9$.

For the preconditioned step—i.e., solving Equation 6, we use $1000$ latent
training samples. We train for $5$ epochs with a batch size of $64$ and a
learning rate $0.001$. Finally, as a comparison, we solve the objective in
Equation 6 for a randomly initialized NF with the same $1000$ latent training
samples for $25$ epochs. We decrease the learning rate each epoch by $0.9$.

###  B.2 Seismic compressed sensing

We use 12 hierarchal coupling blocks, as described above for both
$G_{\phi_{x}}$, $G_{\phi_{y}}$, and we use the same architecture for
$T_{\phi_{x}}$ as $G_{\phi_{x}}$.

For pretraining $G_{\theta}$ according to Equation (3), we use $5282$ low-
fidelity joint training pairs,
$\bm{y},\bm{x}\sim\widehat{\pi}_{y,x}(\bm{y},\bm{x})$. We minimize Equation
(3) for $50$ epochs with a batch size of $4$ and a starting learning rate of
$0.001$. Once again, we decrease the learning rate each epoch by $0.9$.

For the preconditioned step—i.e., solving Equation 6, we use $1000$ latent
training samples. We train for $10$ epochs with a batch size $4$ and a
learning rate of $0.0001$, where we decay the step by $0.9$ in every $5$th
epoch.

##  Appendix C 2D toy example—more results

Here we show the effect $\gamma$ on our proposed method in the 2D toy
experiment. By choosing smaller values for $\gamma$, we make
$\bar{\bm{A}}\mathbin{/}\rho(\bar{\bm{A}})$ with
$\bar{\bm{A}}=\bm{\Gamma}+\gamma\bm{I}$ less close to the identity matrix,
hence enhancing the discrepancy between the low- and high-fidelity posterior.
The first row of Figure 4 shows the low- (purple) and high-fidelity (red) data
densities for decreasing values of $\gamma$ from $2$ down to $0$. The second
row depicts the predicted posterior densities via the preconditioned scheme
(orange contours) and the low-fidelity posterior in green along with MCMC
samples (dark circles). The third row compares the preconditioned posterior
densities to samples obtained via the low-fidelity pretrained NF—i.e.,
Equation (3). Finally, the last row shows the objective function values during
training with (orange) and without (green) preconditioning.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 4: 2D toy example with decreasing values of $\gamma=2,1,0$ from first
to last column, respectively. First row: Low- and high-fidelity data. Second
row: approximated posterior densities via MCMC (dark circles), and objectives
in Equations (2) and (6). Third row: overlaid prior density with predicted
posterior densities via objectives in Equations (3) and (6). Last row:
training objective values during training via Equations (2) and (6).

We observe that by decreasing $\gamma$ from $2$ to $0$, the low-fidelity
posterior approximations become worse. As a result, the objective function for
the preconditioned approach (orange) at the beginning start from a higher
value, indicating more mismatch between low- and high-fidelity posterior
densities. Finally, our preconditioning method consistently improves upon low-
fidelity posterior by training for $5$ epochs.

##  Appendix D Seismic compressed sensing—more results

Here we show more examples to verify the pretraining phase obtained via
solving the objective in Equation (2). Each row in Figure 5 corresponds to a
different testing image. The first column shows the different true seismic
images used to create low-fidelity compressive sensing data, depicted in the
second column. The third and last columns correspond to the conditional mean
and pointwise STD estimates, respectively. Clearly, the pretrained network is
able to successfully recover the true image, and consistently indicates more
uncertainty in areas with low-amplitude events.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

Figure 5: Seismic compressed sensing for four different low-fidelity images.
First column: true seismic images. Second column: low-fidelity observed data.
Third and last columns: conditional mean and pointwise STD estimates obtained
by drawing $1000$ samples from the pretrained conditional NF.


# Field Dispersion and Strong Coupling of Nuclear–Electron Spin Excitation in
MnCO3

Takahiko Makiuchi<EMAIL_ADDRESS>RIKEN Center for Emergent Matter
Science (CEMS), Wako 351–0198, Japan.  Department of Applied Physics, the
University of Tokyo, Tokyo 113-8656, Japan.  Takashi Kikkawa  Department of
Applied Physics, the University of Tokyo, Tokyo 113-8656, Japan.  Thanaporn
Sichanugrist  Department of Applied Physics, the University of Tokyo, Tokyo
113-8656, Japan.  Department of Physics, the University of Tokyo, Tokyo
113-8656, Japan.  Junki Numata  Department of Applied Physics, the University
of Tokyo, Tokyo 113-8656, Japan.  Masaki Imai  Advanced Science Research
Center, Japan Atomic Energy Agency, Tokai 319-1195, Japan.  Hiroyuki Chudo
Advanced Science Research Center, Japan Atomic Energy Agency, Tokai 319-1195,
Japan.  Saburo Takahashi  Advanced Institute for Materials Research, Tohoku
University, Sendai 980-8577, Japan.  Eiji Saitoh  RIKEN Center for Emergent
Matter Science (CEMS), Wako 351–0198, Japan.  Department of Applied Physics,
the University of Tokyo, Tokyo 113-8656, Japan.  Advanced Science Research
Center, Japan Atomic Energy Agency, Tokai 319-1195, Japan.  Advanced Institute
for Materials Research, Tohoku University, Sendai 980-8577, Japan.  Institute
for AI and Beyond, the University of Tokyo, Tokyo 113-8656, Japan

###### Abstract

Hybridized nuclear and electron spin excitation in a MnCO3 crystal, a weakly-
anisotropic antiferromagnet, has been investigated. In this material, the
hyperfine interaction is strong enough to form a nuclear spin wave. We measure
the microwave absorption by a bulk MnCO3 and observe the dispersion
representing strong frequency repulsion between electron and nuclear modes due
to their hybridization, the signature of nuclear spin wave. Additionally, we
observe that the nuclear spin resonance enters a nonlinear regime above a
certain excitation power, attributed to the excitation of finite wavenumber
nuclear spin waves.

nuclear spin, electron spin, antiferromagnet, hyperfine interaction,
Suhl–Nakamura interaction, Dzyaloshinskii–Moriya interaction, magnetization
dynamics, spintronics

††preprint: v3

In condensed matter, nuclear spins exist in isolation and interact with
neighboring electron spins through the hyperfine interaction. The hyperfine
interaction has played a pivotal role in exploring electron spin states
through nuclear magnetic resonance (NMR) [1]. Recently, the field of
spintronics has been extended into the realm of the nucleus, a concept termed
nuclear spintronics [2, 3, 4, 5, 6, 7]. At the core of nuclear spintronics is
the nuclear spin wave enabled by the Suhl–Nakamura interaction—an indirect
interaction among nuclear spins facilitated by the exchange interaction
between electron spins and the hyperfine interaction [8, 9, 10]. The nuclear
spin wave operates within a megahertz frequency range [2], and its fluctuation
persists even at very low temperatures [3], which distinguishes it from the
electron counterpart. The conversion of nuclear spin waves into electron spin
and charge currents may offer unique advantages such as extension of the
frequency range and high coherence regime.

Materials suitable for nuclear spintronics include weakly-anisotropic
antiferromagnets with strong hyperfine interactions [11], such as manganese
carbonate MnCO3 [2, 3]. The angular momenta in such materials can travel as
electron and nuclear spin waves. In this study, we clarify the field
dispersion for the hybridized nuclear and electron spin excitations in MnCO3
as a function of temperature by using broadband microwave spectroscopy. The
observed dispersions of electron and nuclear spins show the impact of the
hybridization and a high cooperativity. Examining the temperature dependence
in these dispersions, we determine various material parameters of our MnCO3
sample. We also found a nonlinear regime in the nuclear spin wave dispersion
indicating a formation of finite wavenumber nuclear spin waves.

Figure 1: (a) Crystal structure of MnCO3. (b) Schematic illustration of the
hybridized electron and nuclear spin resonance. (c) Dispersion relations of
the electron and nuclear spin excitations in MnCO3 [Eqs. (1–4)]. The colored
curves show hybridized dispersions, while the dashed curves show dispersions
without hybridization.

MnCO3 is an insulator whose crystalline structure belongs to the trigonal
crystal system [12, 13, 14]. The rhombohedral unit cell is displayed in Fig.
1(a). A total electron spin $S=5/2$ and a nuclear spin $I=5/2$ (natural
isotopic abundance of 55Mn is 100%) are located at each Mn site and couple
together through the on-site hyperfine interaction. The electron spins order
into the antiferromagnetic state with a finite canting angle due to the
Dzyaloshinskii–Moriya (DM) interaction [15, 16]. The electron spin
orientations at different sublattices can be discussed in terms of various
internal fields: the exchange field $H_{\mathrm{E}}$, the DM field
$H_{\mathrm{D}}$, the hard axis anisotropy field $H_{\mathrm{K}}$ ([111]
direction), the easy axis anisotropy field $H_{\mathrm{K^{\prime}}}$
([1$\bar{2}$1] direction) and the hyperfine field from the nuclear spin
$H_{\mathrm{hf}}$ [17]. Meanwhile, the nuclear spins are paramagnetically
polarized along the electron spin direction by the strong on-site hyperfine
field of $\mu_{0}H_{\mathrm{n}}\approx 60$ T, where $\mu_{0}$ is the vacuum
permeability. The orientations of electron and nuclear magnetic moments,
$\mathbf{M}_{i}\propto-\mathbf{S}_{i}$ and
$\mathbf{m}_{i}\propto\mathbf{I}_{i}$, are as shown in Fig. 1(b). The
subscripts $i=1,\ 2$ represent the sublattice numbers. Here, we define the
coordinate as $x\parallel[1\bar{2}1]$ (easy axis), $y\parallel[111]$ (hard
axis), and $z\parallel[10\bar{1}]$. We apply an external field $H$ along the
$z$ direction. Large $H_{\mathrm{E}}$ and small $H_{\mathrm{K^{\prime}}}$
orient $\mathbf{M}_{i}$ along the $x$ axis, and $H+H_{\mathrm{D}}$ deflects
$\mathbf{M}_{i}$ toward the $z$ direction. The canting angle of
$\mathbf{M}_{i}$ from the $x$ axis is
$\psi\approx\mathrm{arcsin}(H+H_{\mathrm{D}})/(2H_{\mathrm{E}}+H_{\mathrm{K^{\prime}}})$
under the approximation of $H_{\mathrm{E}}\gg$ (Other fields).

Two types of antiferromagnetic electron spin resonance modes can be excited
[18, 4]. One is a precession of the net magnetization called the weak-
ferromagnetic resonance (in-phase mode) and the other is a precession of the
Néel vector (out-of-phase mode). We focus on the weak-ferromagnetic resonance
which significantly hybridizes with the nuclear spin excitation due to its
smaller resonance gap. From a simplified model for MnCO3 [17], one can predict
angular frequencies of uncoupled electron and nuclear spin excitation modes,

$\displaystyle\omega_{\mathrm{e0}}^{2}(k)$ $\displaystyle\simeq$
$\displaystyle\omega_{H}(\omega_{H}+\omega_{\mathrm{D}})+2\omega_{\mathrm{E}}(\omega_{\mathrm{K^{\prime}}}+\omega_{\mathrm{hf}})+\Gamma^{2}(k),$
(1) $\displaystyle\omega_{\mathrm{n0}}$ $\displaystyle\simeq$
$\displaystyle\mu_{0}\gamma_{\mathrm{n}}H_{\mathrm{n}},$ (2)

where $\omega_{H}=\mu_{0}\gamma_{\mathrm{e}}H$,
$\omega_{\mathrm{D}}=\mu_{0}\gamma_{\mathrm{e}}H_{\mathrm{D}}$,
$\omega_{\mathrm{E}}=\mu_{0}\gamma_{\mathrm{e}}H_{\mathrm{E}}$,
$\omega_{\mathrm{K^{\prime}}}=\mu_{0}\gamma_{\mathrm{e}}H_{\mathrm{K^{\prime}}}$,
$\omega_{\mathrm{hf}}=\mu_{0}\gamma_{\mathrm{e}}H_{\mathrm{hf}}$, $k$ is the
wavenumber of magnon, $\gamma_{\mathrm{e}}$ is the gyromagnetic ratio of
electron, and $\gamma_{\mathrm{n}}=2\pi\times 10.553$ MHz/T is the
gyromagnetic ratio of 55Mn nucleus in MnCO3 [19].
$\Gamma^{2}(k)=(1-\gamma_{k})(\mu_{0}\gamma_{\mathrm{e}}H_{\mathrm{E}})^{2}$
is a $k$-dependent term with the form factor
$\gamma_{k}=\frac{1}{z}\sum_{\bm{\delta}}e^{i\bm{\delta}\cdot\mathbf{k}}\simeq
1-(ak)^{2}/6$, the location of nearest neighbors $\bm{\delta}$, and the
lattice constant $a=0.4768$ nm [13, 14]. The hybridized electron and nuclear
spin excitation modes read

$\displaystyle\omega_{\mathrm{e}}^{2}(k)$ $\displaystyle\simeq$
$\displaystyle\omega_{\mathrm{e0}}^{2}(k)+G^{2},$ (3)
$\displaystyle\omega_{\mathrm{n}}^{2}(k)$ $\displaystyle\simeq$
$\displaystyle\omega_{\mathrm{n0}}^{2}-G^{2},$ (4)

where

$G=\sqrt{\frac{2\omega_{\mathrm{E}}\omega_{\mathrm{hf}}\omega_{\mathrm{n0}}^{2}}{\omega_{\mathrm{e0}}^{2}(k)-\omega_{\mathrm{n0}}^{2}}}$
(5)

is a detuning parameter due to the hybridization by the transverse dynamical
components of $\mathbf{S}$ and $\mathbf{I}$. The coupling diminishes as
$\omega_{\mathrm{hf}}$ approaches zero ($G\rightarrow 0$). These
uncoupled/hybridized electron and nuclear spin excitation modes, described by
Eqs. (1)–(4), are plotted in Fig. 1(c) using the material parameters obtained
in this study (as shown later). For smaller $k$, the uncoupled
$\omega_{\mathrm{e0}}(k)$ approaches $\omega_{\mathrm{n0}}$, and
$\omega_{\mathrm{e}}(k)$ and $\omega_{\mathrm{n}}(k)$ repel each other due to
the hybridization.

We experimentally obtain the hybridized nuclear and electron spin dispersions
in MnCO3 using a broadband microwave spectroscopy technique. A bulk MnCO3
crystal (3 $\times$ 3 $\times$ 0.5 mm3) is set on a coplanar waveguide. The
sample is located in a cryostat with a superconducting magnet to apply $H$ to
the $z\parallel$ [10$\bar{1}$] direction. By using a network analyzer, we
irradiate a microwave at the angular frequency $\omega$ and the input power
$P_{\mathrm{in}}=-5$ dBm to the coplanar waveguide to induce a microwave
magnetic field parallel to the $x\parallel$ [1$\bar{2}$1] direction and
measure the ratio of absorbed and input microwave powers
$P_{\mathrm{abs}}/P_{\mathrm{in}}$. We designed two coplanar waveguides. One
is used for most of the measurements and the other improved one is used for
the main data in Fig. 2(a).

Figure 2: (a) Observed microwave absorption by MnCO3 at 1.8 K. The displayed
$P_{\mathrm{abs}}$ below 2 GHz is multiplied by 10 for clarity. (b) Observed
microwave absorption at 1.8 K in the low-frequency region. Figure 3: (a)
Microwave absorptions at different temperatures for MnCO3. (b) Magnified view
of the microwave absorptions in the low-frequency region. (c) Temperature
dependence of the nuclear spin excitation frequency at a large magnetic field
(1.2 T), i.e.
$\omega_{\mathrm{n0}}\approx\lim_{H\rightarrow\infty}\omega_{\mathrm{n}}$. The
arrow indicates the Néel temperature $T_{\mathrm{N}}=34.4$ K. (d) Temperature
dependence of the pulled frequency of the nuclear branch
$\delta_{\mathrm{n}}\equiv\omega_{\mathrm{n0}}-\omega_{\mathrm{n}}(H=0)$. (e)
Temperature dependence of the electron spin excitation gap
$\delta_{\mathrm{e}}\equiv\omega_{\mathrm{e}}(H=0)$.

Figures 2(a) and (b) show the observed $P_{\mathrm{abs}}/P_{\mathrm{in}}$
representing magnetic resonance in MnCO3 at 1.8 K, well below the Néel
temperature of MnCO3. The absorption above 6 GHz is due to the hybridized
electron spin excitation mode, being almost linearly proportional to the
external field $H$. The multiple lines can be attributed to magnetostatic
standing wave modes in the thickness direction of the sample [20, 21]. The
hybridized nuclear spin excitation mode appears in the low-frequency region in
Fig. 2(a). Figure 2(b) is a magnified view around the low-frequency region of
the spectrum shown in Fig. 2(a). The strong frequency drop near $H=0$ is the
signature of the nuclear spin wave formation known as the frequency pulling
effect [10, 17]. The nuclear branch approaches the uncoupled value
$\omega_{\mathrm{n0}}/2\pi\approx$ 640 MHz at large $H$ as the hybridization
weakens due to the increasing frequency difference between nuclear and
electron spin excitations. The linewidth is small at large $H$, whereas it
broadens at small $H$. The increase of linewidth at $H\approx 0$ is attributed
to the significant hybridization with the dissipative electron spin
excitation. Usually, the nuclear spin excitation is invisible in broadband
microwave spectroscopy due to the small nuclear magnetic moment. The large
absorption of the nuclear branch can be attributed to the strong hybridization
with a large electron magnetic moment.

We estimate the damping coefficients, coupling constant, and cooperativity for
the coupling between nuclear and electron spins from the data in Fig. 2. The
damping coefficients of $\kappa_{\mathrm{e}}/2\pi\approx 200$ MHz for the
electron spin excitation mode and $\kappa_{\mathrm{n}}/2\pi\approx 1$ MHz for
the nuclear spin excitation mode are taken from the asymptotic linewidths at
large $H$. We assumed that the four magnetostatic modes are equally separated
and have the same linewidth. The linewidths at $H\approx 0$ are 200 MHz for
the electron spin excitation mode and 100 MHz for the nuclear spin excitation
mode, where the broadening of the nuclear branch can be attributed to the
hybridization. We estimate the coupling constant $g$ from a model Hamiltonian
$\mathcal{H}/\hbar=\omega_{a}a^{\dagger}a+\omega_{b}b^{\dagger}b+g(a^{\dagger}b+ab^{\dagger})$,
where $\omega_{a}$, $a^{\dagger}$ and $a$ are the uncoupled frequency,
creation and annihilation operators for the electron spin excitation,
$\omega_{b}$, $b^{\dagger}$ and $b$ are those for the nuclear spin excitation,
and $\hbar$ is the reduced Planck constant. Diagonalizing $\mathcal{H}$ gives

$\omega_{\mathrm{e}},\omega_{\mathrm{n}}=\frac{1}{2}(\omega_{a}+\omega_{b})\pm\frac{1}{2}\sqrt{(\omega_{a}-\omega_{b})^{2}+4g^{2}}.$
(6)

Using Eqs. (1)–(4), $\omega_{a}=\omega_{\mathrm{e0}}$, and
$\omega_{b}=\omega_{\mathrm{n0}}$, we have $g/2\pi=1.2$ GHz at 1.7 K and 0 T.
As a result, the electron and nuclear spins in MnCO3 are in the strong
coupling regime ($g>\kappa_{\mathrm{e}},\kappa_{\mathrm{n}}$) [22] with the
cooperativity $C=g^{2}/\kappa_{\mathrm{e}}\kappa_{\mathrm{n}}\approx 7000$.
This value is greater than those of magnon–magnon ($C\sim 10$) [23, 24] and
magnon–phonon ($C\sim 3000$) [25] couplings, but is comparable to those of
millimeter-sized cavity–magnon systems [26, 27, 22, 28, 29, 30, 31].

Figures 3(a) and (b) show the microwave absorption spectra at different
temperatures $T$. The electron branch shifts toward lower frequencies as $T$
increases. The frequency pulling and absorption intensity of the hybridized
nuclear spin excitation mode weaken as $T$ increases. The dark spot at 6 GHz
and horizontal lines are attributed to the effect of microwave standing waves
in the coplanar waveguide.

The nuclear spin excitation frequency $\omega_{\mathrm{n0}}$ in Fig. 3(c)
directly gives the electron magnetic order parameter $\langle S_{z}\rangle$
through $\omega_{\mathrm{n0}}=A\langle S_{z}\rangle/\hbar$, where $A$ is the
hyperfine coefficient. A fitting of $\omega_{\mathrm{n0}}$ according to
$\langle S_{z}\rangle\propto[1-T/T_{\mathrm{N}}]^{\beta}$ gives the Néel
temperature $T_{\mathrm{N}}=34.4$ K and the critical exponent $\beta=0.310(3)$
for the antiferromagnetic ordering. The relatively high $T_{\mathrm{N}}$
implies a good quality of the sample [14]. The critical exponent agrees with
values from experiments on weakly anisotropic antiferromagnets [14, 32, 33].

The pulling frequency
$\delta_{\mathrm{n}}\equiv\omega_{\mathrm{n0}}-\omega_{\mathrm{n}}(H=0)$
represents the bandwidth on the nuclear spin wave dispersion. The temperature
dependence of $\delta_{\mathrm{n}}$, plotted in Fig. 3(d), represents that the
hybridization rapidly increases below $\sim 10$ K, driven by the growing
nuclear spin polarization $\langle I_{z}\rangle\propto 1/T$. The gap of the
electron mode $\delta_{\mathrm{e}}$ accordingly increases in Fig. 3(e).
Although the nuclear spin polarization is $\langle I_{z}\rangle/I\approx$ 2%
at $T=1.8$ K, the pulled frequency of the nuclear branch
$\delta_{\mathrm{n}}/2\pi\equiv[\omega_{\mathrm{n0}}-\omega_{\mathrm{n}}(H=0)]/2\pi\approx
400$ MHz is about 60% of the uncoupled nuclear spin excitation frequency
$\omega_{\mathrm{n0}}/2\pi=$ 640 MHz. This large shift implies strong
hybridization due to the exchange amplification in weakly anisotropic
antiferromagnet [34, 35, 11, 36]. This effect is described in the detuning
parameter $G\propto\sqrt{2\omega_{\mathrm{E}}\omega_{\mathrm{hf}}}$ with large
$\omega_{\mathrm{E}}/\pi=1.87$ THz and $\omega_{\mathrm{hf}}\propto\langle
I_{z}\rangle$. The full $T$-dependence data give the internal field values:
$\mu_{0}H_{\mathrm{E}}=33.4$ T, $\mu_{0}H_{\mathrm{D}}=0.461$ T,
$\mu_{0}H_{\mathrm{K^{\prime}}}=0.03$ mT, $\mu_{0}H_{\mathrm{hf}}=(1.1\
\mathrm{mT})/(T/\mathrm{K})$, and $\mu_{0}H_{\mathrm{n}}=60.65$ T.

We found a nonlinear regime of the hybridized nuclear spin excitation mode
under the large input microwave power $P_{\mathrm{in}}$. Figure 4(a) depicts
the microwave absorption at different $P_{\mathrm{in}}$ and $T$. With
$P_{\mathrm{in}}=-5$ dBm and less, the nuclear branch obeys
$\omega_{\mathrm{n}}(k=0)$ in Eq. (4). By contrast, with $P_{\mathrm{in}}=0$
and $+5$ dBm, the absorption spread toward higher frequencies. The spread
appears to be bounded by the upper limit at $\omega_{\mathrm{n0}}/2\pi$ = 640
MHz and the lower limit of $\omega_{\mathrm{n}}(k=0)$. The nonlinear
absorption is prominent at lower $T$’s, at which the hybridization is strong.
Note that we swept the frequency from 0.66 to 0.50 GHz with a long dwell time
to prevent history effects.

Figure 4: (a) Observed nonlinear microwave absorption at different input
microwave powers and temperatures. The arrow indicates the zero-absorption
pocket (5 dBm, 5K). (b) Dispersion relation of the nuclear spin wave. The
diagram shows a process for generating the finite wavenumber nuclear spin
wave. (c) Calculated joint density of state. The arrow indicates the absence
of JDoS.

Several nonlinear responses of hybridized nuclear spin excitation mode in
MnCO3 have already been reported [17, 35, 37]. One possible explanation is an
increase of the nuclear spin temperature $T_{\mathrm{n}}$ from the lattice and
electron temperature $T$ [17, 35]. In this case, the absorption shifts toward
a higher frequency. However, the absorptions in Fig. 4(a) keep the low-
frequency edge at $\omega_{n}(k=0,\ T)$, meaning that $T_{\mathrm{n}}\approx
T$. Another possibility is Duffing (Kerr) nonlinearity of the uniform mode
($k=0$) [37]. In this case, the absorption shows a single triangular peak in
the frequency sweep. However, the absence of the absorption (zero-absorption
pocket) in the middle-frequency region [indicated by the arrow in Fig. 4(a) (5
dBm, 5 K)] contradicts with this single-peak scenario.

We attribute the observed nonlinear absorption to the scattering into finite
wavenumber nuclear spin waves. A process to create finite wavenumber nuclear
magnons is allowed through virtual electron magnon scatterings, described in
Fig. 4(b). The rate of the process is proportional to the joint density of
state (JDoS) of the hybridized nuclear spin wave dispersion,

$J(\omega)=\sum_{\mathbf{k}}\delta(\omega_{\mathrm{n}}(k)-\omega).$ (7)

Figure 4(c) shows calculated JDoS, having a low-frequency limit at
$\omega_{\mathrm{n}}(k=0)$, a spread dense area at high frequencies, and an
upper limit at $\omega_{\mathrm{n0}}$. The location of the JDoS absence
indicated by the arrow in Fig. 4(c) matches well with that of the zero-
absorption pocket, being consistent with the scenario shown in Fig. 4(b).

In conclusion, we measured the field dispersion of the hybridized nuclear and
electron spin excitation in MnCO3 by using the broadband microwave
spectroscopy technique. We quantified the internal hyperfine coupling fields,
the cooperativity, the critical exponent, and the Néel temperature from the
observed dispersions. The nonlinearity in the hybridized nuclear spin
excitation mode implies the generation of the finite wavenumber nuclear spin
waves with relatively small microwave fields. The high cooperativity in MnCO3
may be useful for applications such as electric control of nuclear spin
dynamics.

###### Acknowledgements.

We thank S. Daimon, T. Sugimoto, and T. Hioki for discussions and technical
support. This work was partially supported by JSPS KAKENHI (Nos. JP19H05600,
JP20H02599, JP20K15160, JP22H05114, JP23KJ0678, and JP24K01326), JST CREST
(Nos. JPMJCR20C1 and JPMJCR20T2), MEXT Initiative to Establish Next-generation
Novel Integrated Circuits Centers (X-NICS) (No. JPJ011438), Murata Science
Foundation, and Institute for AI and Beyond of the University of Tokyo.

## References

  * Abragam [1961] A. Abragam,  _The principles of nuclear magnetism_ , 32 (Oxford university press, 1961).
  * Shiomi  _et al._ [2019] Y. Shiomi, J. Lustikova, S. Watanabe, D. Hirobe, S. Takahashi, and E. Saitoh, Nature Physics 15, 22 (2019).
  * Kikkawa  _et al._ [2021] T. Kikkawa, D. Reitz, H. Ito, T. Makiuchi, T. Sugimoto, K. Tsunekawa, S. Daimon, K. Oyanagi, R. Ramos, S. Takahashi, Y. Shiomi, Y. Tserkovnyak, and E. Saitoh, Nature Communications 12, 4356 (2021).
  * Rezende [2022] S. M. Rezende, Journal of Applied Physics 132, 091101 (2022).
  * Maekawa  _et al._ [2023] S. Maekawa, T. Kikkawa, H. Chudo, J. Ieda, and E. Saitoh, Journal of Applied Physics 133, 020902 (2023).
  * Kikkawa and Saitoh [2023] T. Kikkawa and E. Saitoh, Annual Review of Condensed Matter Physics 14, 129 (2023).
  * Chigusa  _et al._ [2023] S. Chigusa, T. Moroi, K. Nakayama, and T. Sichanugrist, Phys. Rev. D 108, 095007 (2023).
  * Suhl [1958] H. Suhl, Phys. Rev. 109, 606 (1958).
  * Nakamura [1958] T. Nakamura, Prog. Theor. Phys. 20, 542 (1958).
  * de Gennes  _et al._ [1963] P. G. de Gennes, P. A. Pincus, F. Hartmann-Boutron, and J. M. Winter, Phys. Rev. 129, 1105 (1963).
  * Andrienko  _et al._ [1991] A. V. Andrienko, V. I. Ozhogin, V. L. Safonov, and A. Y. Yakubovskii, Sov. Phys. Usp. 34, 843 (1991).
  * Sherman [2009] D. M. Sherman, American Mineralogist 94, 166 (2009).
  * Effenberger  _et al._ [1981] H. Effenberger, K. Mereiter, and J. Zemann, Zeitschrift für Kristallographie-Crystalline Materials 156, 233 (1981).
  * Lee  _et al._ [2012] J. B. Lee, W. G. Hong, H. J. Kim, Z. Jagličić, S. Jazbec, M. Wencka, A. Jelen, and J. Dolinšek, Physical Review B 86, 224407 (2012).
  * Dzyaloshinsky [1958] I. Dzyaloshinsky, Journal of Physics and Chemistry of Solids 4, 241 (1958).
  * Moriya [1960] T. Moriya, Physical Review 120, 91 (1960).
  * Fink and Shaltiel [1964] H. Fink and D. Shaltiel, Phys. Rev. 136, 218 (1964).
  * Rezende [2020] S. M. Rezende,  _Fundamentals of Magnonics_ , 1st ed. (Springer, 2020).
  * Mims  _et al._ [1967] W. B. Mims, G. E. Devlin, S. Geschwind, and V. Jaccarino, Phyisics Letters 24, 481 (1967).
  * Beeman [1966] D. E. Beeman, Journal of Applied Physics 37, 1136 (1966).
  * Beeman  _et al._ [1966] D. E. Beeman, H. J. Fink, and D. Shaltiel, Physical Review 147, 454 (1966).
  * Zhang  _et al._ [2014] X. Zhang, C.-L. Zou, L. Jiang, and H. X. Tang, Phys. Rev. Lett. 113, 156401 (2014).
  * MacNeill  _et al._ [2019] D. MacNeill, J. T. Hou, D. R. Klein, P. Zhang, P. Jarillo-Herrero, and L. Liu, Physical Review Letters 123, 047204 (2019).
  * Liensberger  _et al._ [2019] L. Liensberger, A. Kamra, H. Maier-Flaig, S. Geprägs, A. Erb, S. T. Goennenwein, R. Gross, W. Belzig, H. Huebl, and M. Weiler, Physical Review Letters 123, 117204 (2019).
  * Hioki  _et al._ [2022] T. Hioki, Y. Hashimoto, and E. Saitoh, Communications Physics 5, 115 (2022).
  * Huebl  _et al._ [2013] H. Huebl, C. W. Zollitsch, J. Lotze, F. Hocke, M. Greifenstein, A. Marx, R. Gross, and S. T. Goennenwein, Physical Review Letters 111, 127003 (2013).
  * Tabuchi  _et al._ [2014] Y. Tabuchi, S. Ishino, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura, Physical Review Letters 113, 083603 (2014).
  * Bourhill  _et al._ [2016] J. Bourhill, N. Kostylev, M. Goryachev, D. Creedon, and M. Tobar, Physical Review B 93, 144420 (2016).
  * Lachance-Quirion  _et al._ [2019] D. Lachance-Quirion, Y. Tabuchi, A. Gloppe, K. Usami, and Y. Nakamura, Applied Physics Express 12, 070101 (2019).
  * Li  _et al._ [2019] Y. Li, T. Polakovic, Y.-L. Wang, J. Xu, S. Lendinez, Z. Zhang, J. Ding, T. Khaire, H. Saglam, R. Divan, _et al._ , Physical Review Letters 123, 107701 (2019).
  * Hou and Liu [2019] J. T. Hou and L. Liu, Physical Review Letters 123, 107702 (2019).
  * Meijer  _et al._ [1970] H. Meijer, L. Pimmelaar, S. Brouwer, and J. Van den Handel, Physica 46, 279 (1970).
  * Pincini  _et al._ [2018] D. Pincini, F. Fabrizi, G. Beutier, G. Nisbet, H. Elnaggar, V. Dmitrienko, M. Katsnelson, Y. Kvashnin, A. Lichtenstein, V. Mazurenko, _et al._ , Physical Review B 98, 104424 (2018).
  * Turov and Kuleev [1966] E. Turov and V. Kuleev, SOVIET PHYSICS JETP 22, 176 (1966).
  * Shaltiel [1966] D. Shaltiel, Physical Review 142, 300 (1966).
  * Abdurakhimov  _et al._ [2015] L. V. Abdurakhimov, Y. M. Bunkov, and D. Konstantinov, Phys. Rev. Lett. 114, 226402 (2015).
  * Abdurakhimov  _et al._ [2018] L. V. Abdurakhimov, M. A. Borish, Y. M. Bunkov, R. R. Gazizulin, D. Konstantinov, M. I. Kurkin, and A. P. Tankeyev, Phys. Rev. B 97, 0244425 (2018).