# Time-local unraveling of non-Markovian stochastic Schrödinger equations

Antoine Tilloy<EMAIL_ADDRESS>Max-Planck-Institut für
Quantenoptik, Hans-Kopfermann-Straße 1, 85748 Garching, Germany

###### Abstract

Non-Markovian stochastic Schrödinger equations (NMSSE) are important tools in
quantum mechanics, from the theory of open systems to foundations. Yet, in
general, they are but formal objects: their solution can be computed
numerically only in some specific cases or perturbatively. This article is
focused on the NMSSE themselves rather than on the open-system evolution they
unravel and aims at making them less abstract. Namely, we propose to write the
stochastic realizations of linear NMSSE as averages over the solutions of an
auxiliary equation with an additional random field. Our method yields a non-
perturbative numerical simulation algorithm for generic linear NMSSE that can
be made arbitrarily accurate for reasonably short times. For isotropic complex
noises, the method extends from linear to non-linear NMSSE and allows to
sample the solutions of norm-preserving NMSSE directly.

##  1 Introduction

Stochastic pure state representations have a wide range of applications in
quantum mechanics: they serve as computational tools to unravel open-system
evolutions [4, 19, 32], as modelling tools to describe continuous measurement
situations [22, 36] and as foundational tools to solve the measurement problem
in models where the superposition principle breaks down [2, 3]. In the
Markovian limit, the solutions of stochastic Schrödinger equations can be
efficiently computed numerically, justifying their wide use in the three
aforementioned fields. In the non-Markovian case however, there exists no
general purpose method to compute a realization of the random pure state
process: plotting a single trajectory is in general impossible. In this
article, we propose to write the solutions of general stochastic Schrödinger
equations as averages over the solutions of time-local stochastic Schrödinger
equations with an auxiliary noise. We shall use the same trick that usually
allows to replace openness by stochasticity, this time to get rid of non-
Markovianity.

Although stochastic Schrödinger equations have a long history dating back to
the eighties [18, 8], the field took off when their potential in numerics was
understood by Dalibard, Castin and Mølmer [4, 26] (see also Dum, Zoller and
Ritsch [12]) in the jump case and by Gisin and Percival [19] in the diffusive
case. The connection with actual quantum trajectories coming from realistic
measurement setups was understood roughly at the same time by Milburn and
Wiseman [35] (see e.g. [22] for an introduction). Interestingly, the
introduction of continuous stochastic pure state equations in foundations is
slightly anterior with the Continuous Spontaneous Localization (CSL) model of
Pearle, Ghirardi and Rimini [27, 17] and the gravity related collapse model of
Diósi [6], both aimed at solving the measurement problem (see [2, 3] for a
review). The generalization of the formalism to the non-Markovian realm was
carried out by Diósi and Strunz [30, 11]. In that case, the measurement
interpretation exists but is admittedly more subtle: a given trajectory only
has a local point by point measurement interpretation [15] but no real time
meaning [9, 10, 36], at least within orthodox interpretations of the formalism
[16]. Non-Markovian stochastic approaches have also infected foundations with
various extensions of the initial CSL proposal [1, 13]. In all these
applications, the practical problem is always the same: the equations are very
formal, involving functional derivatives, and plotting their solution is
typically as hard as solving a general non-Markovian open-system evolution.
There are, of course, quite a few cases where the equations can be solved [31,
23, 37, 24, 25] as well as perturbative methods [5] to deal with a broader
class of problems, but general non-perturbative approaches have so far been
lacking.

In the following, we will introduce stochastic Schrödinger equations as
unravelings of non-Markovian open-system dynamics making an extensive use of
the general framework introduced recently by Diósi and Ferialdi [7]. Yet, our
objective will not primarily be to find a numerically efficient way to compute
the open-system evolution: we will be interested in the stochastic pure state
trajectories themselves, having in mind their broader range of applications.

##  2 General framework

###  2.1 Gaussian master equations

We consider an open quantum system evolving according to a general Gaussian
Master Equation (GME), one of the broadest non-Markovian generalizations of
the Lindblad equation that is analytically tractable. We suppose the system
has a proper Hamiltonian $H_{0}$ and will use the interaction picture for all
operators $\hat{\mathcal{O}}(t)=e^{iH_{0}t}\hat{\mathcal{O}}e^{-iH_{0}t}$. The
GME for the density matrix in interaction picture can then be written as a
time ordered exponential [7]:

$\displaystyle\rho(t)=$ $\displaystyle\Phi_{t}\cdot\rho(0)$ $\displaystyle=$
$\displaystyle\mathcal{T}\exp\bigg{\\{}\int_{0}^{t}\\!\\!\int_{0}^{t}\\!\\!\mathrm{d}\tau\,\mathrm{d}s\,D_{ij}(\tau,s)\Big{[}A^{j}_{L}(s)A^{i}_{R}(\tau)$
(1) $\displaystyle-\theta_{\tau
s}A^{i}_{L}(\tau)A^{j}_{L}(s)-\theta_{s\tau}A^{j}_{R}(s)A^{i}_{R}(\tau)\Big{]}\bigg{\\}}\,\rho(0)$

where $\mathcal{T}$ is the time ordering operator, $\\{\hat{A}^{k}\\}_{1\leq
k\leq n}$ are arbitrary system Hermitian operators, $D_{ij}(\tau,s)$ is a
complex positive semi-definite kernel, $\theta_{\tau s}=\theta(\tau-s)$ is the
Heaviside function, and we have used summation on repeated indices as well as
the left-right notation for superoperators:

$\begin{array}[]{ll}A_{L}\cdot\rho=\hat{A}\rho\;;&A_{R}\cdot\rho=\rho\hat{A}\end{array}.$

Equation (2.1) is a direct non-Markovian generalization of the Lindblad master
equation and reduces to it in the limit $D_{ij}(\tau,s)\rightarrow
d_{ij}\,\delta(\tau-s)$. It is also an equivalent operator rewriting of the
well known quadratic Feynman-Vernon influence functional for systems linearly
coupled to harmonic oscillators [14], a representation used e.g. by Hu, Paz
and Zhang [20, 21] for the Quantum Brownian Motion and which can be derived
from a wide class of microscopic models [34].

At least if $D$ is time translation invariant, equation (2.1) can be obtained
from the linear coupling of the system of interest with a general bosonic bath
(without the Born-Markov approximation). To avoid being excessively abstract,
let us recall the explicit derivation of [7]. We consider a bath constituted
of a continuum of frequency of $n$ types of harmonic oscillators. We write the
corresponding creation and annihilation operators $a_{k}(\omega)$,
$a_{k}^{\dagger}(\omega)$ which verify
$[a_{i}(\omega),a_{j}^{\dagger}(\omega^{\prime})]=\delta_{ij}\delta(\omega-\omega^{\prime})$.
We then couple our system of interest with the bath through the interaction
Hamiltonian:

$H_{\text{int}}=\hat{A}^{k}\otimes\int\mathrm{d}\omega\,\kappa_{k}^{\ell}a_{\ell}(\omega)+\kappa_{k}^{\ell*}a^{\dagger}_{\ell}(\omega),$
(2)

where $\kappa_{k}^{\ell}(\omega)$ is an arbitrary matrix of complex
coefficients. If the initial state of the total system is a product with the
bath ground state, the reduced system density matrix can be expressed exactly
using Wick’s theorem [7] and is given by formula (2.1). In that case, $D$ is
given by:

$D_{ij}(\tau,s)=\int\mathrm{d}\omega\,\kappa_{i}^{k}(\omega)\kappa_{j}^{\ell*}(\omega)e^{-i\omega(\tau-s)}.$
(3)

Here, we are not interested in the open system dynamics per se and thus step
back from its implementation details to take the very general equation (2.1)
as our starting point.

###  2.2 Stochastic unraveling

A stochastic unraveling of equation (2.1) is a set of pure state trajectories
$|\psi_{\xi}(t)\rangle$ indexed by a set of random fields
$\xi=\\{\xi_{k}\\}_{1\leq k\leq n}$ such that:

$\rho(t)=\mathds{E}_{\xi}\big{[}\,|\psi_{\xi}(t)\rangle\langle\psi_{\xi}(t)|\,\big{]},$
(4)

where $\mathds{E}_{\xi}[\cdot]=\int\cdot\,\mathrm{d}\mu_{\circ}(\xi)$ denotes
averaging over the set of complex stochastic fields $\xi$. The interest of
this rewriting is that it decouples the left and right parts of the time
ordered exponential $\Phi_{t}$ of equation (2.1), a trick which is sometimes
called a Hubbard-Stratonovich transformation. Let us now find an explicit
candidate for $|\psi_{\xi}(t)\rangle$. Following e.g. [7], we first introduce
a set of Gaussian complex fields of zero average and thus fully characterized
by its following two point correlation functions:

$\begin{split}\mathds{E}_{\xi}\left[\xi_{i}(\tau)\xi^{*}_{j}(s)\right]&=C_{ij}(\tau,s)\\\
\mathds{E}_{\xi}\big{[}\xi_{i}(\tau)\xi_{j}\,(s)\big{]}&=S_{ij}(\tau,s)\end{split}$
(5)

where $S$ is the _relation_ function, a symmetric complex kernel which can be
chosen freely provided the full correlation kernel is positive semi-definite.

For two sets of fields $a_{t}^{k}$ and $b_{t}^{k}$, the following generalized
characteristic function is obtained by Gaussian integration:

$\mathds{E}_{\xi}\left[\exp\left\\{-i\\!\\!\int_{0}^{t}\\!\\!\mathrm{d}s(a^{k}_{s}\xi_{k}(s)-b^{k}_{s}\xi_{k}^{*}(s))\right\\}\right]=\exp\left\\{\int_{0}^{t}\\!\\!\int_{0}^{t}\\!\\!\mathrm{d}\tau\mathrm{d}s\,C_{ij}(\tau,s)a^{i}_{\tau}b^{j}_{s}-\frac{S_{ij}(\tau,s)a^{i}_{\tau}a^{j}_{s}+S^{*}_{ij}(\tau,s)b^{i}_{\tau}b^{j}_{s}}{2}\right\\},$
(6)

an equality that can be further generalized for $a$’s and $b$’s promoted to
operators provided both sides are time ordered:

$\mathds{E}_{\xi}\left[\mathcal{T}\exp\left\\{-i\\!\\!\int_{0}^{t}\\!\\!\mathrm{d}s(\hat{a}^{k}_{s}\xi_{k}(s)-\hat{b}^{k}_{s}\xi_{k}^{*}(s))\right\\}\right]=\mathcal{T}\exp\left\\{\int_{0}^{t}\\!\\!\int_{0}^{t}\\!\\!\mathrm{d}\tau\mathrm{d}s\,C_{ij}(\tau,s)\hat{a}^{i}_{\tau}\hat{b}^{j}_{s}-\frac{S_{ij}(\tau,s)\hat{a}^{i}_{\tau}\hat{a}^{j}_{s}+S^{*}_{ij}(\tau,s)\hat{b}^{i}_{\tau}\hat{b}^{j}_{s}}{2}\right\\}.$
(7)

The similarity between the r.h.s of (7) and (2.1) suggests to identify
$\mathds{E}_{\xi}[|\psi_{\xi}\rangle\langle\psi_{\xi}|]$ with the l.h.s of (7)
with $\hat{a}^{k}_{s}=\hat{A}_{L}^{k}(s)$,
$\hat{b}^{k}_{s}=\hat{A}^{k}_{R}(s)$ and $C_{ij}(\tau,s)=D_{ij}(\tau,s)$.
Actually, this does not fully do the trick, a counter term remains, and one
sees following [7] that the correct prescription is:

$\displaystyle|\psi_{\xi}(t)\rangle=$
$\displaystyle\hat{G}_{\xi}(t)|\psi_{0}\rangle$ $\displaystyle:=$
$\displaystyle\mathcal{T}\exp\bigg{\\{}-i\int_{0}^{t}\\!\\!\mathrm{d}s\hat{A}^{k}(s)\xi_{k}(s)$
(8)
$\displaystyle\\!\\!\\!\\!-\\!\\!\int_{0}^{t}\\!\\!\int_{0}^{t}\\!\\!\mathrm{d}\tau\,\mathrm{d}s\,\theta_{\tau
s}\,[D-S]_{ij}(\tau,s)\hat{A}^{i}(\tau)\hat{A}^{j}(s)\bigg{\\}}|\psi_{0}\rangle,$

which fulfills the unraveling condition (4). Equation (2.2) can subsequently
be written in differential form to yield a non-Markovian stochastic
Schrödinger equation:

$\begin{split}\frac{\mathrm{d}}{\mathrm{d}t}|\psi_{\xi}(t)\rangle=&-i\hat{A}^{i}(t)\bigg{[}\xi_{i}(t)\\\
&+\int_{0}^{t}\\!\\!\mathrm{d}s\,[D-S]_{ij}(t,s)\frac{\delta}{\delta\xi_{j}(s)}\bigg{]}|\psi_{\xi}(t)\rangle.\end{split}$
(9)

To our knowledge, this equation is the most general _linear_ non-Markovian
stochastic Schrödinger equation ever proposed [7]. Unfortunately, equation (9)
is quite formal unless it is possible to write the functional derivative111The
reader puzzled by the rigorous definition of the functional derivative can as
well use (2.2) in all situations. as a simple local operator acting on
$|\psi_{\xi}(t)\rangle$.

##  3 Main result

The core objective of this article is to show how the solutions of (9) can
nonetheless be computed in the general case by reusing a stochastic unraveling
trick, i.e. by writing this time:

$|\psi_{\xi}(t)\rangle=\mathds{E}_{\eta}\big{[}\,|\psi_{\xi,\eta}(t)\rangle\,\big{]}$
(10)

where $\eta$ is an auxiliary set of classical Gaussian fields. More precisely,
we can try to write the evolution operator of (2.2) in the following way:

$\hat{G}_{\xi}(t)=\mathds{E}_{\eta}\left[\mathcal{T}\exp\left\\{-i\\!\int_{0}^{t}\\!\\!\mathrm{d}s\hat{A}^{k}(s)[\xi_{k}(s)+\eta_{k}(s)]\right\\}\right].$
(11)

Using again the Gaussian integration formula (7), one sees that the new
unraveling condition (10) is satisfied provided $\eta$ is a set of Gaussian
fields with zero mean and two-point correlation functions given by:

$\begin{split}\mathds{E}_{\eta}[\eta_{i}(\tau)\,\eta_{j}(s)]&=K_{ij}(\tau,s)\\\
\mathds{E}_{\eta}[\eta_{i}(\tau)\,\eta^{*}_{j}(s)]&=J_{ij}(\tau,s)\end{split}$
(12)

with:

$K_{ij}(\tau,s)=\theta_{\tau
s}\,[D-S]_{ij}(\tau,s)+\theta_{s\tau}\,[D-S]_{ji}(s,\tau)$ (13)

and where the other correlation function $J$ can be chosen freely provided the
total kernel $\Gamma$:

$\Gamma=\left(\begin{array}[]{cc}J&K\\\ K^{*}&J^{*}\end{array}\right)$ (14)

is positive semi-definite222For a kernel $K$ given by equation (13), there
always exists a kernel $J$ such that $\Gamma$ is positive semi-definite. The
existence of a random field $\eta$ verifying (12) is then a consequence of the
Bochner-Minlos theorem..

Equation (11) can be written in an explicit linear differential form free of
functional derivatives:

$\frac{\mathrm{d}}{\mathrm{d}t}|\psi_{\xi,\eta}(t)\rangle=-i\hat{A}^{k}(t)[\xi_{k}(t)+\eta_{k}(t)]\,|\psi_{\xi,\eta}(t)\rangle.$
(15)

This latter equation, together with its unraveling interpretation (10) is the
main result of this article. The evolution given by equation (15) is _time-
local_ in the sense that once the random fields $\xi$ and $\eta$ are fixed,
the state can be evolved without reference to the past. Although our interest
was primarily in the stochastic pure state $|\psi_{\xi}\rangle$ itself, we can
write the open-system density matrix as a double average over our successive
unravelings (4) and (10):

$\rho(t)=\mathds{E}_{\xi}\Big{[}\,\mathds{E}_{\eta}\big{[}\,|\psi_{\xi,\eta}(t)\rangle\,\big{]}\;\,\mathds{E}_{\eta}\big{[}\,\langle\psi_{\xi,\eta}(t)|\,\big{]}\,\Big{]}.$
(16)

This can be written equivalently:

$\rho(t)=\mathds{E}_{\xi,\eta^{(1)},\eta^{(2)}}\big{[}\,|\psi_{\xi,\eta^{(1)}}(t)\rangle\langle\psi_{\xi,\eta^{(2)}}(t)|\,\big{]},$
(17)

where $\eta^{(1)}$ and $\eta^{(2)}$ are independent. This is the two state
unraveling proposal of Stockburger and Grabert [29, 28]. We have thus found a
connection between the standard stochastic pure state representations and the
more exotic two state vector method used in numerical approaches to open-
systems. This allows us to identify at least part of the freedom in the noise
present in this latter method as coming from the kernel $S$ which encodes
different stochastic Schrödinger evolutions corresponding to the same open-
system evolution. We can further risk the following heuristic interpretation
of the two noises. The noise $\xi$ is classical in the sense that it is
averaged over at the density matrix level. The noise $\eta$, on the other
hand, is quantum or coherent as all the possible contributions are summed over
coherently at the pure state level.

##  4 Norm preserving unravelings

###  4.1 Context and objectives

The linear stochastic differential equation (9) does not preserve the norm of
the state vector $|\psi_{\xi}\rangle$. Normalized pure states are usually
preferred in foundations but also for numerics: the norm of
$|\psi_{\xi}\rangle$ typically diffuses and for large times, most trajectories
give a vanishing contribution to the average (4) while nearly all the weight
is provided by rare events.

To obtain a norm preserving evolution from the linear one, one needs to define
a normalized state $|\widetilde{\psi}_{\xi}\rangle$:

$|\widetilde{\psi}_{\xi}(t)\rangle=\frac{|\psi_{\xi}(t)\rangle}{\sqrt{\langle\psi_{\xi}(t)|\psi_{\xi}(t)\rangle}}$
(18)

and also to introduce a new “cooked” probability measure:

$\mathrm{d}\mu_{t}(\xi)=\langle\psi_{\xi}(t)|\psi_{\xi}(t)\rangle\,\mathrm{d}\mu_{\circ}(\xi)$
(19)

which insures that the unraveling condition
$\rho(t)=\mathds{E}^{t}_{\xi}\left[\,|\widetilde{\psi}_{\xi}(t)\rangle\langle\widetilde{\psi}_{\xi}(t)|\right]$
–where $\mathds{E}^{t}_{\xi}$ is the expectation value taken with the new
measure $\mu_{t}$– is still valid. It is important to note at that stage that
the state $|\widetilde{\psi}_{\xi}(t)\rangle$ with the probability measure
$\mathrm{d}\mu_{t}(\xi)$ unravels the open-system evolution only for a single
instant $t$. If one wants a trajectory that unravels the open-system evolution
at all times, it is necessary to dynamically redefine the field variable
$\xi$. It can be done by introducing $\xi^{[t]}$, a _complete_ random field
trajectory depending on $\xi$ such that for any functional $f$:

$\mathds{E}^{t}_{\xi}[f(\xi)]=\mathds{E}_{\xi}\Big{[}f\big{(}\xi^{[t]}(\xi)\big{)}\Big{]}.$
(20)

This way, starting from a field realization $\xi$ drawn from the Gaussian
measure $\mathrm{d}\mu_{\circ}$, one can compute a trajectory
$t\mapsto|\widetilde{\psi}_{\xi^{[t]}}(t)\rangle$ that unravels the open
evolution for _all_ times. By “computing the norm preserving unraveling” one
may thus want to sample two different things:

  1. 1.

the state $|\widetilde{\psi}_{\xi^{[t]}}(t)\rangle$ for a single time,

  2. 2.

the full trajectory $t\mapsto|\widetilde{\psi}_{\xi^{[t]}}(t)\rangle$.

While knowing the state at a single time may be enough in most cases, it is
not sufficient to compute correlation functions of the stochastic states at
different times (something one might want to do for collapse models in
foundations). In any case, the ideal solution would be to find a way to write
an explicit stochastic differential equation sampling the states of interest
directly. This is unfortunately not a trivial endeavour and the method we
shall propose is arguably less appealing than in the linear case.

###  4.2 Computation of the full trajectory for S=0

When $S=0$, there exists an efficient way to construct the full norm
preserving trajectory $t\mapsto|\widetilde{\psi}_{\xi^{[t]}}(t)\rangle$,
exploiting the fact that $\xi$ and $\xi^{[t]}$, defined implicitly in (20),
can be connected explicitly in that case. Indeed, after relatively painful
computations (see [19, 35, 33]) one can show that for $S=0$:

$\xi^{[t]}_{k}(u)=\xi_{k}(u)+i\int_{0}^{t}\mathrm{d}s\,D_{k\ell}(u,s)\langle
A^{\ell}(s)\rangle_{s},$ (21)

with $\langle
A_{\ell}(s)\rangle_{s}=\langle\widetilde{\psi}_{\xi^{[s]}}(s)|A^{\ell}(s)|\widetilde{\psi}_{\xi^{[s]}}(s)\rangle$.
Using this result, it is possible to compute
$|\widetilde{\psi}_{\xi^{[t]}}(t)\rangle$ inductively.

We start from $\xi^{[0]}=\xi$ sampled from the Gaussian measure
$\mathrm{d}\mu_{\circ}$. At a time $v\leq t$ we assume that:

  * –

$\forall u\leq v$, $|\widetilde{\psi}_{\xi^{[u]}}(u)\rangle$ is known,

  * –

$\forall u\leq t$, $\xi^{[v]}_{u}$ is known.

We first compute $\xi^{[v+\mathrm{d}v]}$ for $u\leq t$ using (21):

$\xi^{[v+\mathrm{d}v]}_{k}(u)=\xi^{[v]}_{k}(u)+iD_{k,\ell}(u,v)\langle
A^{\ell}(v)\rangle_{v}\mathrm{d}v+O(\mathrm{d}v^{2}),$ (22)

We then compute
$|\widetilde{\psi}_{\xi^{[v+\mathrm{d}v]}}(v+\mathrm{d}v)\rangle$ solving the
linear equation (9) up to time $v+\mathrm{d}v$ with the new complete noise
trajectory $\xi^{[v+\mathrm{d}v]}$ and normalizing the result. This is done
using the unraveling (15). By induction, we thus obtain
$|\widetilde{\psi}_{\xi^{[t]}}(t)\rangle$ for all $t$ starting from a single
realization of the noise $\xi$ (see figure 1). For that matter, it is
necessary to solve the linear stochastic Schrödinger equation from the
beginning with a new field at every time step. Hence this adds a linear a
overhead $\propto t/\mathrm{d}t$ to the time needed to compute linear
trajectories. As the errors coming from the discretization are typically much
smaller than the statistical errors coming from the averaging of (15), there
is however no loss in precision.

Figure 1: Schematic representation of the procedure described in 4.2. The
states and fields are computed from bottom to top for increasing values of
$v$. The final product, the norm-preserving trajectory, is obtained for $u=v$.

Figure 2: Left: Snapshot of a linear trajectory $|\psi_{\xi}\rangle$ for a
random realization of $\xi$. The exact trajectory is obtained from equation
(26) while the approximate trajectory is obtained from (15) –integrated with
$\mathrm{d}t=10^{-2}$–, and the estimation
$|\psi_{\xi}\rangle\simeq\frac{1}{N}\sum_{n=1}^{N}|\psi_{\xi,\eta}\rangle$
with $N=10^{3}$ independent realizations of $\eta$ and the same fixed $\xi$.
Right: Snapshot of a non-linear trajectory
$|\widetilde{\psi}_{\xi^{[t]}}(t)\rangle$ for the same realization $\xi$, and
computed with the method of 4.2. In the algorithm, the linear NMSSE is solved
by integrating (15) with $\mathrm{d}t=10^{-2}$ and averaging over $N=10^{3}$
independent realizations of $\eta$.
The shaded areas indicate the zones where $50\%$ ($\simeq 0.5\sigma$) and
$85\%$ ($\simeq 1\sigma$) of the approximate trajectories would lie (obtained
by computing $100$ samples for the same $\xi$). In both cases, number of
realizations is voluntarily left small to make the estimate of the error
visible as the latter decreases $\propto N^{-1/2}$.

###  4.3 General case, single time

In the general case, the previous method cannot be used as we do not know of a
generalization of formula (20) for $S\neq 0$. Finding a way to sample from
norm-preserving trajectories has eluded us. In most cases however, one is only
interested in the properties of the stochastic state at a single time. In such
a situation, the linear unraveling is sufficient to compute everything.
Indeed, let us consider an arbitrary function $\varphi$ of the norm-preserving
trajectory $|\widetilde{\psi}_{\xi^{[t]}}(t)\rangle$ at a _fixed_ time $t$. We
have:

$\displaystyle\mathds{E}$
$\displaystyle\big{[}\varphi(|\widetilde{\psi}_{\xi^{[t]}}(t)\rangle)\big{]}=\mathds{E}^{t}_{\xi}\big{[}\varphi(|\widetilde{\psi}_{\xi}(t)\rangle)\big{]}$
(23)
$\displaystyle=\mathds{E}^{t}_{\xi}\left[\varphi\left(\frac{|\psi_{\xi}(t)\rangle}{\sqrt{\langle\psi_{\xi}(t)|\psi_{\xi}(t)\rangle}}\right)\right]$
(24)
$\displaystyle=\mathds{E}_{\xi}\left[\varphi\left(\frac{|\psi_{\xi}(t)\rangle}{\sqrt{\langle\psi_{\xi}(t)|\psi_{\xi}(t)\rangle}}\right)\langle\psi_{\xi}(t)|\psi_{\xi}(t)\rangle\right].$
(25)

That is, expectation values of all time-local functionals of the norm
preserving trajectory can be computed by averaging over the solutions of the
linear NMSSE.

##  5 Example

###  5.1 Linear evolution

We now illustrate our method on an example where the stochastic Schrödinger
equation is fully explicit. For that matter, we consider a single Hermitian
operator $\hat{A}$ commuting with $H_{0}$ (s.t. $\hat{A}(t)=\hat{A}$) and
purely isotropic noise, i.e. $S=0$, which gives the following stochastic
Schrödinger equation:

$\frac{\mathrm{d}}{\mathrm{d}t}|\psi_{\xi}(t)\rangle=-i\hat{A}\left[\xi(t)+\int_{0}^{t}\\!\\!\mathrm{d}s\,D(t,s)\frac{\delta}{\delta\xi(s)}\right]\\!|\psi_{\xi}(t)\rangle.$
(26)

For a given realization of $\xi$, $|\psi_{\xi}\rangle$ obeys an explicit
differential equation. Indeed, using the integral form (2.2), we can replace
the functional derivative:

$\frac{\delta}{\delta\xi(s)}|\psi_{\xi}(t)\rangle=-i\hat{A}|\psi_{\xi}(t)\rangle.$
(27)

This example thus offers a nice test bed for our stochastic approach. In this
example, the relaxation timescale responsible for the non-Markovianity and the
timescale of the system-bath coupling are both of order $1$ and we would thus
be deep in the non-Markovian regime had the system Hamiltonian not been
trivial. To compute the linear trajectories, we average over the solutions of
(15) with $\eta$ a real noise such that $\mathds{E}[\eta_{t}\eta_{s}]=D(t,s)$.
Numerical results for $\hat{A}=\mathrm{diag}(1,0,-1)$, $D(t,0)=e^{-t}$, and
$|\psi(0)\rangle=(1,1,1)/\sqrt{3}$ are shown in figure 2.

Our method is efficient for evolution times of the order of the the system-
bath coupling timescale but becomes expensive for larger times. Indeed, we are
summing states $|\psi_{\xi,\eta}(t)\rangle$ with a phase with respect to
$|\psi_{\xi}(t)\rangle$ that is more and more random as time grows. This is an
important limit of this scheme for large time.

###  5.2 Norm preserving evolution

We also sample the norm-preserving non-linear trajectories using the method of
section 4.2. The precision for the sampling of a single trajectory
$|\psi_{\xi^{[t]}}(t)\rangle$ is very similar to that of the linear case (see
figure 2).

To test the robustness of the method, we also compute the probability
distribution of
$\langle\psi_{\xi^{[t]}}(t)|\hat{A}|\psi_{\xi^{[t]}}(t)\rangle$ for two
different times using norm-preserving trajectories computed with the numerical
method of 4.2 and using exact linear trajectories and the re-weighting of 4.3.
Up to statistical and discretization errors, the second method is exact. The
results are displayed in figure 3 and show an excellent agreement between the
two methods.

Figure 3: Top: Snapshot of $80$ normalized trajectories computed with the
method of 4.2. In the algorithm, the linear equation is solved by integrating
(15) with $\mathrm{d}t=10^{-2}$ and averaging over $N=10^{3}$ independent
realizations of $\eta$.
Bottom: Empirical probability densities; green histogram computed by sampling
$10^{3}$ non-linear trajectories (each computed with the same method as in the
snapshot), brown line computed by sampling $10^{5}$ _exact_ linear
trajectories and re-weighting a posteriori as in 4.3

##  6 Conclusion

In this article, we have introduced a rewriting (or “unraveling”) of the
solutions of formal NMSSE as averages over the solutions of explicit time-
local stochastic differential equations. This new formulation may provide a
basis for further analytical studies, notably as it gives an extremely simple
new way to define NMSSE. More importantly, it immediately offers a method to
compute numerically their solutions. Although we have proposed a heuristic
interpretation of the quantum and classical character of the two noises
involved, the auxiliary state $|\psi_{\xi,\eta}\rangle$ does not yet have a
clear operational characterization. The latter would certainly help understand
the freedom in the noise kernel $J$. We have also extended our method to non-
linear norm-preserving NMSSE in the case where $S=0$, i.e. where the complex
noise $\xi$ is isotropic. Computing expectation values of time non-local
functionals of norm-preserving trajectories when $S\neq 0$ is still an open
problem. Our work has been focused on the NMSSE as an end in themselves rather
than on the open-system evolution they unravel. It is however possible that
our methods, especially the one to sample non-linear NMSSE directly, might
ultimately help improve our ability to deal with non-Markovian open-systems
numerically.

I thank Lajos Diósi and Luca Ferialdi for helpful discussions. Part of this
work was supported by the Alexander von Humboldt foundation and the Agence
Nationale de la Recherche (ANR) contract ANR-14-CE25-0003-01.

## References

  * Adler and Bassi [2007] S. L. Adler and A. Bassi.  Collapse models with non-white noises.  _Journal of Physics A: Mathematical and Theoretical_ , 40(50):15083, 2007.  doi: 10.1088/1751-8113/40/50/012.
  * Bassi and Ghirardi [2003] A. Bassi and G. Ghirardi.  Dynamical reduction models.  _Physics Reports_ , 379(5–6):257 – 426, 2003\.  doi: 10.1016/S0370-1573(03)00103-0.
  * Bassi et al. [2013] A. Bassi, K. Lochan, S. Satin, T. P. Singh, and H. Ulbricht.  Models of wave-function collapse, underlying theories, and experimental tests.  _Rev. Mod. Phys._ , 85:471–527, Apr 2013.  doi: 10.1103/RevModPhys.85.471.
  * Dalibard et al. [1992] J. Dalibard, Y. Castin, and K. Mølmer.  Wave-function approach to dissipative processes in quantum optics.  _Phys. Rev. Lett._ , 68:580–583, Feb 1992.  doi: 10.1103/PhysRevLett.68.580.
  * de Vega et al. [2005] I. de Vega, D. Alonso, and P. Gaspard.  Two-level system immersed in a photonic band-gap material: A non-markovian stochastic schrödinger-equation approach.  _Phys. Rev. A_ , 71:023812, Feb 2005.  doi: 10.1103/PhysRevA.71.023812.
  * Diósi [1989] L. Diósi.  Models for universal reduction of macroscopic quantum fluctuations.  _Phys. Rev. A_ , 40:1165–1174, Aug 1989.  doi: 10.1103/PhysRevA.40.1165.
  * Diósi and Ferialdi [2014] L. Diósi and L. Ferialdi.  General non-markovian structure of gaussian master and stochastic schrödinger equations.  _Phys. Rev. Lett._ , 113:200403, Nov 2014.  doi: 10.1103/PhysRevLett.113.200403.
  * Diósi [1986] L. Diósi.  Stochastic pure state representation for open quantum systems.  _Phys. Lett. A_ , 114(8):451 – 454, 1986.  doi: 10.1016/0375-9601(86)90692-4.
  * Diósi [2008a] L. Diósi.  Non-markovian continuous quantum measurement of retarded observables.  _Phys. Rev. Lett._ , 100:080401, Feb 2008a.  doi: 10.1103/PhysRevLett.100.080401.
  * Diósi [2008b] L. Diósi.  Erratum: Non-markovian continuous quantum measurement of retarded observables [phys. rev. lett. 100 , 080401 (2008)].  _Phys. Rev. Lett._ , 101:149902, Oct 2008b.  doi: 10.1103/PhysRevLett.101.149902.
  * Diósi and Strunz [1997] L. Diósi and W. T. Strunz.  The non-markovian stochastic schrödinger equation for open systems.  _Phys. Lett. A_ , 235(6):569–573, 1997.  doi: 10.1016/S0375-9601(97)00717-2.
  * Dum et al. [1992] R. Dum, P. Zoller, and H. Ritsch.  Monte carlo simulation of the atomic master equation for spontaneous emission.  _Phys. Rev. A_ , 45:4879–4887, Apr 1992.  doi: 10.1103/PhysRevA.45.4879.
  * Ferialdi and Bassi [2012] L. Ferialdi and A. Bassi.  Dissipative collapse models with nonwhite noises.  _Phys. Rev. A_ , 86:022108, Aug 2012.  doi: 10.1103/PhysRevA.86.022108.
  * Feynman and Vernon [1963] R. P. Feynman and F. L. Vernon.  The theory of a general quantum system interacting with a linear dissipative system.  _Annals of physics_ , 24:118–173, 1963.  doi: 10.1016/0003-4916(63)90068-X.
  * Gambetta and Wiseman [2002] J. Gambetta and H. M. Wiseman.  Non-markovian stochastic schrödinger equations: Generalization to real-valued noise using quantum-measurement theory.  _Phys. Rev. A_ , 66:012108, Jul 2002.  doi: 10.1103/PhysRevA.66.012108.
  * Gambetta and Wiseman [2003] J. Gambetta and H. M. Wiseman.  Interpretation of non-markovian stochastic schrödinger equations as a hidden-variable theory.  _Phys. Rev. A_ , 68:062104, Dec 2003.  doi: 10.1103/PhysRevA.68.062104.
  * Ghirardi et al. [1990] G. C. Ghirardi, P. Pearle, and A. Rimini.  Markov processes in hilbert space and continuous spontaneous localization of systems of identical particles.  _Phys. Rev. A_ , 42:78–89, Jul 1990.  doi: 10.1103/PhysRevA.42.78.
  * Gisin [1984] N. Gisin.  Quantum measurements and stochastic processes.  _Phys. Rev. Lett._ , 52:1657–1660, May 1984.  doi: 10.1103/PhysRevLett.52.1657.
  * Gisin and Percival [1992] N. Gisin and I. C. Percival.  The quantum-state diffusion model applied to open systems.  _J. Phys. A: Math. Gen._ , 25(21):5677, 1992.  doi: 10.1088/0305-4470/25/21/023.
  * Hu et al. [1992] B. L. Hu, J. P. Paz, and Y. Zhang.  Quantum brownian motion in a general environment: Exact master equation with nonlocal dissipation and colored noise.  _Phys. Rev. D_ , 45:2843–2861, Apr 1992.  doi: 10.1103/PhysRevD.45.2843.
  * Hu et al. [1993] B. L. Hu, J. P. Paz, and Y. Zhang.  Quantum brownian motion in a general environment. ii. nonlinear coupling and perturbative approach.  _Phys. Rev. D_ , 47:1576–1594, Feb 1993.  doi: 10.1103/PhysRevD.47.1576.
  * Jacobs and Steck [2006] K. Jacobs and D. A. Steck.  A straightforward introduction to continuous quantum measurement.  _Contemporary Physics_ , 47(5):279–303, 2006\.  doi: 10.1080/00107510601101934.
  * Jing and Yu [2010] J. Jing and T. Yu.  Non-markovian relaxation of a three-level system: Quantum trajectory approach.  _Phys. Rev. Lett._ , 105:240403, Dec 2010.  doi: 10.1103/PhysRevLett.105.240403.
  * Jing et al. [2012] J. Jing, X. Zhao, J. Q. You, and T. Yu.  Time-local quantum-state-diffusion equation for multilevel quantum systems.  _Phys. Rev. A_ , 85:042106, Apr 2012.  doi: 10.1103/PhysRevA.85.042106.
  * Jing et al. [2013] J. Jing, X. Zhao, J. Q. You, W. T. Strunz, and T. Yu.  Many-body quantum trajectories of non-markovian open systems.  _Phys. Rev. A_ , 88:052122, Nov 2013.  doi: 10.1103/PhysRevA.88.052122.
  * Mølmer et al. [1993] K. Mølmer, Y. Castin, and J. Dalibard.  Monte carlo wave-function method in quantum optics.  _J. Opt. Soc. Am. B_ , 10(3):524–538, Mar 1993\.  doi: 10.1364/JOSAB.10.000524.
  * Pearle [1989] P. Pearle.  Combining stochastic dynamical state-vector reduction with spontaneous localization.  _Phys. Rev. A_ , 39:2277–2289, Mar 1989.  doi: 10.1103/PhysRevA.39.2277.
  * Stockburger [2004] J. T. Stockburger.  Simulating spin-boson dynamics with stochastic liouville–von neumann equations.  _Chemical physics_ , 296(2):159–169, 2004.  doi: 10.1016/j.chemphys.2003.09.014.
  * Stockburger and Grabert [2002] J. T. Stockburger and H. Grabert.  Exact $\mathit{c}$-number representation of non-markovian quantum dissipation.  _Phys. Rev. Lett._ , 88:170407, Apr 2002.  doi: 10.1103/PhysRevLett.88.170407.
  * Strunz [1996] W. T. Strunz.  Linear quantum state diffusion for non-markovian open quantum systems.  _Phys. Lett. A_ , 224(1):25 – 30, 1996.  doi: 10.1016/S0375-9601(96)00805-5.
  * Strunz and Yu [2004] W. T. Strunz and T. Yu.  Convolutionless non-markovian master equations and quantum trajectories: Brownian motion.  _Phys. Rev. A_ , 69:052115, May 2004.  doi: 10.1103/PhysRevA.69.052115.
  * Strunz et al. [1999] W. T. Strunz, L. Diósi, N. Gisin, and T. Yu.  Quantum trajectories for brownian motion.  _Phys. Rev. Lett._ , 83:4909–4913, Dec 1999.  doi: 10.1103/PhysRevLett.83.4909.
  * Tilloy [2017] A. Tilloy.  Interacting quantum field theories as relativistic statistical field theories of local beables.  _arXiv:1702.06325_ , 2017.
  * Weiss [1999] U. Weiss.  _Quantum dissipative systems_ , volume 10.  World Scientific, 1999.
  * Wiseman and Milburn [1993] H. M. Wiseman and G. J. Milburn.  Quantum theory of field-quadrature measurements.  _Phys. Rev. A_ , 47:642–662, Jan 1993.  doi: 10.1103/PhysRevA.47.642.
  * Wiseman and Gambetta [2008] H. M. Wiseman and J. M. Gambetta.  Pure-state quantum trajectories for general non-markovian systems do not exist.  _Phys. Rev. Lett._ , 101:140401, Sep 2008.  doi: 10.1103/PhysRevLett.101.140401.
  * Zhao et al. [2011] X. Zhao, J. Jing, B. Corn, and T. Yu.  Dynamics of interacting qubits coupled to a common bath: Non-markovian quantum-state-diffusion approach.  _Phys. Rev. A_ , 84:032101, Sep 2011.  doi: 10.1103/PhysRevA.84.032101.


# Convective, absolute and global azimuthal magnetorotational instabilities

A. Mishra1<EMAIL_ADDRESS>G. Mamatsashvili1, 2, 3  V. Galindo1  F. Stefani1
1Helmholtz-Zentrum Dresden-Rossendorf, Bautzner Landstraße 400, D-01328
Dresden, Germany 2Abastumani Astrophysical Observatory, Abastumani 0301,
Georgia 3Institute of Geophysics, Tbilisi State University, Tbilisi 0193,
Georgia

###### Abstract

We study the convective and absolute forms of azimuthal magnetorotational
instability (AMRI) in a Taylor-Couette (TC) flow with an imposed azimuthal
magnetic field. We show that the domain of the convective AMRI is wider than
that of the absolute AMRI. Actually, it is the absolute instability which is
the most relevant and important for magnetic TC flow experiments. The absolute
AMRI, unlike the convective one, stays in the device, displaying a sustained
growth that can be experimentally detected. We also study the global AMRI in a
TC flow of finite height using DNS and find that its emerging butterfly-type
structure – a spatio-temporal variation in the form of upward and downward
traveling waves – is in a very good agreement with the linear stability
analysis, which indicates the presence of two dominant absolute AMRI modes in
the flow giving rise to this global butterfly pattern.

###### keywords:

Absolute/convective instability, Taylor–Couette flow, MHD and
electrohydrodynamics

##  1 Introduction

Many processes of geophysical and astrophysical interest entail the
interaction of magnetic fields with conducting fluids or plasma. The dynamo
effect is responsible for the self-excitation of magnetic fields in planets,
stars and galaxies (Rincon, 2019). A key player in cosmic structure formation
is the magnetorotational instability (MRI, Balbus & Hawley, 1998), which
triggers outward transport of angular momentum in accretion disks and mass
concentration onto central objects. While dynamo action and MRI are usually
considered as separate effects, they are treated as nonlinearly interwoven and
mutually reinforcing processes in the modern concept of MRI dynamos (Rincon,
2019; Mamatsashvili et al., 2020).

Over the last two decades, a great deal of theoretical and numerical work on
dynamo action and MRI has been complemented by a number of pertinent
experiments (Stefani et al., 2008). The threshold of self-excitation was
achieved in the liquid sodium experiments in Riga (Gailitis et al., 2018),
Karlsruhe (Müller et al., 2008) and Cadarache (Monchaux et al., 2009). The
helical (Stefani et al., 2006) and azimuthal (Seilmayer et al., 2014) variants
of MRI were observed in the medium-size PROMISE experiment, while an
unequivocal proof of the standard MRI (with an imposed purely axial magnetic
field) is still missing, despite great efforts and promising initial results
(Sisan et al., 2004; Nornberg et al., 2010).

Interestingly, those experimental efforts brought to the fore the distinction
between convective, absolute and global instabilities, which are well known in
hydrodynamics and plasma physics communities (e.g., Huerre & Monkewitz, 1990;
Chomaz, 2005), but seldom addressed in geophysical and astrophysical
magnetohydrodynamics (MHD). It was key for the success of the Riga experiment
to transform the convective instability of the paradigmatic Ponomarenko (1973)
dynamo (a conducting rigid rod, moving in a screw-like manner through an
infinitely extended medium with the same conductivity) into an absolute (and a
global) one by adding a concentric outer return flow to the inner motion of
the rod. The corresponding one-dimensional (1D) eigenvalue problem in the
radial direction and the search for saddle points in the complex wavenumber
plane, representing the absolute instability, are described in Gailitis
(1996). In this case, the absolute instability for the finite length system
coincided quite accurately with the global instability that was later obtained
with a two-dimensional code (Stefani et al., 1999). A similar distinction
between convective, absolute and global instabilities was discussed for the
helical MRI (HMRI) experiment PROMISE (Priede & Gerbeth, 2009; Stefani et al.,
2009). In contrast to the Riga dynamo, where this distinction obviously stems
from the axial flow through the device, the axial propagation of the unstable
HMRI mode is tied to the direction of the background Poynting flux (Liu et
al., 2006).

This paper is concerned with a related MRI problem for which the relevance of
distinguishing among convective, absolute and global instabilities is far less
evident. The azimuthal MRI (AMRI, Hollerbach et al., 2010) is a non-
axisymmetric instability, with dominant azimuthal wavenumbers $m=\pm 1$,
arising in a (hydrodynamically stable) differentially rotating flow in the
presence of a purely azimuthal magnetic field. While any weak additional axial
field would naturally result in a preferred direction and preponderance of
either the $m=1$ or $m=-1$ mode, with no axial field these two modes have the
same weight and it might be naively expected that they add up to a standing
wave. Under this assumption, the two inter-penetrating upward and downward
traveling waves, as observed experimentally in Seilmayer et al. (2014), were
interpreted in terms of the breaking of the axial homogeneity by the endcaps
of the cylindrical Taylor-Couette (TC) device of PROMISE. However, we will
show in this paper that the tendency of AMRI waves to travel away from the
mid-height of the cylinders, exhibiting a butterfly-like spatio-temporal
variation (analogous to the famous “butterfly diagram” in solar physics), has
a deeper rooting as it occurs already in a TC flow that is infinitely extended
in the axial direction.

##  2 Mathematical Formulation

The setup is that of an infinitely long TC flow of a conducting fluid in
cylindrical coordinates $(r,\phi,z)$ with an imposed purely azimuthal magnetic
field, which is usually adopted in AMRI studies. The inner and outer cylinders
with radii $r_{i}$ and $r_{o}$, respectively, rotate with angular velocities
$\Omega_{i}$ and $\Omega_{o}$ around $z$-axis, driving in the gap between them
an azimuthal velocity, ${\bf u}_{0}=r\Omega(r){\bf e}_{\phi}$, with the radial
profile of the angular velocity $\Omega(r)=C_{1}+C_{2}/r^{2}$, where
$C_{1}=(r_{o}^{2}\Omega_{o}-r_{i}^{2}\Omega_{i})/(r_{o}^{2}-r_{i}^{2})$ and
$C_{2}=(\Omega_{i}-\Omega_{o})r_{i}^{2}r_{o}^{2}/(r_{o}^{2}-r_{i}^{2})$. The
ratio of the inner to outer radii is $r_{i}/r_{o}=0.5$, as in PROMISE. The
fluid has constant viscosity $\nu$ and magnetic diffusivity $\eta$
characterized by Reynolds, $Re=\Omega_{i}r_{i}^{2}/\nu$, and magnetic
Reynolds, $Rm=\Omega_{i}r_{i}^{2}/\eta$, numbers. The magnetic Prandtl number
$Pm=Rm/Re=1.4\times 10^{-6}$ is that of GaInSn used in PROMISE (Stefani et
al., 2009). In this experiment, an imposed azimuthal magnetic field is
current-free, ${\bf B}_{0}=B_{0}(r_{i}/r){\bf e}_{\phi}$, where $B_{0}$ is
constant, which we also adopt here. The effect of the field is quantified by
the Hartmann number $Ha=B_{0}r_{i}/\sqrt{\rho\mu\nu\eta}$, where $\rho$ is the
constant density and $\mu$ is the magnetic permeability. We focus on the
Rayleigh stable regime, $\Omega_{o}/\Omega_{i}>(r_{i}/r_{o})^{2}=0.25$, with
the fixed ratio of the cylinders’ rotation rates $\Omega_{o}/\Omega_{i}=0.26$,
so that only magnetic instabilities can develop in the flow.

About the above equilibrium TC flow, we consider small perturbations of
velocity, ${\bf u}$, total (thermal plus magnetic) pressure, $p$, and magnetic
field, ${\bf b}$, which are all functions of the radius $r$ and depend on time
$t$, azimuthal $\phi$ and axial $z$ coordinates as a normal mode $\propto{\rm
exp}(\gamma t+im\phi+ik_{z}z)$, where $\gamma$ is the (complex) eigenvalue
while $k_{z}$ and integer $m$ are the axial and azimuthal wavenumbers,
respectively. The flow is unstable if the real part (growth rate) of any
eigenvalue is positive, i.e., ${\rm Re}(\gamma)>0$. Henceforth, we normalize
length by $r_{i}$, time by $\Omega_{i}^{-1}$, $\gamma$ and $\Omega(r)$ by
$\Omega_{i}$, ${\bf u}$ by $\Omega_{i}r_{i}$, $p$ by $\rho
r_{i}^{2}\Omega_{i}^{2}$, ${\bf B}_{0}$ by $B_{0}$ and ${\bf b}$ by $Rm\cdot
B_{0}$. Substituting these non-dimensional quantities in the basic non-ideal
MHD equations and linearizing them, we get the perturbation equations also in
the non-dimensional form (e.g., Kirillov et al., 2014; Rüdiger et al., 2018):

$(\gamma+im\Omega){\bf u}+2\Omega{\bf e}_{z}\times{\bf
u}+r\Omega^{\prime}u_{r}{\bf e}_{\phi}=-\nabla
p+\frac{Ha^{2}}{Re}\frac{im}{r^{2}}{\bf
b}-\frac{Ha^{2}}{Re}\frac{2}{r^{2}}b_{\phi}{\bf
e}_{r}+\frac{1}{Re}\nabla^{2}{\bf u},$ (1) $Rm(\gamma+im\Omega){\bf
b}=\frac{im}{r^{2}}{\bf u}+Rm\cdot r\Omega^{\prime}b_{r}{\bf
e}_{\phi}+\frac{2}{r^{2}}u_{r}{\bf e}_{\phi}+\nabla^{2}{\bf b},$ (2)
$\nabla\cdot{\bf u}=0,\hskip 28.45274pt\nabla\cdot{\bf b}=0,$ (3)

where $\Omega^{\prime}=d\Omega/dr$. Equations (1)-(3), together with the
adopted no-slip condition for the velocity and perfect-conductor condition for
the magnetic field at the cylinder surfaces, constitute the eigenvalue
problem, describing AMRI in the magnetized TC flow (Hollerbach et al., 2010).
Solving this problem yields the corresponding dispersion relation
$\gamma(m,k_{z})$ and the radial structure of AMRI modes. Since AMRI is
dominated by non-axisymmetric $m=\pm 1$ modes (Hollerbach et al., 2010;
Rüdiger et al., 2018), we focus on these modes in this paper.

Until now, AMRI has been extensively investigated both with global 1D and
local WKB approaches only as a convective instability (e.g., Hollerbach et
al., 2010; Kirillov et al., 2014; Rüdiger et al., 2018). These studies formed
the basis for understanding the first experimental manifestations of AMRI
(Seilmayer et al., 2014). As mentioned above, however, these experiments
revealed a notable feature of the propagation of AMRI modes – the butterfly
diagram – which is not captured by the conventional treatment of AMRI as a
convective instability. Our main goal is therefore to show, using the concept
of absolute instability, that this butterfly-like propagation is in fact
rooted in the dynamics of AMRI itself rather than being induced by the top and
bottom boundaries (endcaps) of the TC device.

##  3 Convective and absolute instabilities – a short overview

In flow systems, one can distinguish two types of instabilities – convective
and absolute (Huerre & Monkewitz, 1990; Chomaz, 2005). An instability is
convective if the perturbation during growth also propagates in the flow as a
traveling wave packet, so that it decays at large times with respect to a
fixed point in any reference frame, but continues to grow in the frame co-
moving with the group velocity of this packet. By contrast, an instability is
absolute when the perturbation grows unlimited at every point in the flow. In
laboratory experiments, it is usually the absolute instability which is more
relevant, and hence, of greater interest than the convective instability. An
experimental device containing the flow, to which the laboratory frame is
attached, has a finite size. As a result, the convective instability can be
more rapidly carried out from the system before attaining the growth
sufficient for detection, whereas the absolute instability always stays within
the system, exhibiting sustained growth detectable in experiments.

The convective instability is determined by a standard procedure of linear
modal stability analysis, that is, by considering the dispersion relation at
real wavenumbers $k$ and finding real (growth rate) and imaginary (frequency)
parts of $\gamma$ (e.g., Rüdiger et al., 2018). In general, this yields some
non-zero group velocity for the unstable modes. Determining the absolute
instability is, however, more delicate, as it requires finding such a
combination of convectively unstable modes, which forms growing perturbations
(wave packets) with zero group velocity in the laboratory frame, and hence
remains in the flow. This is achieved mathematically by considering complex
wavenumbers $k$111$k$ is the wavenumber in the propagation direction of
convective instability, which is along the $z$-axis in the present case.,
analytically extending the dispersion relation $\gamma(k)$ from real $k$-axis
to complex $k$-plane and finding its saddle points $k_{s}$ (Gailitis &
Freibergs, 1980; Huerre & Monkewitz, 1990; Chomaz, 2005), where the complex
derivative

$\frac{\partial\gamma(k)}{\partial k}\Big{|}_{k=k_{s}}=0,\hskip
28.45274ptk_{s}\in\mathbb{C}.$ (4)

If the real part of the complex $\gamma$ at $k=k_{s}$ is positive, ${\rm
Re}(\gamma(k_{s}))>0$, this indicates the presence of the absolute
instability, where ${\rm Re}(\gamma(k_{s}))$ is its growth rate and
$\omega={\rm Im}(\gamma(k_{s}))$ its frequency. The real ${\rm Re}(k_{s})$ and
imaginary ${\rm Im}(k_{s})$ parts of $k_{s}$ describe, respectively, its
oscillations and exponential increase/decrease in space. The condition (4)
implies that the real group velocity of the absolute instability is zero at
$k_{s}$, i.e., $\partial\omega/\partial{\rm Re}(k)|_{k=k_{s}}=0$, although it
can still have a non-zero phase velocity.

##  4 Convective and absolute types of AMRI

The AMRI is essentially the instability of inertial waves (Kirillov et al.,
2014) and therefore it was natural in previous studies to regard it as a
convective instability and apply a standard technique of modal linear
stability analysis, where the axial wavenumber $k_{z}$ is real. In this paper,
we instead address AMRI within the absolute instability framework, which has
not been considered before, and apply the described above methods, where the
wavenumber $k_{z}$ in the axial $z$-direction, along which AMRI waves travel,
is now complex.

###  4.1 WKB analysis

To obtain an initial insight into the absolute AMRI, we start with the local
short-wavelength WKB analysis, in which the radial dependence of the
perturbations is of the form $\exp(ik_{r}r)$, with a large (real) radial
wavenumber $k_{r}r_{i}\gtrsim 1$. Substituting this in (1)-(3) and using the
inductionless limit $Pm\rightarrow 0$, we get an analytical dispersion
relation (Kirillov et al., 2014):

$\gamma=-im-\frac{k^{2}}{Re}-(m^{2}+2\alpha^{2})\frac{\textit{Ha}^{2}}{k^{2}\textit{Re}}+2\left[\frac{\alpha^{2}Ha^{2}}{k^{2}Re^{2}}\left((\alpha^{2}+m^{2})\frac{Ha^{2}}{k^{2}}+im(Ro+2)\right)-(Ro+1)\alpha^{2}\right]^{1/2}$
(5)

where $\alpha=k_{z}/k,~{}k^{2}=k_{r}^{2}+k_{z}^{2}$,
$Ro=r(2\Omega)^{-1}d\Omega/dr$ is the Rossby number. Using the conversion
formula of Stefani & Kirillov (2015) for $\Omega_{o}/\Omega_{i}=0.26$ in the
TC flow, we obtain $Ro=-0.97$ in the local case, which is below the lower Liu
limit, $Ro_{LLL}=-0.83$, ensuring the presence of AMRI (Kirillov et al.,
2014). Typically, AMRI modes extend over the gap width between the cylinders
(see below), so for the effective radial wavenumber we take
$k_{r}=\pi/(r_{o}-r_{i})$ (Ji et al., 2001), though this choice is still
somewhat arbitrary.

Figure 1: (a)-(e) The areas of growth, ${\rm Re}(\gamma)>0$, in the
$k_{z}$-plane obtained from the WKB dispersion relation (5) for fixed
$Re=1480,m=1$, $Ro=-0.97$ and different $Ha$. The red cross in (c) denotes the
saddle point at $k_{z,s}=(5.08,-1.39)$ – the wavenumber of the absolute AMRI
mode with the growth rate ${\rm Re}(\gamma(k_{z,s}))=0.07$. (f) The marginal
stability curves for the convective (dashed) and absolute (red) AMRI in
$(Ha,Re)$-plane.

Figure 1 shows the solution of the dispersion relation (5) for complex $k_{z}$
at $Re=1480$, $m=1$ and different $Ha$.222The $m=-1$ case is similar except
that the two growth areas swap their ${\rm Im}(k_{z})$ values in the
$k_{z}$-plane, so we do not show it here; see Sec. 4.2. At smaller $Ha$, there
are two separate areas of positive ${\rm Re}(\gamma)\geq 0$, however, neither
crosses yet the ${\rm Re}(k_{z})$-axis (red dashed lines in figure 1),
implying that there is no genuine instability in the flow. With increasing
$Ha$, the lower area first crosses the ${\rm Re}(k_{z})$-axis, indicating the
onset of the convective AMRI, and then extends into the upper half of the
$k_{z}$-plane, which is shown at $Ha=20$ in figure 1(a). At $Ha=27.407$, the
boundaries (${\rm Re}(\gamma)=0$) of these two areas touch each other, forming
a saddle point $k_{z,s}$ according to condition (4), which is seen as a cusp
in figure 1(b). This saddle point, still being marginally stable (${\rm
Re}(\gamma(k_{z,s}))=0$), represents the emerging absolute AMRI. It further
develops by increasing its growth rate and spreading towards the lower half of
the $k_{z}$-plane, as the two areas merge more with increasing $Ha$. For
example, at $Ha=110$ in figure 1(c), the saddle point associated with the
absolute AMRI mode is at $k_{z,s}=(5.08,-1.39)$ (red cross). It has the growth
rate ${\rm Re}(\gamma(k_{z,s}))=0.07$ and frequency ${\rm
Im}(\gamma(k_{z,s}))=-0.5$. At higher $Ha$, the saddle point and hence the
absolute AMRI disappears, but the convective instability along the ${\rm
Re}(k_{z})$-axis can still remain, figure 1(d). Eventually, the convective
AMRI also vanishes after some maximum $Ha$ is exceeded, as the growth areas
move away from the ${\rm Re}(k_{z})$-axis, figure 1(e).

The marginal stability (${\rm Re}(\gamma)=0$) curves for convective (optimized
over ${\rm Re}(k_{z})$) and absolute AMRI are shown in figure 1(f), with the
latter being located inside the former. This implies that the critical $Ha$
and $Re$ for the excitation of the convective AMRI are typically lower than
those for the absolute AMRI. This can be graphically understood from figures
1(a)–1(e). Either of the two growth areas, approaching from upper and lower
halves of the $k_{z}$-plane, first overlaps ${\rm Re}(k_{z})$-axis at certain
$Ha$ for a given $Re$, which lies on the black dashed curve in figure 1(f),
indicating the onset of the convective AMRI. Then, they touch each other at
some larger $Ha$ that lies on the red curve in figure 1(f), indicating the
onset of the absolute AMRI. If the saddle point lies on the ${\rm
Re}(k_{z})$-axis, then the marginal stability curves for convective and
absolute AMRI coincide.

Figure 2: Same as in figure 1, but for the 1D stability analysis at the same
$Re=1480$ and $m=1$. Now the saddle point (red cross) in (c), representing the
absolute AMRI, is at $k_{z,s}=(3.24,-0.2)$ with the growth rate ${\rm
Re}(\gamma(k_{z,s}))=1.9\times 10^{-3}$.

Figure 3: The growth areas with ${\rm Re}(\gamma)>0$, in the $k_{z}$-plane at
$Ha=110$ and $Re=1480$ in the 1D case separately for (a) $m=1$ and (b) $m=-1$
modes. Due to the mirror symmetry of the flow, these areas flip with respect
to ${\rm Im}(k_{z})$-axis when $m$ changes sign and, as a result, the saddle
points (red crosses) for the absolute AMRI are symmetric around this axis:
$k_{z,s}=(3.67,\mp 0.42)$ for $m=\pm 1$. The corresponding growth rate is
${\rm Re}(\gamma(k_{z,s}))=0.00356$ and the frequency is ${\rm
Im}(\gamma(k_{z,s}))=\mp 0.2375$ for $m=\pm 1$. (c) The growth rates of the
convective (black) and absolute (red) AMRI vs. $Ha$ at the same $Re$.

Figure 4: Normalised axial velocity $u_{z}$ eigenfunctions in $(r,z)$-slice at
$Ha=110$, $Re=1480$ pertaining to the absolute AMRI with (a) $m=1$ and (b)
$m=-1$ (represented by red crosses in figure 3) and to the convective AMRI
with (c) $m=1$ and (d) $m=-1$.

###  4.2 1D analysis

Now we present the solution of the 1D linear eigenvalue problem in the
$k_{z}$-plane. The eigenvalues $\gamma$ and the radial structure of the
associated eigenmodes are found using a spectral collocation method based on
Chebyshev polynomials (up to $N=30-40$), whereby equations (1)-(3),
supplemented with the boundary conditions as described in Sec. 2, are reduced
to a large ($4N\times 4N$) matrix eigenvalue problem (see code details in
Hollerbach et al., 2010).

The resulting growth areas in the $k_{z}$-plane at the same $Re=1480$, $m=1$
and different $Ha$ are shown in figures 2(a)–2(e), which exhibit a similar
behaviour as the WKB case above. Initially, for lower $Ha$, there is no AMRI,
since both areas are distant from the ${\rm Re}(k_{z})$-axis (red dashed
line), figure 2(a). When $Ha$ is increased to a critical value of 81.6, the
convective AMRI sets in as the upper area crosses the ${\rm Re}(k_{z})$-axis,
figure 2(b). This crossing point has a cusp-like shape and thus appears to be
the saddle point at the same time, implying that the convective and absolute
AMRI are nearly equivalent at the onset. The growth rate of the absolute AMRI
further increases with $Ha$, as these areas merge. For example, at $Ha=90$ in
figure 2(c), the saddle point, corresponding to the absolute AMRI mode, is at
$k_{z,s}=(3.24,-0.2)$ (red cross) with the growth rate ${\rm
Re}(\gamma(k_{z,s}))=1.9\times 10^{-3}$ and frequency ${\rm
Im}(\gamma(k_{z,s}))=-0.25$. On further increasing $Ha$, first vanishes the
absolute AMRI, figure 2(d), and then the convective AMRI, figure 2(e).

In figure 2(f), we also plot the marginal stability curves for the convective
and absolute AMRI, with the latter being located inside the former, as in the
WKB case, but in the 1D analysis, both are excited at almost similar critical
$Ha$ and $Re$. Furthermore, comparing figures 1 and 2, shows that the
distributions of the growth areas in the $k_{z}$-plane in the WKB and 1D
analyses are qualitatively similar. However, the uncertainty in matching the
constant $Ro$ in the WKB case and radially varying $Ro$ in the 1D case for a
moderate gap width $(r_{0}-r_{i})/r_{i}\sim 1$ (Stefani & Kirillov, 2015) as
well as some ambiguity in choosing the effective $k_{r}$ in the WKB analysis,
makes quantitative differences between the two results unavoidable.

Figures 3(a,b) show areas with ${\rm Re}(\gamma)\geq 0$ in the $k_{z}$-plane
at $Ha=110$ for which the absolute AMRI reaches a maximum growth for a given
$Re=1480$. Since the TC flow with a purely azimuthal field is invariant to
reversing the sign of $z$, the saddle points $k_{z,s}$ and hence the absolute
AMRI at $m=\pm 1$ occur symmetrically with respect to the ${\rm
Re}(k_{z})$-axis: their ${\rm Re}(k_{z,s})=3.67$ and growth rate ${\rm
Re}(\gamma(k_{z,s}))=0.00356$ are the same, while ${\rm Im}(k_{z,s})=\mp 0.42$
and frequency ${\rm Im}(\gamma(k_{z,s}))=\mp 0.2375$ differ only in sign. This
implies that although the absolute AMRI modes have zero group velocity, they
have non-zero phase velocities along and opposite the $z$-axis, respectively,
for $m=-1$ and $m=1$. In figure 3(c), we also plot ${\rm Re}(\gamma)$ for the
convective (optimized over ${\rm Re}(k_{z})$) and absolute AMRI. It is seen
that the convective AMRI has a larger growth rate and occurs for a wider range
of $Ha$ than the absolute AMRI. As noted above, however, the convective and
absolute AMRI are nearly identical at the onset for these parameters.

It is now interesting to look at the spatial structure of the absolute and
convective AMRI modes. Figure 4 shows the normalised eigenfunctions for the
axial velocity $u_{z}$ in $(r,z)$-slice belonging to the absolute and
convective AMRI at $m=\pm 1$, $Ha=110$ and $Re=1480$. It is evident that the
eigenfunctions for these two types of AMRI are characteristically similar
except that due to the complex wavenumber $k_{z,s}$ of the absolute AMRI, its
eigenfunctions increase (for $m=1$) and decrease (for $m=-1$) along the
$z$-axis. This property is important for the interpretation of the simulation
results in Sec. 5. By contrast, the wavenumber of the convective AMRI is real
and hence its eigenfunctions retain a periodic structure with spatially
constant amplitude along the $z$-axis.

##  5 Simulations

We also conducted direct numerical simulation (DNS) of AMRI in the same
magnetized TC flow setup and parameters as described in Sec. 2, except that
the cylinders now have a finite height $L=10r_{i}$ as in PROMISE. We solved
the basic equations of non-ideal MHD (in the same non-dimensional form) using
the open source code library OpenFOAM©333https://www.openfoam.com/,
complemented with a Poisson solver for the determination of the induced
electric potential. The code used here has already been validated in our
previous works (Weber et al., 2013; Seilmayer et al., 2014). The initial
conditions are the previous TC profile $\Omega(r)=C_{1}+C_{2}/r^{2}$ with the
azimuthal magnetic field ${\bf B}_{0}=B_{0}(r_{i}/r){\bf e}_{\phi}$. For the
velocity we set ${\bf u}=r\Omega(r){\bf e}_{\phi}$ (i.e., zero perturbation
velocity) at the top and bottom of the cylinders and no-slip condition at the
cylinder walls, while conducting boundary conditions are applied for the
magnetic field, as in the 1D analysis above. Since the TC flow has a finite
height and non-periodic boundary conditions along the $z$-axis, the AMRI
developing in this case is in fact a global instability. For the same values
of $Ha=110$ and $Re=1480$ as used in Sec. 4.2, we analysed the odd $m$
azimuthal modes of this global AMRI and compare it with the 1D linear analysis
results. We note that although DNS solves full nonlinear equations, it is
still in the linear regime, since the emerging global AMRI mode increases
exponentially in time and has not yet reached nonlinear saturation (see
below).

(a)(b)

Figure 5: Spatio-temporal variation (butterfly diagram) of the axial velocity
$u_{z}$ in $(t,z)$-slice at $Ha=110$ and $Re=1480$ for (a) odd azimuthal modes
of the global AMRI in the DNS and (b) constructed from $m=\pm 1$
eigenfunctions of the absolute AMRI in the 1D linear analysis, which are shown
in figures 4(a,b).

In figure 5(a), we show the temporal evolution of the $z$-profile of the axial
velocity $u_{z,o}(z,t)$ at a radius $r_{0}=1.75$ for odd azimuthal modes of
the global AMRI obtained from the DNS. To remove all even $m$ modes from the
flow, we take the difference of the full axial velocity $u_{z}(r,\phi,z,t)$ at
$\phi=0$ and $\phi=\pi$ for $r_{0}=1.75$, i.e.,
$u_{z,o}(z,t)=u_{z}(r_{0},0,z,t)-u_{z}(r_{0},\pi,z,t)$. Since $m=\pm 1$ are
the most dominant non-axisymmetric modes in the flow, it is appropriate to
compare this velocity profile with that of the absolute AMRI obtained by the
1D linear analysis. In figure 5(b), we show the normalised axial velocity
profile $u_{z}(z,t)$ for the absolute AMRI constructed from the linear
superposition of its eigenfunctions at $m=\pm 1$ (figure 4a,b), taken with
equal weights. Figure 5 shows a remarkable similarity in the butterfly
patterns of the axial velocity resulting from the DNS and 1D analysis.
Evidently, it also demonstrates that the global AMRI in the simulations are
dominated by the $m=\pm 1$ modes444We also checked this by directly
calculating the amplitudes of different $m$ modes in the azimuthal Fourier
transform of the global solution. with the flow structure similar to that in
the 1D analysis. Next, we obtain the frequency, growth rate and axial
wavenumbers of this global AMRI structure.

(a)(b)(c)

Figure 6: (a) Time evolution of $u_{z,o}(z_{0},t)$ for the global AMRI at
$z_{0}=-3$ with (b) the corresponding Hilbert transform. (c) Fits to the axial
structure of this dominant global mode (lilac) obtained using POD at
$r_{0}=1.75$ with two functions $f(z)$ (green) and $g(z)$ (blue), which
correspond to the absolute AMRI modes with $m=-1$ and $m=1$, respectively (see
text).

In figure 6(a), we plot $u_{z,o}$ at $z_{0}=-3$ as a function of time and make
a fit with the expressions of the harmonic form $A\exp(\Gamma
t)\cos(\omega(t-t_{0}))$, obtaining in this way for the growth rate $\Gamma$
and frequency $\omega$ the values $\Gamma=0.00337$ and $\omega=0.252$. We also
conducted a Hilbert transform (figure 6b) to obtain envelopes of the velocity
evolution and instantaneous frequency, which yields growth rate $0.0034$ and
frequency $0.251$ of the global AMRI. Both of these results are in a good
agreement with the 1D linear stability analysis for the absolute AMRI that
yielded the growth rate $0.00356$ and frequency $0.2375$, figures 3(a,b).

We also performed Proper Orthogonal Decomposition (POD) of $u_{z,o}(z,t)$ in
order to get the real and imaginary components of the axial wavenumbers
$k_{z}$ of the global AMRI mode structure in figure 5(a) using the model
reduction library modred (Belson et al., 2014). The obtained amplitude of this
dominant global mode as a function of $z$ is plotted in figure 6(c). As the
global AMRI flow actually consists of the two $m=\pm 1$ modes, each prominent
on either side of the symmetry $z=0$-plane, it is reasonable to make two
separate fits using the expressions $f(z)=A_{1}\exp(-{\rm
Im}(k_{z,1})(z+5))\cos({\rm Re}(k_{z,1})(z-z_{1}))$ and $g(z)=A_{2}\exp(-{\rm
Im}(k_{z,2})(z-5))\cos({\rm Re}(k_{z,2})(z-z_{2}))$, which are also shown in
figure 6(c), each corresponding to either of the two absolute AMRI modes with
$m=\pm 1$. Upon such a fitting, we obtain the values $k_{z,1}=(3.26,0.29)$ and
$k_{z,2}=(3.49,-0.392)$, which agrees reasonably well with our 1D analysis
result $k_{z,s}=(3.67,\mp 0.42)$ from figures 3(a,b). However, it should be
noted that the DNS are done for a finite length TC flow, while in the linear
stability analysis, the cylinder height is assumed to be infinite. Hence, the
minor quantitative difference between the DNS and the 1D linear analysis
results may be attributable to the assumption of infinite length of the
cylinders in the latter approach.

##  6 Conclusion

By performing the WKB and 1D linear stability analyses, we have shown that the
absolute variant of AMRI corresponds to a saddle point in the complex axial
wavenumber plane which, in physical space, appears as the characteristic
solution exponentially increasing in the (negative or positive) axial
direction. This is quite different from, and has a more experimental
significance than, the usually considered convective AMRI with a real axial
wavenumber. The corresponding global instability for the finite height TC flow
(with non-periodic axial boundary conditions) was obtained by DNS, which has
been demonstrated to represent the superposition of two absolute AMRI modes
with $m=\pm 1$ from the 1D stability analysis of infinitely long cylinders. In
other words, the butterfly diagram exhibited by the global AMRI mode can be
understood as arising due to the presence of two dominant absolute AMRI wave
modes with zero group velocities but with respective phase velocities of
propagation directed upwards and downwards. These results further refine the
interpretation of the butterfly-type structure observed in the PROMISE
experiment (Seilmayer et al., 2014, 2020), and will be essential for the
design of new liquid sodium MRI experiment in the new DRESDYN project (Stefani
et al., 2019). Finally, we note that understanding the dynamics of AMRI is
also important for establishing its relevance for the “dead zones” of
accretion (protoplanetary) discs (Gammie, 1996; Lesur et al., 2014), which are
a cold, dense, very weakly ionized and therefore magnetically less active
inner regions of discs, where high resistivity, resulting in $Pm\ll 1$, leads
to the suppression of standard MRI.

[Acknowledgments]This project has received funding from the European Union’s
Horizon 2020 research and innovation programme under the ERC Advanced Grant
Agreement No. 787544 and from the Shota Rustaveli National Science Foundation
of Georgia (SRNSFG, Grant No. FR17-107). We thank R. Hollerbach for kindly
providing the linear 1D code used here and J. Ogbonna for carefully reading
the manuscript. We also thank the anonymous Reviewers for constructive reports
that improved the presentation of our work.

[Declaration of Interests]The authors report no conflict of interest.

## References

  * Balbus & Hawley (1998) Balbus, S. & Hawley, J. 1998 Instability, turbulence, and enhanced transport in accretion disks. Rev. Mod. Phys. 70, 1–53.
  * Belson et al. (2014) Belson, B. A., Tu, J. H. & Rowley, C. W. 2014 Algorithm 945: Modred—a parallelized model reduction library. ACM Trans. Math. Softw. 40, 30.
  * Chomaz (2005) Chomaz, J.-M. 2005 Global instabilities in spatially developing flows: Non-normality and nonlinearity. Ann. Rev. Fluid Mech. 37, 357–392.
  * Gailitis (1996) Gailitis, A. 1996 Design of a liquid sodium MHD dynamo experiment. Magnetohydrodynamics 32, 58–62.
  * Gailitis & Freibergs (1980) Gailitis, A. & Freibergs, Ya. 1980 Nature of the instability of a turbulent dynamo. Magnetohydrodynamics 16, 116–121.
  * Gailitis et al. (2018) Gailitis, A., Gerbeth, G., Gundrum, Th., Lielausis, O., Lipsbergs, G., Platacis, E. & Stefani, F. 2018 Self-excitation in a helical liquid metal flow: the Riga dynamo experiments. J. Plasma. Phys. 84, 735840301.
  * Gammie (1996) Gammie, C. F. 1996 Layered Accretion in T Tauri Disks. Astrophys. J. 457, 355.
  * Hollerbach et al. (2010) Hollerbach, R., Teeluck, V. & Rüdiger, G. 2010 Nonaxisymmetric magnetorotational instabilities in cylindrical Taylor-Couette flow. Phys. Rev. Lett. 104, 044502.
  * Huerre & Monkewitz (1990) Huerre, P. & Monkewitz, P. A. 1990 Local and global instabilities in spatially developing flows. Ann. Rev. Fluid Mech. 22, 473–537.
  * Ji et al. (2001) Ji, H., Goodman, J. & Kageyama, A. 2001 Magnetorotational instability in a rotating liquid metal annulus. Mon. Not. R. Astron. Soc. 325, L1–L5.
  * Kirillov et al. (2014) Kirillov, O. N., Stefani, F. & Fukumoto, Y. 2014 Local instabilities in magnetized rotational flows: a short-wavelength approach. J. Fluid Mech. 760, 591–633.
  * Lesur et al. (2014) Lesur, G., Kunz, M. & Fromang, S. 2014 Thanatology in protoplanetary discs. The combined influence of Ohmic, Hall, and ambipolar diffusion on dead zones. Astron. Astrophys. 566, A56.
  * Liu et al. (2006) Liu, W., Goodman, J., Herron, I. & Ji, H. 2006 Helical magnetorotational instability in magnetized Taylor-Couette flow. Phys. Rev. E 74, 056302.
  * Mamatsashvili et al. (2020) Mamatsashvili, G., Chagelishvili, G., Pessah, M. E., Stefani, F. & Bodo, G. 2020 Zero Net Flux MRI Turbulence in Disks: Sustenance Scheme and Magnetic Prandtl Number Dependence. Astrophys. J. 904, 47.
  * Monchaux et al. (2009) Monchaux, R., Berhanu, M., Aumaitre, S., Chiffaudel, A., Daviaud, F., Dubrulle, B., Ravelet, F., Fauve, S., Mordant, N., Petrelis, F., Bourgoin, M., Odier, P., Pinton, J.-F., Plihon, N. & Volk, R. 2009 The von Karman Sodium experiment: Turbulent dynamical dynamos. Phys. Fluids 21, 035108.
  * Müller et al. (2008) Müller, U., Stieglitz, R., Busse, F. H. & Tilgner, A. 2008 The Karlsruhe two-scale dynamo experiment. C. R. Phys 9, 729–740.
  * Nornberg et al. (2010) Nornberg, M. D., Ji, H., Schartman, E., Roach, A. & Goodman, J. 2010 Observation of magnetocoriolis waves in a liquid metal Taylor-Couette experiment. Phys. Rev. Lett. 104, 074501.
  * Ponomarenko (1973) Ponomarenko, Yu. B. 1973 On the theory of hydromagnetic dynamo. J. Appl. Mech. Tech. Phys. 14, 775–779.
  * Priede & Gerbeth (2009) Priede, J. & Gerbeth, G. 2009 Absolute versus convective helical magnetorotational instability in a Taylor-Couette flow. Phys. Rev. E 79, 046310.
  * Rincon (2019) Rincon, F. 2019 Dynamo theories. J. Plasma Phys. 85, 205850401.
  * Rüdiger et al. (2018) Rüdiger, G., Gellert, M., Hollerbach, R., Schultz, M. & Stefani, F. 2018 Stability and instability of hydromagnetic Taylor-Couette flows. Phys. Rep. 741, 1–89.
  * Seilmayer et al. (2014) Seilmayer, M., Galindo, V., Gerbeth, G., Gundrum, Th., Stefani, F., Gellert, M., Rüdiger, G., Schultz, M. & Hollerbach, R. 2014 Experimental evidence for nonaxisymmetric magnetorotational instability in a rotating liquid metal exposed to an azimuthal magnetic field. Phys. Rev. Lett. 113, 024505.
  * Seilmayer et al. (2020) Seilmayer, M., Ogbonna, J. & Stefani, F. 2020 Convection-caused symmetry breaking of azimuthal magnetorotational instability in a liquid metal Taylor-Couette flow. Magnetohydrodynamics 56, 225.
  * Sisan et al. (2004) Sisan, D., Mujica, N., Tillotson, W., Huang, Y., Dorland, W., Hassam, A., Antonsen, T. & Lathrop, D. P. 2004 Experimental observation and characterization of the magnetorotational instability. Phys. Rev. Lett. 93, 114502\.
  * Stefani et al. (2008) Stefani, F., Gailitis, A. & Gerbeth, G. 2008 Magnetohydrodynamic experiments on cosmic magnetic fields. Z. Angew. Math. Mech. 88, 930–954.
  * Stefani et al. (2019) Stefani, F., Gailitis, A., Gerbeth, G., Giesecke, A., Gundrum, Th., Rüdiger, G., Seilmayer, M. & Vogt, T. 2019 The DRESDYN project: liquid metal experiments on dynamo action and magnetorotational instability. Geophys. Astrophys. Fluid Dyn. 113, 51–70.
  * Stefani et al. (1999) Stefani, F., Gerbeth, G. & Gailitis, A. 1999 Velocity profile optimization for the riga dynamo experiment. In Transfer phenomena in magnetohydrodynamics and electroconducting flows (ed. A. Alemany, P. Marty & J. P. Thibault), Fluid mechanics and its applications, vol. 51, pp. 31–44.
  * Stefani et al. (2009) Stefani, F., Gerbeth, G., Gundrum, Th., Hollerbach, R., Priede, J., Rüediger, G. & Szklarski, J. 2009 Helical magnetorotational instability in a Taylor-Couette flow with strongly reduced Ekman pumping. Phys. Rev. E 80, 066303\.
  * Stefani et al. (2006) Stefani, F., Gundrum, Th., Gerbeth, G., Rüdiger, G., Schultz, M., Szklarski, J. & Hollerbach, R. 2006 Experimental evidence for magnetorotational instability in a Taylor-Couette flow under the influence of a helical magnetic field. Phys. Rev. Lett. 97, 184502.
  * Stefani & Kirillov (2015) Stefani, F. & Kirillov, O. N. 2015 Destabilization of rotating flows with positive shear by azimuthal magnetic fields. Phys. Rev. E 92, 051001.
  * Weber et al. (2013) Weber, N., Galindo, V., Stefani, F., Weier, T. & Wondrak, Th. 2013 Numerical simulation of the Tayler instability in liquid metals. New J. Phys. 15, 043034.