# Preparing ground states of the XXZ model using the quantum annealing with
inductively coupled superconducting flux qubits

Takashi Imoto  Research Center for Emerging Computing Technologies, National
Institute of Advanced Industrial Science and Technology (AIST), 1-1-1 Umezono,
Tsukuba, Ibaraki 305-8568, Japan.  Yuya Seki  Graduate School of Science and
Technology, Keio University, Hiyoshi 3-14-1, Kohoku-ku, Yokohama 223-8522,
Japan.  Yuichiro Matsuzaki  Research Center for Emerging Computing
Technologies, National Institute of Advanced Industrial Science and Technology
(AIST), 1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, Japan.
<EMAIL_ADDRESS>
###### Abstract

Preparing ground states of Hamiltonians is important in the condensed matter
physics and the quantum chemistry. The interaction Hamiltonians typically
contain not only diagonal but also off-diagonal elements. Although quantum
annealing provides a way to prepare a ground state of a Hamiltonian, we can
only use the Hamiltonian with Ising interaction by using currently available
commercial quantum annealing devices. In this work, we propose a quantum
annealing for the XXZ model, which contains both Ising interaction and energy-
exchange interaction, by using inductively coupled superconducting flux
qubits. The key idea is to use a recently proposed spin-lock quantum annealing
where the qubits are driven by microwave fields. As long as the rotating wave
approximation is valid, the inductive coupling between the superconducting
flux qubits produces the desired Hamiltonian in the rotating frame, and we can
use such an interaction for the quantum annealing while the microwave fields
driving play a role of the transverse fields. To quantify the performance of
our scheme, we implement numerical simulations, and show that we can prepare
ground states of the two-dimensional Heisenberg model with a high fidelity.

##  1 Introduction

Quantum annealing (QA) is a method to solve the combinatorial optimization
problem [1, 2, 3]. Combinatorial optimization problems can be mapped into
problems of finding the ground state of an Ising model. In this case, the
Ising model is the problem Hamiltonian. On the other hand, we use the
Hamiltonian of transverse fields to prepare an initial state, and this is
called a drive Hamiltonian. In the QA, we prepare a ground state of the drive
Hamiltonian, and then let the system evolve by a time dependent Hamiltonian to
change from the drive Hamiltonian to the problem Hamiltonian in an adiabatic
way. This provides a way to prepare a ground-state of the problem Hamiltonian
as long as an adiabatic condition is satisfied.

There are many other applications of the quantum annealing in addition to
combinatorial optimization problems. Firstly, quantum annealing can be used
for machine learning. A way to learn binary machine learning models using the
quantum annealing is reported [4, 5]. There are applications of the quantum
annealing for clustering [6, 7]. A method of using the quantum annealing to
perform calculations for the topological data analysis(TDA) was reported [8].
Quantum annealing has been applied as a sampling machine for restricted
Boltzmann machines with deep layers [9, 10, 11, 12, 13]. Secondly, quantum
annealing is known to be useful for quantum chemistry calculation. The
Hamiltonian of molecules are written in the second quantized form. There are
known methods to convert the second quantized form of the Hamiltonian into the
spin Hamiltonian [14, 15, 16, 17, 18]. In quantum chemical calculations, it is
important to obtain the energy of the molecules with high accuracy. When the
energy is known to the so called chemical accuracy, a chemical reaction can be
predicted [19]. In fact, calculations to determine the energy of molecules
using quantum annealing have been reported [20, 21, 22, 23]. In addition, the
method to measure similarity of the structures of molecules using the quantum
annealing is proposed [24]. Finally, the quantum annealing can be applied to
simulate condensed matter phenomena The $\mathbb{Z}_{2}$ spin liquid is
emulated with the quantum annealing [25]. It has been reported that the
KT(Kosterlitz–Thouless) phase transition was observed using a quantum
annealing on the d-wave system [26]. Simulations of the Shastry-Sutherland
Ising, a frustrated system, were performed using quantum annealing [27]. A
phase transition of three dimensional transverse Ising model is detected using
the quantum annealing machine [28]. It is worth mentioning that, in the
previous demonstrations, the problem Hamiltonian was mapped into an Ising
model, and this typically requires many ancillary qubits.

Superconducting flux qubits (FQs) are considered as promising systems for
quantum information processing. The FQ is composed of a superconducting loop
with Josephson junctions [29, 30, 31, 32]. It is possible to tune the energy
gap of the FQ and the interaction between FQs [33, 34]. There are many
attempts to improve the coherence time of the FQs [35, 36, 37]. Moreover,
recently, a new design of the FQs called a capacitivly shunted flux qubit with
a long coherence time was proposed and demonstrated [38, 39, 40, 41]. The FQ
can be used as a sensor to detect magnetic fields, electron spins, or defects
[42, 43, 44, 45]. The FQs have advantages in a scalability, and actually it is
possible to fabricate thousands of the FQs even in the current technology
[46].

The company D-Wave Systems, Inc. realizes a quantum annealer with thousands of
the FQs [47, 48]. The energy bias, energy gap, and the coupling between the
FQs can be tuned by changing the applied magnetic flux. With currently
available commercial quantum annealing devices, we can use only transverse
field Ising model as the Hamiltonian where the flux qubits are inductively
coupled. Recently, two flux qubits are coupled by using two degree of
freedoms, charge and flux [49]. In this case, not only Ising coupling but also
energy-exchange interaction such as $\hat{\sigma}_{x}\hat{\sigma}_{x}$ (or
$\hat{\sigma}_{y}\hat{\sigma}_{y}$) can be realized. However, this approach
requires a complicated setup. It is preferable if we could construct the
energy-exchange terms with less degree of freedom of the FQs.

In this paper, we propose a method to prepare a ground-state of the
Hamiltonian that contains both Ising interaction and energy exchange
interaction, which is the XXZ model, by using QA with the inductively coupled
FQs. The key idea is to adopt a spin lock quantum annealing [50, 51, 52]. When
we drive the FQs with the microwave fields, we can rewrite the Hamiltonian in
a rotating frame with a frequency of the microwave fields. By using a rotating
wave approximation, the inductive coupling between the FQs provides not only
Ising interaction but also the energy exchange interaction. The driving with
the microwave fields play a role of the driver Hamiltonian of the QA, and we
can obtain a ground state of a Hamiltonian that contains both Ising and
energy-exchange interaction by adiabatically turning off the driving fields.
Such a preparation of the ground state allows us to calculate an arbitrary
correlation function for the ground state of the XXZ model, which is useful in
the condensed matter physics [53]. To quantify the performance of our scheme,
we implement numerical simulations to prepare ground states of the two-
dimensional Heisenberg model, and shows that we can achieve a high fidelity.

The paper is structured as follow. In Section 2, we review the conventional
quantum annealing and spin lock quantum annealing. In Section 3, we show the
method to create the energy exchange interaction by using the inductivly
coupled FQs for the QA. In Section 4, we perform numerical simulations to
evaluate the performance of our method. In section 5, we summarize and discuss
our results.

##  2 The previous work

In this section, let us review the conventional quantum annealing with DC
fields and recently proposed quantum annealing with spin-lock technique.

###  2.1 Quantum annealing

In this subsection, we review the conventional quantum annealing. The
annealing Hamiltonian is written as

$H(t)=\frac{t}{T}H_{P}+\biggl{(}1-\frac{t}{T}\biggr{)}H_{D}$ (1)

where $T$ is annealing time, $H_{P}$ is the problem Hamiltonian, and $H_{D}$
is the drive Hamiltonian. We define the drive Hamiltonian as follows

$H_{D}=-B\sum_{j=1}^{L}\hat{\sigma}_{j}^{(x)}$ (2)

where $B$ denotes an amplitude of the DC transverse fields. After preparing a
ground state of the drive Hamiltonian, we let the state evolve by the
annealing Hamiltonian. If this dynamics is adiabatic, the ground state of the
drive Hamiltonian is changed into that of the problem Hamiltonian.

###  2.2 Spin lock quantum annealing

We review the spin-lock quantum annealing [50, 51, 52]. Unlike the
conventional QA with the DC transverse fields, we use AC driving fields for
the drive Hamiltonian when we adopt the spin-lock QA. Since the previous
studies focus only on Ising interaction for the problem Hamiltonian, we
consider such a case here. In the case of the spin lock quantum annealing, the
annealing Hamiltonian is given by

$H(\lambda)=B\sum_{j=1}^{L}\cos(\omega
t)\hat{\sigma}_{j}^{(x)}+\sum_{j=1}^{L}\frac{\omega}{2}\hat{\sigma}_{j}^{(z)}+(1-\lambda)H_{P}$
(3)

where $B$ denotes an amplitude of AC field, $\hat{\sigma}_{j}^{(k)}(k=x,y,z)$
denotes the Pauli operator acting on the $j$-th site of the spin chain, $L$ is
the site number, and $H_{P}$ is the problem Hamiltonian. By moving to a
rotating frame defined by $U$, the Hamiltonian is rewritten as

$H^{\prime}=UHU^{{\dagger}}-iU^{{\dagger}}\frac{dU}{dt}$ (4)

Since we drive the qubits with the AC driving fields, we choose the operator
$U$ as follows

$U=\exp\biggl{(}-i\frac{\omega}{2}t\sum_{j=1}^{L}\hat{\sigma}_{j}^{(z)}\biggr{)}.$
(5)

Here, we have

$U^{{\dagger}}H_{P}U=H_{P}$ (6)

and we can calculate as follows.

$\displaystyle U\cos({\omega t})\hat{\sigma}_{j}^{(x)}U^{{\dagger}}$
$\displaystyle=$ $\displaystyle\frac{1}{2}\Bigl{(}(e^{2i\omega
t}+1)\hat{\sigma}_{j}^{(+)}+(e^{-2i\omega t}+1)\hat{\sigma}_{j}^{(-)}\Bigr{)}$
(7) $\displaystyle\approx$ $\displaystyle\frac{1}{2}\hat{\sigma}_{j}^{(x)}.$

where we use the rotating wave approximation. Throughout of this paper, we use
the notation $\approx$ to represent the rotating wave approximation. Using the
relation

$-iU^{{\dagger}}\frac{dU}{dt}=-\frac{\omega}{2}\sum_{j=1}^{L}\hat{\sigma}_{j}^{(z)},$
(8)

we obtain the annealing Hamiltonian in the rotating frame as follows.

$H^{\prime}\approx\sum_{j=1}^{L}\frac{\lambda}{2}\hat{\sigma}_{j}^{(x)}+(1-\lambda)H_{P}$
(9)

Therefore, we obtain the annealing Hamiltonian where the drive Hamiltonian is
the transverse field in the rotating frame. It is worth mentioning that a
violation of the the rotating wave approximation due to strong driving could
cause an error in the preparation of the ground state [52].

##  3 Realization of the XXZ model using flux qubit

Here, we describe our scheme to prepare a ground state of the Hamiltonian with
Ising interaction and energy-exchange interaction by using the spin-lock QA.
Here, we consider a case that the flux qubits are inductively coupled. The
Hamiltonian of the flux qubit [54, 55, 31] or capacitively shunted FQs [39] is
given by

$H=\frac{\epsilon}{2}\sum_{j=1}^{L}\hat{\sigma}_{z}^{(j)}+\frac{\Delta_{\rm{G}}}{2}\sum_{j=1}^{L}\hat{\sigma}_{x}^{(j)}+\lambda\cos{\omega
t}\sum_{j=1}^{L}\hat{\sigma}_{y}^{(j)}+g\sum_{j=1}^{L-1}\hat{\sigma}_{z}^{(j)}\hat{\sigma}_{z}^{(j+1)}.$
(10)

where $\epsilon$ denotes an energy bias, $\Delta_{\rm{G}}$ denotes an energy
gap, $g$ denotes a strength of the inductive interactions, and $\lambda$
denotes a Rabi frequency. First, we use the single qubit operator $U_{y}\equiv
e^{-i\theta\sum_{j=1}^{L}\hat{\sigma}_{y}^{(j)}}$ to rotate this Hamiltonian
by $\theta$ in the y-direction.

$\displaystyle H^{\prime}$ $\displaystyle=$ $\displaystyle U_{y}^{-1}HU_{y}$
$\displaystyle=$
$\displaystyle\biggl{(}\frac{\Delta_{\rm{G}}}{2}\cos{\theta}-\frac{\epsilon}{2}\sin{\theta}\biggr{)}\sum_{j=1}^{L}\hat{\sigma}_{x}^{(j)}+\biggl{(}\frac{\epsilon}{2}\cos{\theta}+\frac{\Delta_{\rm{G}}}{2}\sin{\theta}\biggr{)}\sum_{j=1}^{L}\hat{\sigma}_{z}^{(j)}+\lambda\cos{\omega
t}\sum_{j=1}^{L}\hat{\sigma}_{y}^{(j)}$ $\displaystyle+$ $\displaystyle
g\sum_{j=1}^{L-1}\biggl{(}\cos^{2}{\theta}\hat{\sigma}_{z}^{(j)}\hat{\sigma}_{z}^{(j+1)}+\sin^{2}{\theta}\hat{\sigma}_{x}^{(j)}\hat{\sigma}_{x}^{(j+1)}-\cos{\theta}\sin{\theta}(\hat{\sigma}_{z}^{(j)}\hat{\sigma}_{x}^{(j+1)}+\hat{\sigma}_{x}^{(j)}\hat{\sigma}_{z}^{(j+1)})\biggr{)}$

In order to remove the term of $\sum_{j=1}^{L}\hat{\sigma}_{x}^{(j)}$, we
determine the $\theta$ as follows

$\displaystyle\sin{\theta}=\frac{\Delta_{\rm{G}}}{\sqrt{\epsilon^{2}+\Delta_{\rm{G}}^{2}}},$
(12)
$\displaystyle\cos{\theta}=\frac{\epsilon}{\sqrt{\epsilon^{2}+\Delta_{\rm{G}}^{2}}}.$
(13)

Next, the Hamiltonian in the interaction picture is given by

$\displaystyle H^{\prime\prime}=U^{-1}H^{\prime}U-iU^{-1}\frac{dU}{dt}$ (14)

where the single qubit operator $U$ is defined as

$\displaystyle U=e^{-\frac{i\omega
t}{2}\sum_{j=1}^{L}\hat{\sigma}_{z}^{(j)}}.$ (15)

Using the relation

$\displaystyle-
iU^{-1}\frac{dU}{dt}=-\frac{\omega}{2}\sum_{j=1}^{L}\hat{\sigma}_{z}^{(j)},$

we obtain

$\displaystyle H^{\prime\prime}$ $\displaystyle=$
$\displaystyle\lambda\cos{\omega
t}U^{-1}\sum_{j=1}^{L}\hat{\sigma}_{y}^{(j)}U+\biggl{(}\frac{\epsilon}{2}\cos{\theta}+\frac{\Delta_{\rm{G}}}{2}\sin{\theta}-\frac{\omega}{2}\biggr{)}U^{-1}\sum_{j=1}^{L}\sigma_{z}^{(j)}U$
(16) $\displaystyle+$
$\displaystyle\sum_{j=1}^{L-1}g\biggl{(}\cos^{2}{\theta}U^{-1}\hat{\sigma}_{z}^{(j)}\hat{\sigma}_{z}^{(j+1)}U+\sin^{2}{\theta}U^{-1}\hat{\sigma}_{x}^{(j)}\hat{\sigma}_{x}^{(j+1)}U$
$\displaystyle-$
$\displaystyle\cos{\theta}\sin{\theta}U^{-1}(\hat{\sigma}_{z}^{(j)}\hat{\sigma}_{x}^{(j+1)}+\hat{\sigma}_{x}^{(j)}\hat{\sigma}_{z}^{(j+1)})U\biggr{)}.$

We calculate each term of this Hamiltonian. The first term of this Hamiltonian
is

$\displaystyle\lambda\cos{\omega t}U^{-1}\hat{\sigma}_{y}^{(j)}U$
$\displaystyle=$ $\displaystyle\frac{\lambda}{4i}\biggl{(}e^{2i\omega
t}\hat{\sigma}_{+}^{(j)}+\hat{\sigma}_{+}^{(j)}-e^{-2i\omega
t}\hat{\sigma}_{-}^{(j)}-\hat{\sigma}_{-}^{(j)}\biggr{)}$ (17)
$\displaystyle\approx$ $\displaystyle\frac{\lambda}{2}\hat{\sigma}_{y}^{(j)}$

by rotating wave approximation. We obtain the following

$\displaystyle U^{-1}\hat{\sigma}_{x}^{(j)}\hat{\sigma}_{x}^{(j+1)}U$
$\displaystyle=$ $\displaystyle e^{2i\omega
t}\hat{\sigma}_{+}^{(j)}\hat{\sigma}_{+}^{(j+1)}+e^{-2i\omega
t}\hat{\sigma}_{-}^{(j)}\hat{\sigma}_{-}^{(j+1)}+\hat{\sigma}_{+}^{(j)}\hat{\sigma}_{-}^{(j+1)}+\hat{\sigma}_{-}^{(j)}\hat{\sigma}_{+}^{(j+1)}$
(18) $\displaystyle\approx$
$\displaystyle\hat{\sigma}_{+}^{(j)}\hat{\sigma}_{-}^{(j+1)}+\hat{\sigma}_{-}^{(j)}\hat{\sigma}_{+}^{(j+1)}.$

Also, we obtain

$\displaystyle
U^{-1}(\hat{\sigma}_{z}^{(j)}\hat{\sigma}_{x}^{(j+1)}+\hat{\sigma}_{x}^{(j)}\hat{\sigma}_{z}^{(j+1)})U$
$\displaystyle=$ $\displaystyle\hat{\sigma}_{z}^{(j)}(e^{i\omega
t}\hat{\sigma}_{+}^{(j+1)}+e^{-i\omega t}\hat{\sigma}_{-}^{(j+1)})+(e^{i\omega
t}\hat{\sigma}_{+}^{(j)}+e^{-i\omega
t}\hat{\sigma}_{-}^{(j)})\hat{\sigma}_{z}^{(j+1)}$ (19) $\displaystyle\approx$
$\displaystyle 0.$

Finally, the other terms are unaffected by the operator $U$. Therefore,
substituting the Eq. (17), (18) and, (19) into the Eq. (16), we obtain the

$\displaystyle H^{\prime\prime}\approx\frac{\lambda}{2}$
$\displaystyle\sum_{j=1}^{L}\hat{\sigma}_{y}^{(j)}+\Delta\omega\sum_{j=1}^{L}\hat{\sigma}_{z}^{(j)}$
(20)
$\displaystyle+g\sum_{j=1}^{L-1}\biggl{(}\frac{\epsilon^{2}}{\epsilon^{2}+\Delta_{\rm{G}}^{2}}\hat{\sigma}_{z}^{(j)}\hat{\sigma}_{z}^{(j+1)}+\frac{\Delta_{\rm{G}}^{2}}{\epsilon^{2}+\Delta_{\rm{G}}^{2}}(\hat{\sigma}_{+}^{(j)}\hat{\sigma}_{-}^{(j+1)}+\hat{\sigma}_{-}^{(j)}\hat{\sigma}_{+}^{(j+1)})\biggr{)}\
\ \ \ \ \ \ \ \ \ \ \ \ $

where the detuning is defined as
$\delta\omega=\biggl{(}\frac{\epsilon^{2}}{2\sqrt{\epsilon^{2}+\Delta_{\rm{G}}^{2}}}+\frac{\Delta_{\rm{G}}^{2}}{2\sqrt{\epsilon^{2}+\Delta_{\rm{G}}^{2}}}-\frac{\omega}{2}\biggr{)}$.
By rotating the Hamiltonian by $\pi/2$ along the z axis, we obtain the
Hamiltonian with the transverse field parallel to the x axis. This means that,
by changing $\lambda$ and $g$, we can perform the QA where the driver
Hamiltonian is the transverse field and the problem Hamiltonian is the XXZ
model. Throughout of our paper, we take $\lambda(t)=\lambda(0)\frac{t}{T}$ and
$g(t)=g(0)(1-\frac{t}{T})$. Moreover, since we have not only the transverse
field parallel to the x axis but also the longitudinal field, we could
construct a driver Hamiltonian with twisted fields, which is shown to be
useful to obtain a ground state with the QA [56, 57]. However, in our paper,
we focus on the resonant case with $\delta\omega=0$.

##  4 Application and numerical simulation

In this section, in order to evaluate the performance of our scheme, we
perform numerical simulations. For the problem Hamiltonians, we consider an
example of a two dimensional anisotropic Heisenberg model in condensed matter
physics.

The problem Hamiltonian is given by

$H_{\rm{P}}=J\sum_{\braket{j,k}}\biggl{(}\sigma_{j}^{(x)}\sigma_{k}^{(x)}+\sigma_{j}^{(y)}\sigma_{k}^{(y)}+\Delta\sigma_{j}^{(z)}\sigma_{k}^{(z)}\biggr{)}.$
(21)

where $\Delta$ denotes the anisotropic parameter, and $J$ denotes the coupling
constant. This is called the ferro magnetic model when $J<0$, and anti-ferro
magnetic model when $J>0$. To realize this problem Hamiltonian in the spin
lock we set $J=g\frac{\epsilon^{2}}{\epsilon^{2}+\Delta_{\rm{G}}^{2}}$ and
$\Delta=\frac{\Delta_{\rm{G}}^{2}}{\epsilon^{2}}$. We consider a square
lattice of size $2\times 3$, as shown in Figure.1.

Figure 1: Schematic of our model.

About the coupling constants $J$ and anisotropic parameters $\Delta$, we
consider the following four cases.

$(J,\Delta)=(1,1.7),(1,0.7,),(-1,1.7,),(-1,0.7).$ (22)

where the unit of the value of $J$ is GHz, as shown in Figure.2.

Figure 2:  A conceptual diagram of the anisotropic Heisenberg model.

The ground state was found to be doubly degenerate in the case of ferro-
magnetic coupling ($J=-1$) and $\Delta=1.7$. On the other hand, the ground
state is not degenerate in the other cases. Throughout of our manuscript, when
we plot a fidelity $F$, we calculate $F=\sum_{n=1}^{m}\langle
g_{n}|\rho(t)|g_{n}\rangle$ where $m$ denotes the degeneracy,
$\\{|g_{n}\rangle\\}_{n=1}^{m}$ denotes the ground state, and $\rho(t)$
denotes a quantum state after the QA.

(a) $J=-1$ and $\Delta=0.7$

(b) $J=-1$ and $\Delta=1.7$

(c) $J=1$ and $\Delta=0.7$

(d) $J=1$ and $\Delta=1.7$

Figure 3:  Fidelity between the ground state of the two-dimensional
anisotropic Heisenberg model and the state during QA against time. The
parameters for the simulation are chosen as $T=1000$ ns, $\lambda=1$ and,
$\omega=2.0,5.0,10.0,20.0$ GHz. We adopt the uniform microwave driving (DC
transverse fields) as the drive Hamiltonian for our (conventional) scheme.

We plot the fidelity against time for each frequency $\omega$ with numerical
simulations in the Figure 3. For these simulations, the annealing time $T$ is
set to be $1000$ [ns], the driving strength $\lambda$ is set to be 1 GHz, and
the frequencies $\omega$ are set to be $2.0$, $5.0$, $10.0$, $20.0$ GHz. For
the ferro magnetic case, as we increase the frequency $\omega$, the fidelity
increases. This is consistent with the fact that the rotating wave
approximation becomes more accurate as the qubit frequency increases. Thus,
after the preparation of the ground state with a high fidelity, we can measure
an arbitrary correlation function for the ground state of the Hamiltonian,
which is useful in the condensed matter physics.

In the case of anti-ferro magnetic ($J=1$ GHz), for $\Delta=0.7$ and
$\Delta=1.7$, the fidelities are zero up to numerical errors. Our
considerations suggest that there is a symmetry during QA for these
parameters. In general, when there is a symmetry, the total Hamiltonian
commutes with some observables. In this case, the Hamiltonian can be block
diagonalized into sectors by applying a suitable unitary operator. As long as
the initial state belongs to a different sector from that to which the target
ground state belongs, we cannot prepare a ground state even after the QA by
keeping the adiabatic condition because a level crossing occurs. Actually, in
our case, we confirm that an operator
$U_{swap}^{(1,2)}U_{swap}^{(3,4)}U_{swap}^{(5,6)}$ commutes with both the
transverse field and the ferromagnetic XXZ model Hamiltonian where
$U_{swap}^{(j,k)}$ denotes a swap gate between the $j$-th qubit and the $k$-th
qubit. We calculate the sector of
$U_{swap}^{(1,2)}U_{swap}^{(3,4)}U_{swap}^{(5,6)}$ for the ground states of
the transvers field and the ferromagnetic XXZ model, and we obtain
$\bra{E_{0}^{(TF)}}U_{swap}^{(1,2)}U_{swap}^{(3,4)}U_{swap}^{(5,6)}\ket{E_{0}^{(TF)}}=1$
and
$\bra{E_{0}^{(XXZ)}}U_{swap}^{(1,2)}U_{swap}^{(3,4)}U_{swap}^{(5,6)}\ket{E_{0}^{(XXZ)}}=-1$
where $\ket{E_{0}^{(TF)}}$ denotes the ground state of the transvers field and
$\ket{E_{0}^{(XXZ)}}$ denotes the ground state of the XXZ model. Therefore, we
can confirm that $\ket{E_{0}^{(TF)}}$ and $\ket{E_{0}^{(XXZ)}}$ belong to
different sectors, and this is the reason why the QA fails in our numerical
simulations.

$J_{1}$ | $-0.9052919617126958$
---|---
$J_{2}$ | $0.2500243810835503$
$J_{3}$ | $-1.931378367720707$
$J_{4}$ | $-0.007622719480759765$
$J_{5}$ | $-1.259154537693434$
$J_{6}$ | $-1.261510983981174$
Table 1: Coefficients of random transverse field $\\{J_{j}\\}_{j=1}^{6}$ used
in Fig. 4. The unit of these values is GHz, as described in the main text.

In such a case, we apply not homogeneous but randomized microwave driving (DC
transverse fields) to break the symmetry for our (conventional) scheme.
Moreover, we plot the fidelities against time when we apply the randomized
driving Hamiltonian in the Figure. 4. This demonstrates that we can prepare a
ground state with a high fidelity by using the randomized fields during QA
wihtout level crossing. Also, our results show that a proper choice of the
driver Hamiltonian to break a symmetry of the Hamiltonian is important in QA,
which has been overlooked in previous works.

(a) $J=1$ and $\Delta=0.7$

(b) $J=1$ and $\Delta=1.7$

Figure 4: We plot the fidelity between the ground state of the two-dimensional
anisotropic Heisenberg model and the state during QA against time. The
parameters for the simulation are chosen as , $T=1000$ ns, and
$\omega=2.0,5.0,10.0,20.0$ GHz. In order to break the symmetry, we adopt the
random microwave driving (DC transverse field) as the drive Hamiltonian for
our (conventional) scheme.

##  5 Conclusion

In conclusion, we propose the way to prepare a ground state of the XXY model
with the QA by using inductively coupled superconducting flux qubits. The key
idea is to use a recently proposed spin-lock quantum annealing. We drive the
flux qubits by the microwave driving, and we can construct the Hamiltonian of
QA for the XXY model in the rotating frame. By choosing the problem
Hamiltonian as the anisotropic Heisenberg model, we can prepare a ground state
of the anisotropic Heisenberg model with our scheme unless there is a level
crossing. Importantly, when there is a symmetry during QA, the Hamiltonian can
be block diagonalized into sectors. As long as the initial state belongs to a
different sector from that to which the target ground state belongs, we cannot
prepare the ground state with QA because of the level crossing. In this case,
we show that, by adding randomized fields to break the symmetry, we can avoid
the problem of the level crossing. Our numerical simulation shows that our
spin-lock QA actually provides a practical way to prepare a ground state. Our
scheme is useful to know the correlation function of the relevant Hamiltonians
in the condensed matter physics.

Note added.—While preparing our manuscript, we became aware of a related work
that also proposes an adiabatic scheme to prepare a ground state of a
Hamiltonian (that contains non-diagonal terms) by applying additional fields
to break the symmetry [58].

We thank S. Kawabata and A. Yoshinaga for a useful advice. This work was
supported by MEXT’s Leading Initiative for Excellent Young Researchers and JST
PRESTO (Grant No. JPMJPR1919), Japan. This paper is partly based on the
results obtained from a project, JPNP16007, commissioned by the New Energy and
Industrial Technology Development Organization (NEDO), Japan.

## References

## References

  * [1] Tadashi Kadowaki and Hidetoshi Nishimori.  Quantum annealing in the transverse ising model.  Physical Review E, 58(5):5355, 1998.
  * [2] Edward Farhi, Jeffrey Goldstone, Sam Gutmann, and Michael Sipser.  Quantum computation by adiabatic evolution.  arXiv preprint quant-ph/0001106, 2000.
  * [3] Edward Farhi, Jeffrey Goldstone, Sam Gutmann, Joshua Lapan, Andrew Lundgren, and Daniel Preda.  A quantum adiabatic evolution algorithm applied to random instances of an np-complete problem.  Science, 292(5516):472–475, 2001.
  * [4] Prasanna Date and Thomas Potok.  Adiabatic quantum linear regression.  2020\.
  * [5] Michele Sasdelli and Tat-Jun Chin.  Quantum annealing formulation for binary neural networks.  arXiv preprint arXiv:2107.02751, 2021.
  * [6] Kenichi Kurihara, Shu Tanaka, and Seiji Miyashita.  Quantum annealing for clustering.  In 25th Conference on Uncertainty in Artificial Intelligence, UAI 2009, 2009.
  * [7] Vaibhaw Kumar, Gideon Bass, Casey Tomlin, and Joseph Dulny.  Quantum annealing for combinatorial clustering.  Quantum Information Processing, 17(2):1–14, 2018.
  * [8] Jesse J Berwald, Joel M Gottlieb, and Elizabeth Munch.  Computing wasserstein distance for persistence diagrams on a quantum computer.  arXiv preprint arXiv:1809.06433, 2018.
  * [9] Steven H Adachi and Maxwell P Henderson.  Application of quantum annealing to training of deep neural networks.  arXiv preprint arXiv:1510.06356, 2015.
  * [10] Max Wilson, Thomas Vandal, Tad Hogg, and Eleanor G Rieffel.  Quantum-assisted associative adversarial network: Applying quantum annealing in deep learning.  Quantum Machine Intelligence, 3(1):1–14, 2021.
  * [11] Richard Y Li, Tameem Albash, and Daniel A Lidar.  Limitations of error corrected quantum annealing in improving the performance of boltzmann machines.  Quantum Science and Technology, 5(4):045010, 2020.
  * [12] Hartmut Neven, Vasil S Denchev, Geordie Rose, and William G Macready.  Training a binary classifier with the quantum adiabatic algorithm.  arXiv preprint arXiv:0811.0416, 2008.
  * [13] Walter Winci, Lorenzo Buffoni, Hossein Sadeghi, Amir Khoshaman, Evgeny Andriyash, and Mohammad H Amin.  A path towards quantum advantage in training deep generative models with quantum annealers.  Machine Learning: Science and Technology, 1(4):045028, 2020.
  * [14] Sergey B Bravyi and Alexei Yu Kitaev.  Fermionic quantum computation.  Annals of Physics, 298(1):210–226, 2002.
  * [15] Frank Verstraete and J Ignacio Cirac.  Mapping local hamiltonians of fermions to local hamiltonians of spins.  Journal of Statistical Mechanics: Theory and Experiment, 2005(09):P09012, 2005.
  * [16] Jacob T Seeley, Martin J Richard, and Peter J Love.  The bravyi-kitaev transformation for quantum computation of electronic structure.  The Journal of chemical physics, 137(22):224109, 2012.
  * [17] Andrew Tranter, Sarah Sofia, Jake Seeley, Michael Kaicher, Jarrod McClean, Ryan Babbush, Peter V Coveney, Florian Mintert, Frank Wilhelm, and Peter J Love.  The b ravyi–k itaev transformation: Properties and applications.  International Journal of Quantum Chemistry, 115(19):1431–1441, 2015\.
  * [18] Rongxin Xia, Teng Bian, and Sabre Kais.  Electronic structure calculations and the ising hamiltonian.  The Journal of Physical Chemistry B, 122(13):3384–3395, 2017.
  * [19] Henry Eyring.  The activated complex in chemical reactions.  The Journal of Chemical Physics, 3(2):107–115, 1935.
  * [20] Justin Copenhaver, Adam Wasserman, and Birgit Wehefritz-Kaufmann.  Using quantum annealers to calculate ground state properties of molecules.  The Journal of Chemical Physics, 154(3):034105, 2021.
  * [21] Guglielmo Mazzola, Vadim N Smelyanskiy, and Matthias Troyer.  Quantum monte carlo tunneling from quantum chemistry to quantum annealing.  Physical Review B, 96(13):134305, 2017.
  * [22] Scott N Genin, Ilya G Ryabinkin, and Artur F Izmaylov.  Quantum chemistry on quantum annealers.  arXiv preprint arXiv:1901.04715, 2019.
  * [23] Michael Streif, Florian Neukart, and Martin Leib.  Solving quantum chemistry problems with a d-wave quantum annealer.  In International Workshop on Quantum Technology and Optimization Problems, pages 111–122. Springer, 2019.
  * [24] Maritza Hernandez and Maliheh Aramon.  Enhancing quantum annealing performance for the molecular similarity problem.  Quantum Information Processing, 16(5):133, 2017.
  * [25] Shiyu Zhou, Dmitry Green, Edward D Dahl, and Claudio Chamon.  Experimental realization of classical z 2 spin liquids in a programmable quantum device.  Physical Review B, 104(8):L081107, 2021.
  * [26] Andrew D King, Juan Carrasquilla, Jack Raymond, Isil Ozfidan, Evgeny Andriyash, Andrew Berkley, Mauricio Reis, Trevor Lanting, Richard Harris, Fabio Altomare, et al.  Observation of topological phenomena in a programmable lattice of 1,800 qubits.  Nature, 560(7719):456–460, 2018.
  * [27] Paul Kairys, Andrew D King, Isil Ozfidan, Kelly Boothby, Jack Raymond, Arnab Banerjee, and Travis S Humble.  Simulating the shastry-sutherland ising model using quantum annealing.  PRX Quantum, 1(2):020320, 2020.
  * [28] R Harris, Y Sato, AJ Berkley, M Reis, F Altomare, MH Amin, K Boothby, P Bunyk, C Deng, C Enderud, et al.  Phase transitions in a programmable quantum spin glass simulator.  Science, 361(6398):162–165, 2018.
  * [29] JE Mooij, TP Orlando, L Levitov, Lin Tian, Caspar H Van der Wal, and Seth Lloyd.  Josephson persistent-current qubit.  Science, 285(5430):1036–1039, 1999.
  * [30] TP Orlando, JE Mooij, Lin Tian, Caspar H Van Der Wal, LS Levitov, Seth Lloyd, and JJ Mazo.  Superconducting persistent-current qubit.  Physical Review B, 60(22):15398, 1999.
  * [31] John Clarke and Frank K Wilhelm.  Superconducting quantum bits.  Nature, 453(7198):1031–1042, 2008.
  * [32] Yuriy Makhlin, Gerd Schön, and Alexander Shnirman.  Quantum-state engineering with josephson-junction devices.  Reviews of modern physics, 73(2):357, 2001.
  * [33] FG Paauw, A Fedorov, CJP M Harmans, and JE Mooij.  Tuning the gap of a superconducting flux qubit.  Physical review letters, 102(9):090501, 2009.
  * [34] AO Niskanen, K Harrabi, F Yoshihara, Y Nakamura, S Lloyd, and Jaw Shen Tsai.  Quantum coherent tunable coupling of superconducting qubits.  Science, 316(5825):723–726, 2007.
  * [35] Guido Burkard, David P DiVincenzo, P Bertet, I Chiorescu, and JE Mooij.  Asymmetry and decoherence in a double-layer persistent-current qubit.  Physical Review B, 71(13):134504, 2005.
  * [36] F Yoshihara, K Harrabi, AO Niskanen, Y Nakamura, and Jaw Shen Tsai.  Decoherence of flux qubits due to 1/f flux noise.  Physical review letters, 97(16):167001, 2006.
  * [37] Jonas Bylander, Simon Gustavsson, Fei Yan, Fumiki Yoshihara, Khalil Harrabi, George Fitch, David G Cory, Yasunobu Nakamura, Jaw-Shen Tsai, and William D Oliver.  Noise spectroscopy through dynamical decoupling with a superconducting flux qubit.  Nature Physics, 7(7):565–570, 2011.
  * [38] JQ You, Xuedong Hu, S Ashhab, and Franco Nori.  Low-decoherence flux qubit.  Physical Review B, 75(14):140515, 2007.
  * [39] Fei Yan, Simon Gustavsson, Archana Kamal, Jeffrey Birenbaum, Adam P Sears, David Hover, Ted J Gudmundsen, Danna Rosenberg, Gabriel Samach, Steven Weber, et al.  The flux qubit revisited to enhance coherence and reproducibility.  Nature communications, 7(1):1–9, 2016.
  * [40] Antonio D Córcoles, Jerry M Chow, Jay M Gambetta, Chad Rigetti, James R Rozen, George A Keefe, Mary Beth Rothwell, Mark B Ketchen, and Matthias Steffen.  Protecting superconducting qubits from radiation.  Applied Physics Letters, 99(18):181906, 2011.
  * [41] Matthias Steffen, Shwetank Kumar, David P DiVincenzo, JR Rozen, George A Keefe, Mary Beth Rothwell, and Mark B Ketchen.  High-coherence hybrid superconducting qubit.  Physical review letters, 105(10):100502, 2010.
  * [42] Mustafa Bal, Chunqing Deng, Jean-Luc Orgiazzi, FR Ong, and Adrian Lupascu.  Ultrasensitive magnetic field detection using a single artificial atom.  Nature communications, 3(1):1–8, 2012.
  * [43] Hiraku Toida, Yuichiro Matsuzaki, Kosuke Kakuyanagi, Xiaobo Zhu, William J Munro, Hiroshi Yamaguchi, and Shiro Saito.  Electron paramagnetic resonance spectroscopy using a single artificial atom.  Communications Physics, 2(1):1–7, 2019.
  * [44] Rangga P Budoyo, Kosuke Kakuyanagi, Hiraku Toida, Yuichiro Matsuzaki, and Shiro Saito.  Electron spin resonance with up to 20 spin sensitivity measured using a superconducting flux qubit.  Applied Physics Letters, 116(19):194001, 2020.
  * [45] Leonid V Abdurakhimov, Imran Mahboob, Hiraku Toida, Kosuke Kakuyanagi, Yuichiro Matsuzaki, and Shiro Saito.  Driven-state relaxation of a coupled qubit-defect system in spin-locking measurements.  Physical Review B, 102(10):100502, 2020.
  * [46] Kosuke Kakuyanagi, Yuichiro Matsuzaki, Corentin Déprez, Hiraku Toida, Kouichi Semba, Hiroshi Yamaguchi, William J Munro, and Shiro Saito.  Observation of collective coupling between an engineered ensemble of macroscopic artificial atoms and a superconducting resonator.  Physical review letters, 117(21):210503, 2016.
  * [47] Mark W Johnson, Mohammad HS Amin, Suzanne Gildert, Trevor Lanting, Firas Hamze, Neil Dickson, Richard Harris, Andrew J Berkley, Jan Johansson, Paul Bunyk, et al.  Quantum annealing with manufactured spins.  Nature, 473(7346):194–198, 2011.
  * [48] Sergio Boixo, Troels F Rønnow, Sergei V Isakov, Zhihui Wang, David Wecker, Daniel A Lidar, John M Martinis, and Matthias Troyer.  Evidence for quantum annealing with more than one hundred qubits.  Nature physics, 10(3):218–224, 2014.
  * [49] Isil Ozfidan, Chunqing Deng, AY Smirnov, T Lanting, R Harris, L Swenson, J Whittaker, F Altomare, M Babcock, C Baron, et al.  Demonstration of a nonstoquastic hamiltonian in coupled superconducting flux qubits.  Physical Review Applied, 13(3):034037, 2020.
  * [50] Hongwei Chen, Xi Kong, Bo Chong, Gan Qin, Xianyi Zhou, Xinhua Peng, and Jiangfeng Du.  Experimental demonstration of a quantum annealing algorithm for the traveling salesman problem in a nuclear-magnetic-resonance quantum simulator.  Physical Review A, 83(3):032314, 2011.
  * [51] Mikio Nakahara.  Lectures on quantum computing, thermodynamics and statistical physics, volume 8.  World Scientific, 2013.
  * [52] Yuichiro Matsuzaki, Hideaki Hakoshima, Yuya Seki, and Shiro Kawabata.  Quantum annealing with capacitive-shunted flux qubits.  Japanese Journal of Applied Physics, 59(SG):SGGI06, 2020.
  * [53] Subir Sachdev.  Quantum phase transitions.  Cambridge university press, 2011.
  * [54] Jonathan R Friedman, Vijay Patel, Wei Chen, SK Tolpygo, and James E Lukens.  Quantum superposition of distinct macroscopic states.  nature, 406(6791):43–46, 2000.
  * [55] Caspar H Van Der Wal, ACJ Ter Haar, FK Wilhelm, RN Schouten, CJPM Harmans, TP Orlando, Seth Lloyd, and JE Mooij.  Quantum superposition of macroscopic persistent-current states.  Science, 290(5492):773–777, 2000.
  * [56] Takashi Imoto, Yuya Seki, Yuichiro Matsuzaki, and Shiro Kawabata.  Quantum annealing with twisted fields.  arXiv preprint arXiv:2111.15283, 2021.
  * [57] Tadashi Kadowaki and Hidetoshi Nishimori.  Greedy parameter optimization for diabatic quantum annealing.  arXiv preprint arXiv:2111.13287, 2021.
  * [58] Akhil Francis, Ephrata Zelleke, Ziyue Zhang, Alexander F Kemper, and JK Freericks.  Determining ground-state phase diagrams on quantum computers via a generalized application of adiabatic state preparation.  arXiv preprint arXiv:2112.04625, 2021.


# MRI-driven dynamo at very high magnetic Prandtl numbers

Jérôme Guilet1, Alexis Reboul-Salze1,2, Raphaël Raynaud3, Matteo Bugli1,
Basile Gallet4
1 Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191, Gif-
sur-Yvette, France
2 Max Planck Institute for Gravitational Physics (Albert Einstein Institute),
D-14476 Potsdam, Germany
3Université Paris Cité, Université Paris-Saclay, CEA, CNRS, AIM, F-91191, Gif-
sur-Yvette, France
4 Université Paris-Saclay, CNRS, CEA, Service de Physique de l’Etat Condensé,
91191, Gif-sur-Yvette, France

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The dynamo driven by the magnetorotational instability (MRI) is believed to
play an important role in the dynamics of accretion discs and may also explain
the origin of the extreme magnetic fields present in magnetars. Its saturation
level is an important open question known to be particularly sensitive to the
diffusive processes through the magnetic Prandtl number Pm (the ratio of
viscosity to resistivity). Despite its relevance to proto-neutron stars and
neutron star merger remnants, the numerically challenging regime of high Pm is
still largely unknown. Using zero-net flux shearing box simulations in the
incompressible approximation, we studied MRI-driven dynamos at unprecedentedly
high values of Pm reaching 256. The simulations show that the stress and
turbulent energies are proportional to Pm up to moderately high values
($\mathrm{Pm}\sim 50$). At higher Pm, they transition to a new regime
consistent with a plateau independent of Pm for $\rm Pm\gtrsim 100$. This
trend is independent of the Reynolds number, which may suggest an asymptotic
regime where the energy injection and dissipation are independent of the
diffusive processes. Interestingly, large values of Pm not only lead to
intense small-scale magnetic fields but also to a more efficient dynamo at the
largest scales of the box.

###### keywords:

stars: magnetars – neutron star mergers – supernovae: general – MHD –
instabilities – magnetic fields

††pubyear: 2015††pagerange: MRI-driven dynamo at very high magnetic Prandtl
numbers–References

##  1 Introduction

The magnetorotational instability (MRI) is believed to play a crucial role in
a large number of astrophysical objects. This includes accretion discs around
a variety of objects (Balbus & Hawley, 1998), neutron star mergers (e.g.
Siegel et al., 2013; Kiuchi et al., 2014, 2018; Guilet et al., 2017), stellar
mergers (e.g. Schneider et al., 2019) and core-collapse supernovae (e.g.
Akiyama et al., 2003; Obergaulinger et al., 2009; Masada et al., 2007; Guilet
& Müller, 2015; Mösta et al., 2015; Reboul-Salze et al., 2021a). In neutron
star mergers and core-collapse supernovae, the MRI may generate extreme
magnetic fields reaching up to $10^{15}-10^{16}\,\mathrm{G}$, potentially
leading to the formation of a magnetar (Reboul-Salze et al., 2021a, b). In
combination with fast rotation, such extreme magnetic fields can trigger
powerful magnetorotational explosions111Note that with _magnetorotational
explosion_ we refer to the physical explosion mechanism of massive magnetised
stars in rotation, not to be confounded with the MRI. (e.g. Takiwaki et al.,
2009; Kuroda et al., 2020; Bugli et al., 2020; Bugli et al., 2021). This so-
called millisecond magnetar scenario may provide an explanation for
outstanding explosions such as long gamma-ray bursts, hypernovae and
superluminous supernovae (e.g. Woosley, 2010; Metzger et al., 2011). The
formation of a magnetar after a neutron star merger has also been proposed as
a central engine of short gamma-ray bursts (e.g. Mösta et al., 2020).

The origin of the extreme magnetic fields of magnetars is still an important
open question as several scenarios have been proposed, such as a convective
dynamo in a fast rotating proto-neutron star (PNS) (Thompson & Duncan, 1993;
Raynaud et al., 2020, 2022), the Tayler-Spruit dynamo following fallback
(Barrere et al, in prep), amplification in main sequence stellar mergers
(Schneider et al., 2019) or the fossil field scenario (Ferrario &
Wickramasinghe, 2006). To shed light on this question, it is important to
assess the efficiency of each of the dynamo mechanisms in the conditions
specific to a PNS. The impact of different ingredients on the MRI has been
studied in recent years: the shear parameter (Masada et al., 2012), neutrino
viscosity and drag (Guilet et al., 2015), stable stratification (Guilet &
Müller, 2015; Reboul-Salze et al., 2021b) and spherical geometry (Reboul-Salze
et al., 2021a). A remaining open question is the dependence on diffusion
coefficients. Local disc simulations have shown that the efficiency of the MRI
is strongly correlated to the magnetic Prandtl number $\mathrm{Pm}$ (the ratio
of viscosity to resistivity) when $\mathrm{Pm}\sim 0.1-16$ (Fromang et al.,
2007; Lesur & Longaretti, 2007; Simon & Hawley, 2009; Longaretti & Lesur,
2010; Shi et al., 2016; Potter & Balbus, 2017) and plateaus at
$\mathrm{Pm}<0.1$ with an imposed external magnetic field (Meheut et al.,
2015). Low and high Pm regimes are computationally challenging because a high
numerical resolution is needed to resolve the small viscous (low Pm) or
resistive (high Pm) scale. The regime of low Pm has attracted particular
attention because it is relevant to most regions of accretion discs as well as
to liquid metals (laboratory experiments, Earth core…), but the opposite
regime of large $\mathrm{Pm}$ has not been targeted specifically by previous
studies of the MRI. This regime is relevant to PNSs (Thompson & Duncan, 1993;
Guilet et al., 2015; Lander, 2021), neutron star merger remnants (Rossi et
al., 2008), which exhibit physical conditions similar to proto-neutron stars
(Guilet et al., 2017), interstellar and intergalactic media (Schekochihin et
al., 2004), and to the inner parts of some accretion discs (Balbus & Henri,
2008; Potter & Balbus, 2014, 2017). In the latter case, the Pm-dependence of
the MRI may drive an instability leading to variability in the accretion rate
(Potter & Balbus, 2014, 2017; Kawanaka & Masada, 2019).

Simulations of the MRI in a PNS including explicit diffusion coefficients are
able to describe a realistic parameter regime for the (large) viscosity and
thermal diffusion induced by neutrinos (Guilet et al., 2015; Reboul-Salze et
al., 2021a, b). By contrast, the physical resistivity in a PNS is much smaller
than the values that can be afforded in these simulations. This translates
into a large discrepancy between the huge values of $\mathrm{Pm}\sim 10^{13}$
relevant to a PNS and the modest values $\rm Pm=4-16$ of these simulations.
The regime of high magnetic Prandtl numbers is very challenging for numerical
simulations because a high resolution is required to resolve the resistive
length scale. Local models are best suited to tackle this problem owing to
their simplicity and comparatively low computational cost. We probe
unprecedentedly large $\mathrm{Pm}$ values (up to 256), using a zero-net flux
shearing-box model inspired by a comparison between local and spherical models
of the MRI (Reboul-Salze et al., 2021a).

##  2 Numerical setup

Our simulations are designed to represent a small region in the equatorial
plane of a fast rotating PNS. The local dynamics is described in the framework
of a Cartesian shearing box (e.g. Goldreich & Lynden-Bell, 1965). The
coordinates $x$, $y$, and $z$ represent the radial, azimuthal and vertical
directions, respectively, and the corresponding unit vectors are $\bm{e_{x}}$,
$\bm{e_{y}}$, and $\bm{e_{z}}$. The angular frequency vector points in the $z$
direction $\bm{\Omega}=\Omega\,\bm{e_{z}}$, while gravity and shear are in the
$x$ direction. Neutrinos are assumed to be in the diffusive regime such that
their effects on the dynamics can be described by a viscosity (Guilet et al.,
2015). We used the incompressible approximation for the following reasons.
Soundproof approximations are well justified for the study of the MRI in a PNS
because the fluid velocity and the Alfvén speed are small compared to the
sound speed ($v/c_{s}<v_{A}/c_{s}\lesssim 10^{-2}$, see e.g. Guilet et al.,
2015; Reboul-Salze et al., 2021a, b). Limits to the applicability of the
incompressible approximation in a PNS come mostly from the density gradient
and the buoyancy driven by entropy and composition gradients. Neglecting the
density gradient is an essential part of the local approximation, which is
necessary in order to reach high magnetic Prandl numbers. On the other hand,
while buoyancy could have been included with the Boussinesq approximation, it
was shown by Reboul-Salze et al. (2021b) to have relatively minor effects on
the MRI in a PNS except for a small region near the equator. For this reason
and for the sake of simplicity we adopted the incompressible approximation.
The incompressible MHD equations in the shearing box approximation read

$\displaystyle\begin{split}\partial_{t}\bm{v}+\bm{v}\cdot\bm{\nabla}\bm{v}&=-\frac{1}{\rho_{0}}\bm{\nabla}\Pi+\frac{1}{\mu_{0}\rho_{0}}(\bm{\nabla}\times\bm{B})\times\bm{B}\\\
&\quad+q\Omega x\partial_{y}\bm{v}+q\Omega
v_{x}\bm{e_{y}}-2\bm{\Omega}\times\bm{v}+\nu\bm{\Delta}\bm{v},\end{split}$ (1)
$\displaystyle\partial_{t}\bm{B}$
$\displaystyle=\bm{\nabla}\times(\bm{v}\times\bm{B})+q\Omega
x\partial_{y}\bm{B}-q\Omega B_{x}\bm{e_{y}}+\eta\bm{\Delta{B}},$ (2)
$\displaystyle\nabla\cdot\bm{v}$ $\displaystyle=0,$ (3)
$\displaystyle\nabla\cdot\bm{B}$ $\displaystyle=0,$ (4)

where $\bm{B}$ is the magnetic field, $\rho_{0}$ the density, $\nu$ the
kinematic viscosity and $\eta$ the magnetic diffusivity. The gradient of the
pressure perturbation $\bm{\nabla}\Pi$ is obtained from the constraint of a
divergence-free flow field (Eq. 3). $\bm{v}$ is the velocity fluctuation with
respect to the mean shear profile $\bm{U}=-q\Omega x\bm{e_{y}}$, where the
shear parameter $q\equiv-{\rm d}\log\Omega/{\rm d}\log r$ is assumed to have a
sub-Keplerian value of $q=0.8$ in all simulations (such value is relevant in
particular for proto-neutron stars, e.g. Reboul-Salze et al., 2021a; Bugli et
al., 2020). The box dimensions $(L_{x},L_{y},L_{z})=(1,3,3)L$ are chosen in
light of the comparison with global simulations performed by Reboul-Salze et
al. (2021a).

Table 1: Overview of the numerical simulations. The second and third columns show the dimensionless control parameters that have been varied in this study: the Reynolds number $\mathrm{Re}\equiv L^{2}\Omega/\nu$ and the magnetic Prandtl number $\mathrm{Pm}\equiv\nu/\eta$. The time and volume averages of the kinetic ($E_{\rm kin}\equiv v^{2}/2$) and magnetic ($E_{\rm mag}\equiv B^{2}/2\mu_{0}$) energy densities and the total stress ($\alpha\equiv\rho_{0}v_{x}v_{y}-B_{x}B_{y}/\mu_{0}$) are in units of $\rho_{0}L^{2}\Omega^{2}$, with $L$ the radial size of the box with dimensions $(L_{x},L_{y},L_{z})=(1,3,3)L$. The resistive length scale $l_{\rm res}$ is defined as the wavelength such that $50\%$ of the resistive dissipation takes place at shorter/longer scales. The ratio of the resistive scale to the radial size of the grid cells $\Delta x$ is used to check that the simulation is sufficiently resolved. Name | ${\rm Re}$ | ${\rm Pm}$ | $E_{\rm kin}$ | $E_{\rm mag}$ | $E_{\rm mag}^{k<3\pi/L}$ | $E_{\rm mag}^{k<\pi/L}$ | $\alpha$ | $n_{x}$ | $n_{y}$ | $n_{z}$ | $l_{\rm res}/\Delta x$
---|---|---|---|---|---|---|---|---|---|---|---
Re400Pm24n64 | $400$ | $24$ | $2.89\times 10^{-3}$ | $4.87\times 10^{-2}$ | $9.14\times 10^{-3}$ | $3.58\times 10^{-3}$ | $1.90\times 10^{-2}$ | $64$ | $96$ | $192$ | $8.4$
Re400Pm24n96 | $400$ | $24$ | $3.85\times 10^{-3}$ | $6.27\times 10^{-2}$ | $1.42\times 10^{-2}$ | $6.48\times 10^{-3}$ | $2.45\times 10^{-2}$ | $96$ | $128$ | $256$ | $12.1$
Re400Pm32n96 | $400$ | $32$ | $4.14\times 10^{-3}$ | $7.01\times 10^{-2}$ | $1.38\times 10^{-2}$ | $5.91\times 10^{-3}$ | $2.75\times 10^{-2}$ | $96$ | $128$ | $256$ | $10.4$
Re400Pm32n96hyz | $400$ | $32$ | $4.41\times 10^{-3}$ | $7.51\times 10^{-2}$ | $1.61\times 10^{-2}$ | $7.27\times 10^{-3}$ | $2.90\times 10^{-2}$ | $96$ | $256$ | $512$ | $10.3$
Re400Pm32n128 | $400$ | $32$ | $4.62\times 10^{-3}$ | $7.70\times 10^{-2}$ | $1.64\times 10^{-2}$ | $7.77\times 10^{-3}$ | $3.01\times 10^{-2}$ | $128$ | $192$ | $384$ | $13.6$
Re400Pm48n96 | $400$ | $48$ | $5.57\times 10^{-3}$ | $1.03\times 10^{-1}$ | $2.68\times 10^{-2}$ | $1.79\times 10^{-2}$ | $3.70\times 10^{-2}$ | $96$ | $128$ | $256$ | $8.4$
Re400Pm48n128 | $400$ | $48$ | $5.40\times 10^{-3}$ | $9.77\times 10^{-2}$ | $1.88\times 10^{-2}$ | $9.57\times 10^{-3}$ | $3.75\times 10^{-2}$ | $128$ | $192$ | $384$ | $11.0$
Re400Pm64n96 | $400$ | $64$ | $6.05\times 10^{-3}$ | $1.10\times 10^{-1}$ | $2.46\times 10^{-2}$ | $1.30\times 10^{-2}$ | $4.03\times 10^{-2}$ | $96$ | $128$ | $256$ | $7.5$
Re400Pm64n128 | $400$ | $64$ | $5.98\times 10^{-3}$ | $1.17\times 10^{-1}$ | $2.54\times 10^{-2}$ | $1.51\times 10^{-2}$ | $4.28\times 10^{-2}$ | $128$ | $192$ | $384$ | $9.5$
Re400Pm80n128 | $400$ | $80$ | $6.55\times 10^{-3}$ | $1.33\times 10^{-1}$ | $2.87\times 10^{-2}$ | $1.69\times 10^{-2}$ | $4.75\times 10^{-2}$ | $128$ | $192$ | $384$ | $8.5$
Re400Pm80n192 | $400$ | $80$ | $6.24\times 10^{-3}$ | $1.35\times 10^{-1}$ | $3.00\times 10^{-2}$ | $1.70\times 10^{-2}$ | $4.69\times 10^{-2}$ | $192$ | $256$ | $512$ | $12.6$
Re400Pm96n128 | $400$ | $96$ | $6.75\times 10^{-3}$ | $1.41\times 10^{-1}$ | $3.15\times 10^{-2}$ | $1.76\times 10^{-2}$ | $4.93\times 10^{-2}$ | $128$ | $192$ | $384$ | $7.9$
Re400Pm96n192 | $400$ | $96$ | $6.69\times 10^{-3}$ | $1.45\times 10^{-1}$ | $2.58\times 10^{-2}$ | $1.06\times 10^{-2}$ | $5.23\times 10^{-2}$ | $192$ | $256$ | $512$ | $11.5$
Re400Pm128n128 | $400$ | $128$ | $6.76\times 10^{-3}$ | $1.58\times 10^{-1}$ | $4.42\times 10^{-2}$ | $2.94\times 10^{-2}$ | $4.98\times 10^{-2}$ | $128$ | $192$ | $384$ | $7.1$
Re400Pm128n192 | $400$ | $128$ | $7.04\times 10^{-3}$ | $1.61\times 10^{-1}$ | $3.28\times 10^{-2}$ | $2.14\times 10^{-2}$ | $5.56\times 10^{-2}$ | $192$ | $256$ | $512$ | $10.0$
Re400Pm192n256 | $400$ | $192$ | $7.56\times 10^{-3}$ | $2.03\times 10^{-1}$ | $3.69\times 10^{-2}$ | $1.80\times 10^{-2}$ | $6.69\times 10^{-2}$ | $256$ | $384$ | $768$ | $10.8$
Re400Pm256n256 | $400$ | $256$ | $7.23\times 10^{-3}$ | $2.13\times 10^{-1}$ | $4.93\times 10^{-2}$ | $3.59\times 10^{-2}$ | $6.57\times 10^{-2}$ | $256$ | $384$ | $768$ | $9.6$
Re800Pm10n96 | $800$ | $10$ | $2.31\times 10^{-3}$ | $2.87\times 10^{-2}$ | $5.99\times 10^{-3}$ | $2.74\times 10^{-3}$ | $1.17\times 10^{-2}$ | $96$ | $128$ | $256$ | $13.0$
Re800Pm12.5n96 | $800$ | $12$ | $2.43\times 10^{-3}$ | $3.10\times 10^{-2}$ | $5.26\times 10^{-3}$ | $2.28\times 10^{-3}$ | $1.30\times 10^{-2}$ | $96$ | $128$ | $256$ | $11.5$
Re800Pm16n128 | $800$ | $16$ | $3.46\times 10^{-3}$ | $4.35\times 10^{-2}$ | $7.66\times 10^{-3}$ | $3.52\times 10^{-3}$ | $1.86\times 10^{-2}$ | $128$ | $192$ | $384$ | $13.0$
Re800Pm22.5n128 | $800$ | $22$ | $4.79\times 10^{-3}$ | $6.06\times 10^{-2}$ | $1.05\times 10^{-2}$ | $4.64\times 10^{-3}$ | $2.62\times 10^{-2}$ | $128$ | $192$ | $384$ | $10.5$
Re800Pm32n128 | $800$ | $32$ | $6.26\times 10^{-3}$ | $8.03\times 10^{-2}$ | $1.47\times 10^{-2}$ | $7.45\times 10^{-3}$ | $3.42\times 10^{-2}$ | $128$ | $192$ | $384$ | $8.7$
Re800Pm48n192 | $800$ | $48$ | $7.18\times 10^{-3}$ | $1.01\times 10^{-1}$ | $1.84\times 10^{-2}$ | $8.27\times 10^{-3}$ | $4.20\times 10^{-2}$ | $192$ | $256$ | $512$ | $10.4$
Re800Pm64n256 | $800$ | $64$ | $8.17\times 10^{-3}$ | $1.22\times 10^{-1}$ | $2.26\times 10^{-2}$ | $1.35\times 10^{-2}$ | $4.98\times 10^{-2}$ | $256$ | $384$ | $768$ | $11.8$
Re800Pm80n256 | $800$ | $80$ | $8.10\times 10^{-3}$ | $1.25\times 10^{-1}$ | $1.35\times 10^{-2}$ | $3.82\times 10^{-3}$ | $5.12\times 10^{-2}$ | $256$ | $384$ | $768$ | $10.5$
Re800Pm96n256 | $800$ | $96$ | $8.69\times 10^{-3}$ | $1.47\times 10^{-1}$ | $2.83\times 10^{-2}$ | $1.51\times 10^{-2}$ | $5.52\times 10^{-2}$ | $256$ | $384$ | $768$ | $9.7$
Re1600Pm8n96 | $1600$ | $8$ | $1.03\times 10^{-3}$ | $1.25\times 10^{-2}$ | $1.71\times 10^{-3}$ | $7.30\times 10^{-4}$ | $4.93\times 10^{-3}$ | $96$ | $128$ | $256$ | $10.4$
Re1600Pm10n128 | $1600$ | $10$ | $1.37\times 10^{-3}$ | $1.60\times 10^{-2}$ | $1.89\times 10^{-3}$ | $7.83\times 10^{-4}$ | $6.74\times 10^{-3}$ | $128$ | $192$ | $384$ | $11.8$
Re1600Pm12.5n128 | $1600$ | $12$ | $3.26\times 10^{-3}$ | $3.34\times 10^{-2}$ | $5.39\times 10^{-3}$ | $2.76\times 10^{-3}$ | $1.49\times 10^{-2}$ | $128$ | $192$ | $384$ | $9.7$
Re1600Pm16n128 | $1600$ | $16$ | $3.98\times 10^{-3}$ | $4.06\times 10^{-2}$ | $6.36\times 10^{-3}$ | $3.17\times 10^{-3}$ | $1.83\times 10^{-2}$ | $128$ | $192$ | $384$ | $8.5$
Re1600Pm24n192 | $1600$ | $24$ | $6.14\times 10^{-3}$ | $6.39\times 10^{-2}$ | $1.04\times 10^{-2}$ | $4.93\times 10^{-3}$ | $2.87\times 10^{-2}$ | $192$ | $256$ | $512$ | $9.8$
Re1600Pm32n256 | $1600$ | $32$ | $5.40\times 10^{-3}$ | $6.27\times 10^{-2}$ | $7.47\times 10^{-3}$ | $2.68\times 10^{-3}$ | $2.76\times 10^{-2}$ | $256$ | $384$ | $768$ | $11.2$
Re1600Pm48n256 | $1600$ | $48$ | $8.84\times 10^{-3}$ | $1.10\times 10^{-1}$ | $2.50\times 10^{-2}$ | $1.51\times 10^{-2}$ | $4.32\times 10^{-2}$ | $256$ | $384$ | $768$ | $8.9$

Figure 1: Top panel: Time evolution of the magnetic energy for a subset of
simulations with $\mathrm{Re}=400$ and $\mathrm{Pm}\in[24,256]$. Lower panels:
Time and volume averaged magnetic energy ($E_{\rm mag}\equiv B^{2}/2\mu_{0}$),
kinetic energy ($E_{\rm kin}\equiv v^{2}/2$) and stress
($\alpha\equiv\rho_{0}v_{x}v_{y}-B_{x}B_{y}/\mu_{0}$) as functions of the
magnetic Prandtl number for three different values of the Reynolds number
$\mathrm{Re}=400$ (red circles), $\mathrm{Re}=800$ (blue crosses) and
$\mathrm{Re}=1600$ (black triangles). The error bars show the standard
deviation, which probably overestimates the actual error on the averages. The
dotted line shows a linear fit valid at moderate $\mathrm{Pm}$, while the
dashed line shows the fit with Eq. 5. The stress and energies are in units of
$\rho_{0}L^{2}\Omega^{2}$ (with $L$ the radial size of the box).

Figure 2: 3D rendering of the azimuthal component of the velocity field (left
panels) and of the magnetic field (middle panels). The right panels represent
the time-averaged spectra of the force balance. The upper row corresponds to a
simulation with a moderate value of $\mathrm{Pm}=24$, while the lower row
corresponds to a large value of $\mathrm{Pm}=256$. Both simulations share the
same Reynolds number $\mathrm{Re}=400$. In the force balance spectrum, only
the divergence-free part of the forces is included as it is appropriate in the
incompressible approximation.

The shear rate and the box aspect ratio being fixed, this setup is governed by
only two dimensionless numbers: the Reynolds number $\mathrm{Re}\equiv
L^{2}\Omega/\nu$ and the magnetic Reynolds number $\mathrm{Rm}\equiv
L^{2}\Omega/\eta$. We obtained self-sustained MRI-driven turbulence in 30
simulations with $\mathrm{Rm}$ ranging from $8\,000$ to $102\,400$ and with
three different values of the Reynolds number ($\mathrm{Re}=400$, $800$ and
$1\,600$). In this set of simulations, the magnetic Prandtl number
$\mathrm{Pm}\equiv\nu/\eta=\mathrm{Rm/Re}$ lies in the range
$\mathrm{Pm}=8-256$. Note that for each Reynolds number, we ran simulations
with lower magnetic Reynolds number where the dynamo was not self-sustained,
which are therefore not included in the present analysis.

Our simulations are initialised with a random superposition of large-scale
magnetic modes, with exactly zero net magnetic flux. We checked in a few cases
that, after a transient phase, the turbulent state was statistically
independent of the initial conditions, provided that the initial magnetic
field was strong enough to initiate an MRI-driven dynamo. The simulations were
run for a duration varying between $1500$ and $3000\,\Omega^{-1}$, which is
longer than the typical timescale of the fluctuations (see upper panel of Fig.
1). The time-averaged results shown in this paper are performed at times
$t>500\,\Omega^{-1}$ in order to exclude any initial transient behaviour.

Throughout the paper, our results are normalised using the radial size of the
domain $L$, the angular frequency $\Omega$, and the density $\rho_{0}$. With
parameters typical of a proto-neutron star $L=10\,{\rm km}$,
$\Omega=10^{3}\,{\rm s^{-1}}$, and $\rho_{0}=2\times 10^{13}\,{\rm
g\,cm^{-3}}$, time would therefore be measured in units of $1\,{\rm ms}$,
velocity in units of $10^{9}\,{\rm cm\,s^{-1}}$, the magnetic field in units
of $1.6\times 10^{16}\,{\rm G}$, and the energy density in units of $2\times
10^{31}\,{\rm erg\,cm^{-3}}$. Our range of Reynolds numbers translates into a
viscosity ranging from $\nu=6\times 10^{11}$ to $2.5\times 10^{12}\,{\rm
cm^{2}\,s^{-1}}$, which is comparable to estimates of the neutrino viscosity
in the outer parts of a PNS (Guilet et al., 2015).

###  2.1 Numerical methods and convergence tests

In order to solve the incompressible MHD equations (1)–(4), we use the pseudo-
spectral code snoopy (Lesur & Longaretti, 2005, 2007), which has been used in
numerous studies of the MRI (e.g. Lesur & Longaretti, 2011; Guilet et al.,
2015; Walker & Boldyrev, 2017). Our simulations were performed using a grid
resolution varying from $(n_{x},n_{y},n_{z})=(64,128,192)$ to
$(n_{x},n_{y},n_{z})=(256,512,768)$. In order to ensure that the resistive
scale is resolved, the number of grid points is increased as Pm is increased.
In all simulations presented in the figures of this paper (except for some of
the low resolution simulations performed for the convergence study), the
resistive scale is resolved by at least 8.5 grid cells (see Table 1).

In order to check the convergence of our results, we ran a subset of 8
additional simulations at lower resolution for $\mathrm{Re}=400$,
$\mathrm{Pm}=24-128$. We define a convergence criterion $l_{\rm res}/\Delta x$
as the ratio of the resistive scale (measured such that half of the resistive
dissipation takes place at higher/lower scales) to our grid scale. The precise
threshold needed for convergence is expected to depend on the numerical
scheme222Pseudo-spectral methods like the one we use are known for their low
dissipation and for necessitating fewer grid points than non-spectral grid-
based methods (e.g. Fromang et al., 2007) and potentially on the physical
problem considered through the width of the dissipation peak. As a
consequence, it is necessary to calibrate it through dedicated convergence
tests as it is done here. In all our simulations, the wavenumber at which the
dissipation peak has decreased by a factor 2 is roughly twice as large as the
wavenumber of the dissipation maximum independently of Pm. This suggests that
a fixed threshold in $l_{\rm res}/\Delta x$ can be a meaningful measure of
convergence in the parameter space explored. When comparing the low resolution
runs with their higher resolution counterparts (Table 1 and Fig. 1), we
observe that the stress and turbulent energies of the low resolution
simulations have a tendency to be slightly smaller (with the exception of a
few cases) but they are all consistent within the statistical error bars
defined with the standard deviation. Most simulations show differences of at
most a few percent, while the least resolved test (Re400Pm128n128 with $l_{\rm
res}/\Delta x=7.1$, which we would deem slightly under-resolved) has a larger
difference of about $10\%$ in the stress compared to its high resolution
counterpart (however, the difference is smaller for the turbulent energies).
We stress that all our runs are better resolved than this under-resolved test
according to our resolution criterion ($l_{\rm res}/\Delta x>8.5$ for all
simulations, and the three highest Pm runs have $l_{\rm res}/\Delta x>9.6$).
Our convergence tests therefore suggest that all the simulations included in
the analysis give reliable results on the average energies and stress with
errors of at most a few percent due to numerical artefacts. Such errors are
within the statistical errors shown in Fig. 1 and do not compromise our main
conclusions.

In some of the spectra, a numerical artefact can be discerned at wavenumbers
larger than the dissipation peak. We checked that, in all the simulations
included in the analysis, the resistive dissipation rate at this artefact is
smaller than the resistive dissipation peak by a factor of at least 10 (e.g.
40 and 15 for the simulations at $\mathrm{Pm}=192$ and 256, respectively). In
four of the low resolution tests this artefact is more pronounced and leads to
a breakdown of this criterion (Re400Pm48n96, Re400Pm64n96, Re400Pm96n128,
Re400Pm128n128 with a ratio of the artefact to peak resistive dissipation rate
of 6, 3, 7 and 4, respectively), which allows to check the potential influence
of such artefact on the overall dynamics and energetics. Therefore, the
convergence tests discussed above give confidence that the influence of the
sub-resistive-scale artefact on the averaged energies and stress is small (at
most a few percent) in the simulations included in the analysis.

Finally, note that the resolution in the azimuthal direction is twice as low
as in the radial or vertical direction, because the structures are more
elongated in the azimuthal direction due to the shear (Fig. 2). We checked in
a few cases that our results are not affected by the lower azimuthal
resolution.

##  3 Results

Figure 3: Time-averaged spectra of a subset of simulations with
$\mathrm{Re}=400$ and varying $\mathrm{Pm}$. Upper panel: kinetic (dotted
lines) and magnetic (solid lines) energy spectra. Lower panel: viscous (dotted
lines) and resistive (solid lines) dissipation rate spectra. All spectra are
multiplied by $k$ so that the typical scale containing most of the energy or
dissipation is clearly visualised as the maximum of the curve. Figure 4:
Typical wavenumbers of the resistive dissipation (red symbols, defined such
that half of the resistive dissipation takes place at higher/lower wavenumber)
and of the kinetic energy (blue, defined such that half of the kinetic energy
is located at higher/lower wavenumbers). The wavenumbers have been normalised
by $\sqrt{\mathrm{Re}}$ so that the dependence on the Reynolds number is
approximately scaled out. The symbol shapes correspond to different Reynolds
numbers as in Fig. 1: $\mathrm{Re}=400$ (circles), $\mathrm{Re}=800$ (crosses)
and $\mathrm{Re}=1600$ (triangles). The low resolution tests at
$\mathrm{Re}=400$ are represented with the black plus signs. The dashed line
represent a scaling of the resistive wavenumber as
$\propto\sqrt{\mathrm{Pm}}$. Figure 5: Magnetic energy contained in the
largest scales of the box as a function of the magnetic Prandtl number. Blue
and red symbols correspond to the energy at wavenumbers smaller than $\pi/L$
and $3\pi/L$ respectively, while black symbols show the total magnetic energy.
The symbol shapes correspond to different Reynolds numbers as in Fig. 1:
$\mathrm{Re}=400$ (circles), $\mathrm{Re}=800$ (crosses) and
$\mathrm{Re}=1600$ (triangles). The low resolution tests at $\mathrm{Re}=400$
are represented with the plus signs. The dashed lines represent the fit of the
total magnetic energy (black), and a fraction of $20\%$ (red) and $10\%$
(blue) of this fit, respectively.

The time-evolution of the magnetic energy for a subset of our simulations with
varying $\mathrm{Pm}$ shows a clear trend of increasing magnetic energy for
larger Pm (Fig. 1, upper panel). This trend nevertheless stops towards the
largest values of $\mathrm{Pm}$, as the two simulations $\mathrm{Pm}=192$ and
$\mathrm{Pm}=256$ appear to converge on the same value. This conclusion is
confirmed quantitatively by the lower panels of Fig. 1, showing the magnetic
energy, kinetic energy and total stress time-averaged over the quasi-
stationary turbulent phase. All three quantities increase approximately
linearly with Pm up to $\mathrm{Pm}\lesssim 50$. Interestingly, at still
larger values of $\mathrm{Pm}$, they transition to a less steep dependence and
are consistent with a plateau for $\mathrm{Pm}\gtrsim 100$ (or slightly more
for the magnetic energy). It is noteworthy that all three sets of simulations
with different Reynolds numbers fall on the same universal curves as a
function of $\mathrm{Pm}$. This dependence can be well fitted by a functional
form

$A(\mathrm{Pm})=A^{\infty}\frac{\mathrm{Pm}/\mathrm{Pm_{c}}}{\sqrt{1+(\mathrm{Pm}/\mathrm{Pm_{c}})^{2}}},$
(5)

as represented by dotted lines in Fig. 1. The best fit parameters give
critical magnetic Prandtl numbers for the start of the plateau of
$\mathrm{Pm_{c}}=91$, $37$ and $61$ for the magnetic energy, kinetic energy
and the stress respectively. The asymptotic values at large $\mathrm{Pm}$ are
$E^{\infty}_{\mathrm{mag}}=0.21$, $E^{\infty}_{\mathrm{kin}}=7.8\times
10^{-3}$ and $\alpha^{\infty}=6.4\times 10^{-2}$. Because of the higher value
of $\mathrm{Pm_{c}}$ for the magnetic energy, Pm-independent plateaus are more
clearly visible for the stress and kinetic energy. As a consequence, we cannot
rule out a slight increase of the magnetic energy with Pm at high Pm, either
logarithmically or in a very weak power-law (similarly to Alexakis, 2011).

A striking feature of Fig. 1 is that, for a given Pm, the stress and magnetic
energy are independent of the Reynolds number and that the kinetic energy has
only a weak dependence with Re. This may seem surprising given the relatively
low values of Re considered, which could have suggested that the viscosity
would impact the MRI. The impact of viscosity on the MRI linear growth is
controlled by the viscous Elsasser number $E_{\nu}\equiv v_{A}^{2}/\nu\Omega$,
with strong viscous effects for $E_{\nu}<1$. Although there is no linear phase
of the MRI for the zero-net flux case, we can try to estimate the impact of
viscosity on the saturated state with the Elsasser number computed with the
turbulent magnetic field. In our simulations, $E_{\nu}$ is found to range from
40 to 350. Such relatively high values may give a hint as to why the stress
and energies are roughly independent of the Reynolds number.

In MRI-driven turbulence, the turbulent magnetic energy is usually
significantly larger than the kinetic energy, and similarly the Maxwell stress
is larger than the Reynolds stress. In a Keplerian disk, the magnetic to
kinetic ratios are often in the range 3 to 6, though they may depend on the
precise setup. In our simulations, the kinetic to magnetic ratios are higher
(between 10 and 20), which is likely due to the sub-Keplerian shear (see
Figure 4 of Pessah et al., 2006).

The structures of the velocity and magnetic field are illustrated in Fig. 2
with snapshots representative of the regimes of moderate $\mathrm{Pm}$ (upper
row) and large $\mathrm{Pm}$ (lower row). The first striking feature is the
highly entangled small-scale structure of the magnetic field in the high
$\mathrm{Pm}$ regime. By contrast, in both the moderate and high $\mathrm{Pm}$
regimes the velocity field is dominated by structures at much larger scales
comparable to the radial size of the box. One can nonetheless note the
appearance of subdominant small-scales structures in the velocity field in the
high $\mathrm{Pm}$ regime.

A more quantitative description of the velocity and magnetic field structures
is provided by their energy and dissipation spectra shown in Fig. 3 for
varying values of $\mathrm{Pm}$ at $\mathrm{Re}=400$. As $\mathrm{Pm}$
increases, both the magnetic and kinetic energy spectra initially increase at
all wavenumbers333The only exception being the intermediate range of
wavenumbers for $\mathrm{Pm}=256$.. The kinetic energy peaks at a small
wavenumber $k\sim 10$ independently of Pm, while the peak of the magnetic
energy spectrum increases with Pm from $k\sim 30$ to $k\sim 70$. Although the
peak of the kinetic energy does not change with Pm, a significant tail at
higher wavenumbers appears in the high Pm regime. This can be understood as a
result of a balance between the Lorentz and viscous forces, which drives fluid
motions (Fig. 2, right panels). The viscous dissipation in this tail broadens
significantly the peak of viscous dissipation across larger wavenumbers (lower
panel of Fig. 3). Nevertheless, the majority of the dissipation is due to the
resistivity acting at still higher wavenumbers, with a peak of the dissipation
rate ranging from $k\sim 60$ to $k\sim 200$ as $\mathrm{Pm}$ increases.

Fig. 4 shows the Pm dependence of the characteristic scales of resistive
dissipation (red symbols) and kinetic energy (blue symbols), where the
dependence on the Reynolds number has been approximately scaled out by
dividing the wavenumbers by $\sqrt{\mathrm{Re}}$. The characteristic resistive
wavenumber is proportional to $\sqrt{\mathrm{Pm}}$ while the kinetic energy
scale is approximately independent of Pm. The scaling of the resistive
wavenumber may be explained by equating the shearing term that generates
magnetic field at a rate $q\Omega$ and the resistive decay rate $\eta k^{2}$,
which leads to an estimate of the resistive wavenumber
$kL\sim\sqrt{\mathrm{Rm}}$ or equivalently
$kL/\sqrt{\mathrm{Re}}\sim\sqrt{\mathrm{Pm}}$.

The analysis of the spectra and characteristic scales may suggest the
following interpretation for the transition to a plateau independent of Pm. As
Pm is increased (with fixed Re), the resistive scale becomes shorter and
shorter as $\propto\sqrt{\mathrm{Pm}}$ while the velocity scale stays constant
(see Fig. 3 and Fig. 4). At very high $\mathrm{Pm}>100$, the resistive scale
is at least an order of magnitude shorter than the velocity scale. The
dynamics may then become independent of the resistivity because the scale
separation between the velocity and resistive scales prevents an efficient
feedback of the resistive scales on the much larger scales of the flow. This
situation would then be analogous to hydrodynamic turbulence, whose large-
scale dynamics and overall energy budget become independent of the Reynolds
number when a sufficient scale separation prevents feedback of the viscous
scales on the injection scales.

We stress again that increasing $\mathrm{Pm}$ not only leads to more intense
smaller-scale magnetic fields but also to stronger large-scale magnetic
fields. Fig. 5 shows that the magnetic energy contained in large-scale
structures increases with $\mathrm{Pm}$ proportionally to the total magnetic
energy. Structures at some of the largest allowed scales of the box with
$k<\pi/L$ and $k<3\pi/L$ represent respectively about $10\%$ and $20\%$ of the
total magnetic energy independently of $\mathrm{Pm}$.

##  4 Conclusion

We performed direct numerical simulations of MRI-driven dynamos with explicit
viscosity and resistivity reaching unprecedentedly large values of the
magnetic Prandtl number Pm. In the quasi-stationary state, the magnetic
energy, kinetic energy and the angular momentum transport are approximately
independent of the Reynolds number and follow a universal curve as a function
of Pm. They first increase linearly with Pm up to moderately large values of
this parameter (for $\mathrm{Pm}\lesssim 50$), and smoothly transition to a
regime consistent with a plateau independent of $\mathrm{Pm}$ at
$\mathrm{Pm}\gtrsim 100$. As $\mathrm{Pm}$ is increased, the peak of the
magnetic energy shifts to larger wavenumbers. Interestingly, however, the
energy contained in the largest scales of the magnetic field increases
proportionally to the total magnetic energy, suggesting the presence of a
large-scale dynamo whose efficiency increases with Pm.

These results are particularly important for the formation of magnetars in
fast rotating proto-neutron stars and in neutron star merger remnants. The
increase with Pm of the magnetic energy contained at large scales would
suggest that the dipolar component of the magnetic field should likewise
increase in a spherical model. With this assumption, we may use our results to
extrapolate the results of Reboul-Salze et al. (2021a) obtained at Pm=16 to
the asymptotic regime of very high Pm relevant in a PNS. The fitting formula
(Eq. 5) predicts that the magnetic energy obtained by Reboul-Salze et al.
(2021a) is underestimated due to their moderate value of Pm by a factor
$\mathrm{Pm_{c}}/16\simeq 6$.

In parallel to this work, Held & Mamatsashvili (2022) have also recently
explored the MRI-driven dynamo in the regime of high magnetic Prandtl numbers.
Several aspects of their results are qualitatively consistent with ours: the
stress increases as a power-law with Pm at moderate Pm and transitions to a
weaker dependence, possibly a plateau, at $\mathrm{Pm}>50-100$. Some
differences can nonetheless be noted and they may be explained by the
different setup and parameter space exploration. First, they find a shallower
power-law dependence of the stress with Pm (power law index $0.5-0.7$) than
the linear dependence that we find for moderate $\mathrm{Pm}$. This different
slope is most likely primarily due to their choice of a Keplerian shear rate
($q=1.5$), while we studied a sub-Keplerian shear rate ($q=0.8$)444the three
simulations performed by Held & Mamatsashvili (2022) at $q=0.8$ give a steeper
Pm dependence with a power law index $\simeq 1.4$.. Their slope is consistent
with previously published studies assuming Keplerian shear (see their Figure
17), although one should note that there is also a dependence on the box
aspect ratio, such that Simon & Hawley (2009) for example is consistent with a
linear dependence. The dependence on the shear rate and box aspect ratio
should therefore be studied in more details in the future. Another difference
lies in the parameter space exploration: Held & Mamatsashvili (2022) ran
series of simulations where Pm was varied at fixed Rm (by varying Re) while in
our series of simulations Re was kept fixed and Rm varied. This might play a
role in the other main difference with our results: while in our simulations
the stress and magnetic energy are independent of the Reynolds number (for
fixed Pm), Held & Mamatsashvili (2022) find an additional dependence on the
magnetic Reynolds number at fixed Pm.

The dependence of the MRI saturated state on $\mathrm{Pm}$ (rather than
$\mathrm{Re}$ or $\mathrm{Rm}$) in our simulations highlights the importance
for numerical simulations to describe explicitly the diffusive processes. The
effective $\mathrm{Pm}$ of implicit large eddy simulations (e.g. Mösta et al.,
2015; Kiuchi et al., 2018, for the context of core collapse supernovae and
neutron star mergers) is of order unity regardless of the resolution they may
reach. They are therefore very far from the high-$\mathrm{Pm}$ regime
described in this paper and our results suggest that they may underestimate
the magnetic energy by a factor up to 100.

Beyond the clear increase of the turbulent energies and stress with Pm, our
results provide the first evidence that an asymptotic regime independent of Pm
may exist for $\mathrm{Pm}\gtrsim 100$ (at least for the stress and kinetic
energy). Such a regime had not been obtained by previously published MRI
simulations, because they were restricted to relatively low Pm values
($\mathrm{Pm}\lesssim 16$). We propose that the plateau at very high Pm
probably originates from the scale separation between the velocity and
resistive scales, which prevents an efficient feedback of the resistive scales
on the much larger scales of the flow. One should however caution that our
simulations did not reach a very large scale separation between velocity and
resistive scales, such that a weak dependence of the magnetic energy at high
$\mathrm{Pm}$ cannot be excluded with the present data. An extrapolation of
our results to the much higher values of Pm relevant to proto-neutron stars is
therefore still uncertain and will require a deeper analysis and physical
understanding. In this perspective, it will be important to provide a robust
physical explanation of the $\mathrm{Pm}$ dependence (Riols et al., 2017;
Mamatsashvili et al., 2020) and the asymptotic high-$\mathrm{Pm}$ regime with
a detailed analysis of the energy transfers between different scales.

## Acknowledgements

JG, RR, and MB acknowledge support from the European Research Council
(MagBURST grant 715368). Numerical simulations have been carried out at the
CINES on the Occigen supercomputer (DARI projects A0070410317, A0090410317 and
A0110410317).

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Akiyama et al. (2003) Akiyama S., Wheeler J. C., Meier D. L., Lichtenstadt I., 2003, ApJ, 584, 954
  * Alexakis (2011) Alexakis A., 2011, Phys. Rev. E, 83, 036301
  * Balbus & Hawley (1998) Balbus S. A., Hawley J. F., 1998, Rev. Mod. Phys., 70, 53
  * Balbus & Henri (2008) Balbus S. A., Henri P., 2008, ApJ, 674, 408
  * Bugli et al. (2021) Bugli M., Guilet J., Obergaulinger M., 2021, MNRAS, 507, 443
  * Bugli et al. (2020) Bugli M., Guilet J., Obergaulinger M., Cerdá-Durán P., Aloy M. A., 2020, MNRAS, 492, 58
  * Ferrario & Wickramasinghe (2006) Ferrario L., Wickramasinghe D., 2006, MNRAS, 367, 1323
  * Fromang et al. (2007) Fromang S., Papaloizou J., Lesur G., Heinemann T., 2007, A&A, 476, 1123
  * Goldreich & Lynden-Bell (1965) Goldreich P., Lynden-Bell D., 1965, MNRAS, 130, 125
  * Guilet et al. (2017) Guilet J., Bauswein A., Just O., Janka H.-T., 2017, MNRAS, 471, 1879
  * Guilet & Müller (2015) Guilet J., Müller E., 2015, MNRAS, 450, 2153
  * Guilet et al. (2015) Guilet J., Müller E., Janka H.-T., 2015, MNRAS, 447, 3992
  * Held & Mamatsashvili (2022) Held L. E., Mamatsashvili G., 2022, arXiv e-prints, p. arXiv:2206.00497
  * Kawanaka & Masada (2019) Kawanaka N., Masada Y., 2019, ApJ, 881, 138
  * Kiuchi et al. (2018) Kiuchi K., Kyutoku K., Sekiguchi Y., Shibata M., 2018, Phys. Rev. D, 97, 124039
  * Kiuchi et al. (2014) Kiuchi K., Kyutoku K., Sekiguchi Y., Shibata M., Wada T., 2014, Phys. Rev. D, 90, 041502
  * Kuroda et al. (2020) Kuroda T., Arcones A., Takiwaki T., Kotake K., 2020, ApJ, 896, 102
  * Lander (2021) Lander S. K., 2021, MNRAS, 507, L36
  * Lesur & Longaretti (2005) Lesur G., Longaretti P.-Y., 2005, A&A, 444, 25
  * Lesur & Longaretti (2007) Lesur G., Longaretti P.-Y., 2007, MNRAS, 378, 1471
  * Lesur & Longaretti (2011) Lesur G., Longaretti P.-Y., 2011, A&A, 528, A17
  * Longaretti & Lesur (2010) Longaretti P.-Y., Lesur G., 2010, A&A, 516, A51
  * Mamatsashvili et al. (2020) Mamatsashvili G., Chagelishvili G., Pessah M. E., Stefani F., Bodo G., 2020, ApJ, 904, 47
  * Masada et al. (2007) Masada Y., Sano T., Shibata K., 2007, ApJ, 655, 447
  * Masada et al. (2012) Masada Y., Takiwaki T., Kotake K., Sano T., 2012, ApJ, 759, 110
  * Meheut et al. (2015) Meheut H., Fromang S., Lesur G., Joos M., Longaretti P.-Y., 2015, A&A, 579, A117
  * Metzger et al. (2011) Metzger B. D., Giannios D., Thompson T. A., Bucciantini N., Quataert E., 2011, MNRAS, 413, 2031
  * Mösta et al. (2015) Mösta P., Ott C. D., Radice D., Roberts L. F., Schnetter E., Haas R., 2015, Nature, 528, 376Ð379
  * Mösta et al. (2020) Mösta P., Radice D., Haas R., Schnetter E., Bernuzzi S., 2020, ApJ, 901, L37
  * Obergaulinger et al. (2009) Obergaulinger M., Cerdá-Durán P., Müller E., Aloy M. A., 2009, A&A, 498, 241
  * Pessah et al. (2006) Pessah M. E., Chan C.-K., Psaltis D., 2006, MNRAS, 372, 183
  * Potter & Balbus (2014) Potter W. J., Balbus S. A., 2014, MNRAS, 441, 681
  * Potter & Balbus (2017) Potter W. J., Balbus S. A., 2017, MNRAS, 472, 3021
  * Raynaud et al. (2022) Raynaud R., Cerdá-Durán P., Guilet J., 2022, MNRAS, 509, 3410
  * Raynaud et al. (2020) Raynaud R., Guilet J., Janka H.-T., Gastine T., 2020, Science Advances, 6, eaay2732
  * Reboul-Salze et al. (2021a) Reboul-Salze A., Guilet J., Raynaud R., Bugli M., 2021a, A&A, 645, A109
  * Reboul-Salze et al. (2021b) Reboul-Salze A., Guilet J., Raynaud R., Bugli M., 2021b, arXiv e-prints, p. arXiv:2111.02148
  * Riols et al. (2017) Riols A., Rincon F., Cossu C., Lesur G., Ogilvie G. I., Longaretti P.-Y., 2017, A&A, 598, A87
  * Rossi et al. (2008) Rossi E. M., Armitage P. J., Menou K., 2008, MNRAS, 391, 922
  * Schekochihin et al. (2004) Schekochihin A. A., Cowley S. C., Taylor S. F., Maron J. L., McWilliams J. C., 2004, ApJ, 612, 276
  * Schneider et al. (2019) Schneider F. R. N., Ohlmann S. T., Podsiadlowski P., Röpke F. K., Balbus S. A., Pakmor R., Springel V., 2019, Nature, 574, 211
  * Shi et al. (2016) Shi J.-M., Stone J. M., Huang C. X., 2016, MNRAS, 456, 2273
  * Siegel et al. (2013) Siegel D. M., Ciolfi R., Harte A. I., Rezzolla L., 2013, Phys. Rev. D, 87, 121302
  * Simon & Hawley (2009) Simon J. B., Hawley J. F., 2009, ApJ, 707, 833
  * Takiwaki et al. (2009) Takiwaki T., Kotake K., Sato K., 2009, ApJ, 691, 1360
  * Thompson & Duncan (1993) Thompson C., Duncan R. C., 1993, ApJ, 408, 194
  * Walker & Boldyrev (2017) Walker J., Boldyrev S., 2017, MNRAS, 470, 2653
  * Woosley (2010) Woosley S. E., 2010, ApJ, 719, L204