# Informative Path Planning of Autonomous Vehicle for Parking Occupancy
Estimation

Yunze Hu1†, Jiaao Chen1†, Kangjie Zhou1, Han Gao1, Yutong Li2, and Chang Liu1∗
${\dagger}$Yunze Hu and Jiaao Chen contributed equally to this work.1Yunze Hu,
Jiaao Chen, Kangjie Zhou, Han Gao and Chang Liu are with the Department of
Advanced Manufacturing and Robotics, College of Engineering, Peking
University, Beijing 100871, China<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>(Corresponding author: Chang Liu.).2Yutong Li is with
the Department of Aerospace Engineering, University of Michigan, Ann Arbor, MI
48109, USA<EMAIL_ADDRESS>
###### Abstract

Parking occupancy estimation holds significant potential in facilitating
parking resource management and mitigating traffic congestion. Existing
approaches employ robotic systems to detect the occupancy status of individual
parking spaces and primarily focus on enhancing detection accuracy through
perception pipelines. However, these methods often overlook the crucial aspect
of robot path planning, which can hinder the accurate estimation of the entire
parking area. In light of these limitations, we introduce the problem of
informative path planning for parking occupancy estimation using autonomous
vehicles and formulate it as a Partially Observable Markov Decision Process
(POMDP) task. Then, we develop an occupancy state transition model and
introduce a Bayes filter to estimate occupancy based on noisy sensor
measurements. Subsequently, we propose the Monte Carlo Bayes Filter Tree, a
computationally efficient algorithm that leverages progressive widening to
generate informative paths. We demonstrate that the proposed approach
outperforms the benchmark methods in diverse simulation environments,
effectively striking a balance between optimality and computational
efficiency.

###### Index Terms:

Occupancy Estimation, Path Planning, Active Sensing, POMDP, Monte Carlo Tree
Search.

##  I Introduction

The development of parking occupancy estimation techniques has emerged as a
promising endeavor to enhance parking resource management and mitigate traffic
congestion [1, 2, 3]. A significant research focus has been directed towards
the development of parking occupancy detection systems for pre-installed
surveillance cameras, leveraging advancements in computer vision
technologies[4, 5, 6]. Nonetheless, the reliance on pre-installed cameras
severely constrains the occupancy detection to parking spaces within the
cameras’ field of view, rendering it unable to provide accurate estimates for
entire parking areas, let alone facilities without the presence of such
cameras.

To overcome this limitation, recent studies start employing robotic systems,
such as drones, to obtain accurate occupancy estimation on an ad hoc basis [7,
8]. This approach is particularly suitable for parking areas that lack pre-
installed sensors, such as street parking spaces and temporary parking areas
for events. For instance, Wang and Ren [7] developed an autonomous parking
occupancy detection system that utilized a quadrotor to navigate the parking
area and provide real-time occupancy data. Zhou et al. [8] proposed a vehicle
counting method based on UAV images, which converts vehicle counting into
estimating density values across image pixels.

While these methods offer increased mobility and flexibility for occupancy
estimation, their primary focus lies in improving systems’ perception pipeline
to accurately identify the occupancy states of individual parking spaces [4,
5, 90611552020Quadrotor, 8]. However, the lack of consideration given to the
path planning of robots in these methods can result in limited estimation
accuracy of the entire parking area, since the path taken by the robot greatly
influences the informativeness of sensor observations, which in turn affects
the state estimation accuracy. Furthermore, the utilization of drones requires
dedicated labor and additional expenses for operational maintenance, further
restricting the applicability of existing approaches.

To deal with these problems, this work focuses on informative path planning
for an autonomous vehicle to accurately estimate the occupancy status of a
parking area. The objective is to enable the vehicle to consistently and
accurately estimate the occupancy states of all parking spaces by utilizing
real-time measurements from the onboard sensor, and to dynamically plan an
optimal path that maximizes the informativeness of sensor measurements,
thereby enhancing the overall estimation accuracy across the entire parking
area. Importantly, as the population of intelligent vehicles with advanced
sensing and planning capabilities continues to grow, utilizing such vehicles
for occupancy estimation can present a commercially and technically feasible
alternative to relying solely on drones.

State estimation is a vital component for accurately estimating parking
occupancy, as it involves fusing real-time sensor measurements while
considering the arrival and departure processes of vehicles at each parking
space to determine the current occupancy status. Filtering techniques, such as
the Kalman filter, particle filter, and the more general Bayes filter
framework [9], are effective tools for state estimation in dynamic systems. We
propose modeling the state transition processes of parking spaces and
estimating the parking occupancy status using the Bayes filter framework. This
approach allows us to effectively capture the dynamics of parking space
occupancy and provide reliable estimates.

Planning informative paths to enhance state estimation accuracy presents a
general challenge [10, 11, 12, 13]. One popular solution is formulating the
path planning as a Partially Observable Markov Decision Process (POMDP).
However, obtaining exact solutions to POMDPs is computationally intractable
[14]. Therefore, many approximation algorithms have been proposed to obtain
near-optimal solutions. In recent years, sampling-based online algorithms
utilizing Monte Carlo tree search (MCTS) [15] to solve POMDPs have achieved
noteworthy success, particularly for large-scale problems. For example, Silver
and Veness proposed POMCP that incorporates observation simulations and
updates to MCTS during tree construction [16]. To tackle POMDPs with
continuous or large action and observation spaces, Sunberg and Kochenderfer
proposed POMCPOW and PFT-DPW, utilizing double progressive widening to
restrict the number of child nodes, ensuring that simulations traverse the
same node multiple times to prevent the search tree from being shallow [17].
Based on PFT-DPW, Fischer and Tas developed the Information Particle Filter
Tree that utilizes the information-theoretic reward function to generate
informative paths for information gathering [18] . Despite the progresses
made, these methods rely on specific approximation of belief states, such as
particles, and cannot directly apply to informative path planning for parking
occupancy estimation due to the substantial computational expenses associated
with maintaining a large particle set to achieve satisfactory performance.

To address these limitations, this study proposes a systematic approach to
informative path planning based on MCTS for an autonomous vehicle, aiming to
accurately estimate the occupancy states of all parking spaces in a given
parking area. The main contributions of the work are threefold:

  1. 1.

We define the problem of informative path planning for parking occupancy
estimation and formulate it as a POMDP task. Additionally, We develop an
occupancy state transition model and present a Bayes filter for occupancy
estimation based on noisy sensor measurements.

  2. 2.

We propose the Monte Carlo Bayes Filter Tree (MCBFT) algorithm for informative
path planning. Specifically, we employ progressive widening in the observation
space to avoid constructing a shallow policy tree. Furthermore, we design a
heuristic policy to guide the rollout procedure of MCTS, enabling the
computationally efficient generation of informative paths.

  3. 3.

Through extensive simulations conducted on various parking lot scenarios, we
demonstrate that the solutions obtained by the MCBFT algorithm are consistent
with the optimal solutions in over $90\%$ of scenarios. Moreover, the
computational time of MCBFT is only around $30\%$ of the optimal algorithm,
showcasing a favorable trade-off between optimality and computational
efficiency.

##  II Problem Formulation

Consider a two-dimensional parking lot with $N$ parking spaces. The autonomous
vehicle is assumed to have access to the parking lot map, and is equipped with
a sensor to perceive the occupancy status of the parking spaces within the
field of view (FOV). The autonomous vehicle is tasked with planning the
optimal paths to observe and estimate the occupancy states of the whole
parking area continually and accurately.

###  II-A Environment Model

Figure 1: (a) Environment model and observation model. The coordinate system
is established with the origin $O$ and two coordinate axes
$\overrightarrow{OX}$ and $\overrightarrow{OY}$. The parking lot is arranged
in a $2$-by-$2$ zone configuration, with each zone containing $10$ parking
spaces. The numbers $1$, $5$, $6$, $10$, $11$, $21$ represent the indices of
parking spaces, and $L^{7}$ denotes the coverage set of the $7$th parking
space. The vehicle is at pose $\boldsymbol{p_{k}}$ with the sensor FOV
highlighted in orange, and
$\mathcal{O}(\boldsymbol{p_{k}})=\\{31,32,36,37\\}$. The black solid rectangle
represents that the corresponding parking space is “occupied”, and the dark
gray star symbols depict the discrete locations associated with the pose
points in $V$. (b) Motion model. The vehicle’s current pose is
$\boldsymbol{p_{k}}$ with three available actions: turning left, turning
right, and moving forward. The potential future poses are depicted by vehicle
icons with higher transparency.

A two-dimensional Cartesian coordinate system is established in the parking
lot (Fig. 1(a)) that has a width of $w_{l}$ and a height of $h_{l}$, and is
arranged in a grid of $r$ rows and $c$ columns, with a corridor of length
$L_{c}$ located between each column. The parking lot is divided into $r\cdot
c$ zones, each containing $n_{z}$ parking spaces. Each parking space is of
rectangular shape, with the width denoted as $w_{p}$ and the height denoted as
$h_{p}$. The coordinates of the four vertices of the $i$th parking space are
denoted as $(u^{i}_{1},w^{i}_{1})$, $(u^{i}_{2},w^{i}_{2})$,
$(u^{i}_{3},w^{i}_{3})$, $(u^{i}_{4},w^{i}_{4})$, respectively, where
$i\in\\{1,\cdots,N\\}$. The 2D coverage set of the $i$th parking space is
defined as $L^{i}\triangleq\\{(u,w)\in\mathbb{R}^{2}\ |\
u=\sum_{j=1}^{4}\epsilon_{j}u^{i}_{j},w=\sum_{j=1}^{4}\epsilon_{j}w^{i}_{j},\sum_{j=1}^{4}\epsilon_{j}=1,\epsilon_{j}\in[0,1],j=1,\cdots,4\\}\subset\mathbb{R}^{2}$.
Let $\Delta t>0$ represent the discretization time interval, and the $k$th
time step corresponds to time $t_{k}$. The autonomous vehicle’s pose at $k$th
time step is denoted as $\boldsymbol{p}_{k}=(x_{k},y_{k},\theta_{k})$, where
$x_{k}$ and $y_{k}$ represent the position and $\theta_{k}$ denotes the
orientation. A directed graph $G=(V,E)$ with $m$ vertices can be constructed
such that $V=\\{\boldsymbol{v}^{1},\boldsymbol{v}^{2},\
\cdots,\boldsymbol{v}^{m}\\}$ represents the vehicle’s discrete pose space,
where $\boldsymbol{v}^{j}\in\mathbb{R}^{3}$ corresponds to the $j$th pose
point including the position and orientation with respect to the parking lot
frame, and $E=V\times V$ denotes the directed edge set. If
$\langle\boldsymbol{v}^{i},\boldsymbol{v}^{j}\rangle\in E$, it indicates that
the autonomous vehicle can travel from pose $\boldsymbol{v}^{i}$ to pose
$\boldsymbol{v}^{j}$ within $\Delta t$.

###  II-B Observation Model

Let $\boldsymbol{x}_{k}=(x^{1}_{k},\cdots,x^{N}_{k})\in\\{0,1\\}^{N}$ denote
the ground-truth occupancy states, and
$\boldsymbol{z}_{k}=(z^{1}_{k},\cdots,z^{N}_{k})\in\\{0,1,\emptyset\\}^{N}$
represent the sensor measurements of parking spaces. Here “$1$” means
“occupied”, “$0$” means “unoccupied”, and “$\emptyset$” means the parking
space is not within the sensor’s FOV. We adopt a rectangle-shaped area to
model the sensor’s FOV, noted as
$\mathcal{F}(\boldsymbol{p}_{k})\subset\mathbb{R}^{2}$. Let
$\mathcal{S}(\cdot)$ denote the area covered by a 2D set, and define
$\mathcal{O}(\boldsymbol{p}_{k})\subset\\{1,\dots,N\\}$ as the index set of
parking spaces observable by the vehicle at pose $\boldsymbol{p}_{k}$. The
measurement $z_{k}^{i}$ of occupancy state $x^{i}_{k}$ is obtained when more
than half of the parking space is within the sensor’s FOV, i.e.,

$\mathcal{S}(L^{i}\cap\mathcal{F}(\boldsymbol{p}_{k}))\geq
0.5\mathcal{S}(L^{i})\Leftrightarrow
i\in\mathcal{O}(\boldsymbol{p}_{k})\Leftrightarrow\;z_{k}^{i}\neq\emptyset.$

The sensor follows a probabilistic observation model,

$p(z^{i}_{k}\ |\ x^{i}_{k})=\begin{cases}1,&\text{ $\text{if
}z^{i}_{k}=\emptyset,i\notin\mathcal{O}(\boldsymbol{p}_{k}),$ }\\\
p_{1},&\text{ $\text{if
}z^{i}_{k}=1,x^{i}_{k}=1,i\in\mathcal{O}(\boldsymbol{p}_{k}),$ }\\\
1-p_{1},&\text{ $\text{if
}z^{i}_{k}=0,x^{i}_{k}=1,i\in\mathcal{O}(\boldsymbol{p}_{k}),$ }\\\
p_{2},&\text{ $\text{if
}z^{i}_{k}=0,x^{i}_{k}=0,i\in\mathcal{O}(\boldsymbol{p}_{k}),$ }\\\
1-p_{2},&\text{ $\text{if
}z^{i}_{k}=1,x^{i}_{k}=0,i\in\mathcal{O}(\boldsymbol{p}_{k}),$ }\end{cases}$
(1)

where parameters $p_{1}\in[0,1]$ and $p_{2}\in[0,1]$ represent the probability
of correct observation of occupied and unoccupied parking spaces,
respectively111Note that the actual observation model is sensor dependent and
can be different from the one used in this work. However, the proposed method
can be easily extended to sensors with different observation models..

###  II-C POMDP Model for Parking Occupancy Estimation

The POMDP offers a powerful framework for sequential decision-making under
uncertainty. A finite-horizon POMDP is defined by the 9-tuple
$(D,S,A,\Omega,T,O,R,\gamma,b_{0})$ [19], where the $S$, $A$, $\Omega$ are the
state, action, and observation set, respectively. Function $T$, $O$, $R$
represent the transition model, observation model, and reward function,
respectively. Here $D$ is the planning horizon, $\gamma$ is the discount
factor, and $b_{0}$ is the initial belief. We aim to find a policy
$\boldsymbol{\pi^{*}}$ that can maximize the expected cumulative reward of the
POMDP.

Our work focuses on the sequential decision making for the autonomous vehicle
with the objective of maximizing the acquisition of occupancy information
within a predefined planning horizon, which can be modeled as a finite-horizon
POMDP. The state and observation can be defined as
$\boldsymbol{s}_{k}=(\boldsymbol{p}_{k},\boldsymbol{x}_{k})\in S$,
$\boldsymbol{o}_{k}=\boldsymbol{z}_{k}\in\Omega$, respectively. Each pose
$\boldsymbol{p}_{k}$ has a corresponding discrete action space
$A(\boldsymbol{p}_{k})$, in which each action $\boldsymbol{a}_{k}\in
A(\boldsymbol{p}_{k})$ drives the vehicle from pose $\boldsymbol{p}_{k}$ to
another pose $\boldsymbol{p}_{k+1}\in\\{\boldsymbol{p}\ |\
\langle\boldsymbol{p}_{k},\boldsymbol{p}\rangle\in E\\}$ following the vehicle
motion model $\boldsymbol{p}_{k+1}=f(\boldsymbol{p}_{k},\boldsymbol{a}_{k})$,
as illustrated in Fig. 1(b). The transition model $T$ is defined as

$\boldsymbol{s}_{k+1}=T(\boldsymbol{s}_{k},\boldsymbol{a}_{k}),$ (2)

where $\boldsymbol{p}_{k+1}$ is propagated using the vehicle motion model $f$,
and the transition from $\boldsymbol{x}_{k}$ to $\boldsymbol{x}_{k+1}$ is
based on the parking space state transition model $\sigma$, which will be
specified in Section III-A. The observation
$\boldsymbol{o}_{k}=O(\boldsymbol{s}_{k})$ is obtained based on Eq. 1.

Let $B_{k}=\\{b^{1}_{k},\cdots,b^{N}_{k}\\}$ denote the vehicle’s estimation
of the occupancy state, with the belief state $b^{i}_{k}$ representing the
probability of the $i$th parking space being occupied at time step $k$. An
optimal policy sequence
$\boldsymbol{\pi}^{*}=(\pi^{*}_{k},\cdots,\pi^{*}_{k+D-1})$ is calculated so
that the autonomous vehicle can plan an optimal path to efficiently estimate
the parking occupancy,

$\boldsymbol{\pi^{*}}=\mathop{\arg\max}\limits_{\boldsymbol{\pi}}\mathbb{E}\left[\sum_{t=k}^{k+D-1}\gamma^{t-k}R\
\bigg{|}\ B_{k},\boldsymbol{p}_{k},\boldsymbol{\pi}\right],$ (3)

where $\mathbb{E}$ is the expectation operator, and the reward function $R$ is
related with the estimation of the occupancy state, which will be specified in
Section IV-A. The vehicle then executes the optimal action
$\boldsymbol{a}_{k}^{*}=\pi^{*}_{k}(B_{k},\boldsymbol{p}_{k})\in
A(\boldsymbol{p}_{k})$ and repeats this process to replan at the next time
step in the receding horizon manner, so that the vehicle can proactively
generate informative paths based on the newly received observations.

###  II-D MCTS-Based Solutions to POMDPs

MCTS provides a computationally efficient way to solve POMDPs [15]. MCTS
incrementally constructs an asymmetric policy tree that comprises alternating
layers of state nodes and action nodes through the repetition of four key
phases: selection, expansion, simulation, and back-propagation. During the
selection phase, a strategy that balances exploration and exploitation is
employed to navigate through the current policy tree. Once a leaf node is
reached, new nodes are added to the tree during the expansion phase. Then a
simulation rollout is conducted from the leaf node to estimate the
corresponding reward. Finally, in the back-propagation phase, the parameters
of the visited nodes in this iteration are updated. Once the policy tree has
been constructed, we can choose the optimal solution from the first layer of
action nodes.

##  III Bayes Filtering for Occupancy Estimation

###  III-A Parking Space State Transition Model

We propose a discrete-time model for the transition of individual parking
space states based on existing parking lot traffic models. We assume the
states of all parking spaces are independent of each other.

####  III-A1 Arrival

The arrival of vehicles in a parking lot is usually modeled as a Poisson
process [20]. In a similar vein, we propose to adopt a Poisson distribution to
model the arrival process for each parking space. Let $l$ denote the number of
vehicles that pass the parking space during a time interval and are willing to
park in the parking space if it is unoccupied. The Poisson distribution with
mean parameter $\lambda$ can be expressed as:

$P_{a}(l,\lambda)=\left\\{\begin{array}[]{rcl}\frac{\lambda^{l}}{l!}e^{-\lambda}&&{\text{if}\
l=0,1,2,...}\\\ 0&&{\ \ \ \ \ \ \ \ \text{else}}\end{array}\right..$ (4)

If a parking space is unoccupied at time $t_{k}$, then during time interval
$[t_{k},t_{k+1})$, the probability of no vehicle parking in the space is
$P_{a}(l=0,\lambda)$, and the transition probabilities can be defined as

$\displaystyle p(x^{i}_{k+1}=0\ |\
x^{i}_{k}=0)=P_{a}(l=0,\lambda)=e^{-\lambda},$ (5a) $\displaystyle
p(x^{i}_{k+1}=1\ |\ x^{i}_{k}=0)=1-e^{-\lambda}.$ (5b)

For narrative simplicity, we use $p_{3}$ to represent the transition
probability Eq. 5b, i.e., $p_{3}\triangleq p(x^{i}_{k+1}=1\ |\ x^{i}_{k}=0)$.

####  III-A2 Departure

For each parked vehicle, the parking time is usually modeled as an exponential
distribution [20]. Let $\mu$ denote the departure rate for a parked vehicle,
and the probability density function of the exponential distribution is

$P_{d}(t)=\left\\{\begin{array}[]{rcl}\mu e^{-\mu t}&&{\text{if}\ t>0}\\\
0&&{\ \ \text{else}}\end{array}\right..$ (6)

In the case that a vehicle is parked in the parking space at time $t_{k}$, we
can calculate the probability of the vehicle leaving the space during
$[t_{k},t_{k+1})$ through Eq. 6. The memoryless property of the exponential
distribution leads to the following transition probabilities,

$\displaystyle p(x^{i}_{k+1}=0\ |\ x^{i}_{k}=1)$
$\displaystyle=1-e^{-\mu\Delta t},$ (7a) $\displaystyle p(x^{i}_{k+1}=1\ |\
x^{i}_{k}=1)$ $\displaystyle=e^{-\mu\Delta t}.$ (7b)

For narrative simplicity, we use $p_{4}$ to represent the transition
probability Eq. 7b, i.e., $p_{4}\triangleq p(x^{i}_{k+1}=1\ |\ x^{i}_{k}=1)$.

Eqs. 5 and 7 comprise the parking space state transition model, which we
compactly denote as

$\boldsymbol{x}_{k+1}=\sigma(\boldsymbol{x}_{k}).$ (8)

###  III-B State Estimation with Bayes Filter

We use the Bayes filter to recursively fuse sensor measurements to estimate
the occupancy state of parking spaces. The Bayes filter consists of two steps:
prediction and update.

Prediction. At time step $k$, each parking space’s belief is forward predicted
using the parking space state transition model Eq. 8 in the prediction step,

$\overline{b}^{i}_{k+1}=p_{3}(1-b^{i}_{k})+p_{4}b^{i}_{k},i=1,\cdots,N.$ (9)

Update. The posterior distribution of the occupancy state is calculated by
fusing sensor measurements via the Bayes rule,

${b}^{i}_{k+1}=\eta
p(z_{k+1}^{i}|x_{k+1}^{i})\overline{b}^{i}_{k+1},i=1,\cdots,N,$ (10)

where $p(z_{k+1}^{i}|x_{k+1}^{i})$ is the observation model defined in Eq. 1
and $\eta$ is the normalization factor.

##  IV Monte Carlo Bayes Filter Tree for
Path Planning

###  IV-A Reward Function Design

A belief-based reward function $R$ is essential for quantifying the
information gain associated with executing a specific action. In order to
measure information uncertainty of parking occupancy status, we consider the
entropy of the $i$th parking space, which can be written as

$\mathcal{H}(b_{k}^{i})=-b_{k}^{i}\mathrm{log_{2}}b_{k}^{i}-(1-b_{k}^{i})\mathrm{log_{2}}(1-b_{k}^{i}).$
(11)

Then we select the negation of entropy $-\mathcal{H}$ as the information
measure $\mathcal{I}$ and subsequently define the information gain as
$\Delta\mathcal{I}(b,b^{\prime})=\mathcal{I}(b^{\prime})-\mathcal{I}(b)$,
which has been widely used in the field of active sensing [21]. Following this
idea, the reward function is formulated as

$R=\Delta\mathcal{I}(B_{k},B_{k+1})=\sum_{i=1}^{N}(\mathcal{H}(b_{k}^{i})-\mathcal{H}(b_{k+1}^{i})),$
(12)

where $B_{k+1}$ is calculated through the Bayes filter Eqs. 9 and 10 from the
occupancy estimation at the previous step $B_{k}$, action
$\boldsymbol{a}_{k}=\pi_{k}(B_{k},\boldsymbol{p}_{k})$ and observation
$\boldsymbol{z}_{k+1}$.

###  IV-B Traversal Path Planning

0: $B_{k}:$ the set of beliefs, $\Psi:$ the set of all possible paths

0: $\psi^{*}:$ the path selected

1: for $j=1:N_{T}$ do

2: $\Delta H_{j}=0$

3: for $i=1:N$ do

4: $N_{0}^{i,1}\leftarrow(1,b^{i}_{k})$

5: $H_{j}^{i,0}=\mathcal{H}(b_{k}^{i})$

6: for $d=1:D$ do

7: for $e=1:m_{d-1}$ do

8: $({\phi}_{d-1}^{i,e},b_{k+d-1}^{i,e})\leftarrow N_{d-1}^{i,e}$

9: $\overline{b}_{k+d}^{i,e}=p_{3}(1-b_{k+d-1}^{i,e})+p_{4}b_{k+d-1}^{i,e}$

10: if $i\in\mathcal{O}(\boldsymbol{p}_{k+d}^{j})$ then

11: $b_{k+d}^{i,e}=\eta p_{1}\overline{b}_{k+d}^{i,e}$

112: $\xi=p_{1}\overline{b}_{k+d}^{i,e}+(1-p_{2})(1-\overline{b}_{k+d}^{i,e})$

213: $N_{d}^{i}\leftarrow
N_{d}^{i}\cup\\{({\phi}_{d-1}^{i,e}\xi,b_{k+d}^{i,e}),({\phi}_{d-1}^{i,e}(1-\xi),1-b_{k+d}^{i,e})\\}$

14: else

15: $N_{d}^{i}\leftarrow
N_{d}^{i}\cup({\phi}_{d-1}^{i,e},\overline{b}_{k+d}^{i,e})$

16: end if

17: end for

18: end for

19: $\Delta H_{j}^{i}=0$

20: for $d=1:D$ do

21: $H_{j}^{i,d}=0$

22: for $e=1:m_{d}$ do

23: $H_{j}^{i,d}=H_{j}^{i,d}+{\phi}_{d}^{i,e}\mathcal{H}(b_{k+d}^{i,e})$

24: end for

25: $\Delta H_{j}^{i}=\Delta
H_{j}^{i}+\gamma^{d-1}(H_{j}^{i,d-1}-H_{j}^{i,d})$

26: end for

27: $\Delta H_{j}=\Delta H_{j}+\Delta H_{j}^{i}$

28: end for

29: end for

30: $j^{*}=\mathop{\arg\max}\limits_{j}\Delta H_{j}$

31: $\psi^{*}=\psi_{j^{*}}$

Algorithm 1 Traversal Algorithm

A straightforward strategy for path planning is to enumerate all feasible
paths and compare their associated cumulative rewards. The path with the
highest cumulative reward, as defined in Eq. 3, is selected. This Traversal
path planning algorithm is presented in Alg. 1.

First we generate all $N_{T}$ feasible paths for the upcoming $D$ planning
steps according to the current pose $\boldsymbol{p}_{k}$ and the environment
model. Denote $\Psi=\\{\psi_{1},\cdots,\psi_{N_{T}}\\}$ as the set of all the
feasible paths, and each path
$\psi_{j}=\\{\boldsymbol{p}_{k}^{j},\boldsymbol{p}_{k+1}^{j},\cdots,\boldsymbol{p}_{k+D}^{j}\\}$
consists of the pose sequence of the vehicle within the planning horizon $D$.
For the $i$th parking space, a search tree for simulating possible future
observations is produced. The nodes in the $d$th-layer of the tree are
described as $N_{d}^{i}=\\{N_{d}^{i,1},\cdots,N_{d}^{i,m_{d}}\\}$, where the
$e$th node $N_{d}^{i,e}=({\phi}_{d}^{i,e},b_{k+d}^{i,e})$ consists of the
probability ${\phi}_{d}^{i,e}$ of the realized observation sequence from
planning step $1$ to $d$ and the corresponding belief $b_{k+d}^{i,e}$, and
$m_{d}$ is the total number of distinct observation sequences between step $1$
and $d$. Note that the growth of new nodes in $N_{d}^{i}$ from $N_{d-1}^{i}$
is determined by the observation of the $i$th parking space at planning step
$d$.

For each path, the search tree for the $i$th parking space is built
recursively based on Bayes filter considering all possible future
observations. We first set the initial node $N_{0}^{i,1}=(1,b_{k}^{i})$ (Line
4). At planning step $d$, we extract the information
$(\phi_{d-1}^{i,e},b_{d+k-1}^{i,e})$ from the $e$th node $N_{d-1}^{i,e}$ that
has been calculated in the $(d-1)$th layer (Line 8) , and use them to compute
the belief $\overline{b}_{d+k}^{i,e}$ based on the prediction step Eq. 9 (Line
9). If the vehicle can observe the parking space, the sensor observation model
is utilized in the update step to calculate the belief $b_{d+k}^{i,e}$
according to Eq. 10 (Line 11). Furthermore, we add two new nodes to consider
the two possible observation results at planning step $d$ (Line 12-13), where
$\xi$ denotes the probability of the “occupied” observation. If the vehicle
cannot observe the parking space, only one node is created at planning step
$d$, and the update step of the Bayes filter is skipped (Line 15). Once all
the nodes are expanded to the $D$th layer, we compute the cumulative entropy
reduction $\Delta H^{i}_{j}$ of the $i$th parking space for the whole
prediction horizon (Line 19-26). Then the entropy reduction of all parking
spaces are subsequently summed to compute the overall entropy reduction of the
current path $\Delta H_{j}$ (Line 27). Finally, we select the optimal path
$\psi^{*}=\\{\boldsymbol{p}_{k}^{*},\boldsymbol{p}_{k+1}^{*},\cdots,\boldsymbol{p}_{k+D}^{*}\\}$
with the highest cumulative entropy reduction of the parking lot $\Delta
H_{j^{*}}$ (Line 30-31) and choose the associated pose
$\boldsymbol{p}^{*}=\boldsymbol{p}_{k+1}^{*}$ as the next pose point for the
vehicle to head for.

###  IV-C Monte Carlo Bayes Filter Tree for Path Planning

While the Traversal algorithm can obtain optimal solutions through
enumeration, the time complexity exponentially increases with growing planning
horizon $D$. Inspired by [17, 18], we propose MCBFT, an MCTS-based algorithm
for informative path planning, which enables us to evaluate parking occupancy
under low computational burden with minor accuracy loss compared to the
Traversal algorithm, as shown in Alg. 2. The Search procedure serves as the
entry point for the decision-making process. It takes the initial belief
$B_{k}$ and the vehicle’s initial pose $\boldsymbol{p}_{k}$ as inputs,
constructs a policy tree, and selects the optimal action $\boldsymbol{a}^{*}$
when the simulation times reach the preset number $I$.

0: Search$(B_{k},\boldsymbol{p}_{k})$

1: for $i\in 1:I$ do

2: Simulate($B_{k},\boldsymbol{p}_{k},0$)

3: end for

14: return
$\boldsymbol{a}^{*}=\mathop{\arg\max}\limits_{\boldsymbol{a}^{\prime}}Q(B_{k},\boldsymbol{p}_{k},\boldsymbol{a}^{\prime})$

4: Simulate$(B,\boldsymbol{p},h)$

5: if $h=D$ then

6: return $0$

7: else

8: if $|C(B,\boldsymbol{p})|<|A(\boldsymbol{p})|$ then

9: $\boldsymbol{a}\leftarrow$ sample from
$A(\boldsymbol{p})-C(B,\boldsymbol{p})$

10: $C(B,\boldsymbol{p})\leftarrow
C(B,\boldsymbol{p})\cup\\{\boldsymbol{a}\\}$

11: end if

12:
$\boldsymbol{a}^{*}=\mathop{\arg\max}\limits_{\boldsymbol{a}^{\prime}}Q(B,\boldsymbol{p},\boldsymbol{a}^{\prime})+c\sqrt{\frac{\log
q(B,\boldsymbol{p})}{q(\boldsymbol{a}^{\prime})}}$

13: if $|C(\boldsymbol{a}^{*})|<\kappa q(\boldsymbol{a}^{*})^{\delta}$ then

214: $\boldsymbol{x}\leftarrow$ sample from $B$

315: $(\boldsymbol{p}^{\prime},\boldsymbol{x}^{\prime})\leftarrow
T((\boldsymbol{p},\boldsymbol{x}),\boldsymbol{a}^{*})$

416: $\boldsymbol{o}\leftarrow
O(\boldsymbol{p}^{\prime},\boldsymbol{x}^{\prime})$

517: $B^{\prime}\leftarrow$ Bayes Filter$(B,\boldsymbol{o})$

618: $\omega(B^{\prime},\boldsymbol{p}^{\prime})\leftarrow p(\boldsymbol{o}\
|\ B,\boldsymbol{p}^{\prime})$

19: if $(B^{\prime},\boldsymbol{p}^{\prime})\notin C(\boldsymbol{a}^{*})$ then

20: $C(\boldsymbol{a}^{*})\leftarrow
C(\boldsymbol{a}^{*})\cup\\{(B^{\prime},\boldsymbol{p}^{\prime})\\}$

21: $J\leftarrow\Delta
I(B,B^{\prime})+\gamma$Rollout($B^{\prime},\boldsymbol{p}^{\prime},h+1$)

22: else

23: $J\leftarrow\Delta
I(B,B^{\prime})+\gamma$Simulate($B^{\prime},\boldsymbol{p}^{\prime},h+1$)

24: end if

25: else

26: $(B^{\prime},\boldsymbol{p}^{\prime})\leftarrow$ sample
$(B^{\prime},\boldsymbol{p}^{\prime})$ w.p.
$\frac{\omega(B^{\prime},\boldsymbol{p}^{\prime})}{\sum_{\overline{B^{\prime}}}{\omega(\overline{B^{\prime}},\boldsymbol{p}^{\prime})}}$

27: $J\leftarrow\Delta
I(B,B^{\prime})+\gamma$Simulate($B^{\prime},\boldsymbol{p}^{\prime},h+1$)

28: end if

29: $q(B,\boldsymbol{p})\leftarrow q(B,\boldsymbol{p})+1$

30: $q(\boldsymbol{a}^{*})\leftarrow q(\boldsymbol{a}^{*})+1$

31: $Q(B,\boldsymbol{p},\boldsymbol{a}^{*})\leftarrow
Q(B,\boldsymbol{p},\boldsymbol{a}^{*})+\frac{J-Q(B,\boldsymbol{p},\boldsymbol{a}^{*})}{q(\boldsymbol{a}^{*})}$

32: end if

733: return $J$

33: Rollout$(B,\boldsymbol{p},h)$

34: if $h=D$ then

35: return $0$

36: else

837: $\boldsymbol{a}^{*}\leftarrow$ RolloutPolicy($B,\boldsymbol{p}$)

938: $\boldsymbol{x}\leftarrow$ sample from $B$

1039: $(\boldsymbol{p}^{\prime},\boldsymbol{x}^{\prime})\leftarrow
T((\boldsymbol{p},\boldsymbol{x}),\boldsymbol{a}^{*})$

1140: $\boldsymbol{o}\leftarrow
O(\boldsymbol{p}^{\prime},\boldsymbol{x}^{\prime})$

1241: $B^{\prime}\leftarrow$ Bayes Filter$(B,\boldsymbol{o})$

42: $J\leftarrow\Delta
I(B,B^{\prime})+\gamma$Rollout($B^{\prime},\boldsymbol{p}^{\prime},h+1$)

43: end if

44: return $J$

Algorithm 2 Monte Carlo Bayes Filter Tree

The policy tree comprises two distinct types of nodes: belief nodes and action
nodes. The former contains a tuple of the occupancy estimation and the
vehicle’s pose, denoted as $(B,\boldsymbol{p})$, while the latter consists of
the vehicle’s action $\boldsymbol{a}$. Compared to PFT-DPW [17] and IPFT [18]
that utilize the particle filter to approximate the belief, MCBFT stores
precise belief in each belief node $(B,\boldsymbol{p})$, which allows for the
exact calculation of the total belief-based reward $J$. We maintain a record
$\langle q,C,\boldsymbol{o},\omega\rangle$ for each belief node
$(B,\boldsymbol{p})$ in the tree, which consists of its number of visits $q$,
set of child nodes $C$, corresponding observation $\boldsymbol{o}$, and
probability $\omega$ of generating $\boldsymbol{o}$ from the observation
space. Likewise, for each action node $\boldsymbol{a}$, we track $\langle
q,C,Q\rangle$, which includes its number of visits $q$, set of child nodes
$C$, and its estimated reward $Q$ based on belief $B$ and pose
$\boldsymbol{p}$, where $B$ and $\boldsymbol{p}$ are obtained from the parent
node of action node $\boldsymbol{a}$.

The Simulate function progressively constructs the search tree of depth
$h\in\\{0,1,\cdots,D\\}$. An action node $\boldsymbol{a}^{*}$ is selected
according to the upper confidence bound (UCB) (Line 12) [15]. As the number of
parking spaces that the autonomous vehicle can observe increases, the
observation space increases exponentially, and the probability of sampling the
same observation decreases exponentially, leading to constructing a shallow
tree. To prevent this issue, we use progressive widening to limit the number
of belief nodes linked to the same parent node $\boldsymbol{a}^{*}$ to $\kappa
q(\boldsymbol{a}^{*})^{\delta}$, where $\kappa$ and $\delta$ are preset hyper-
parameters (Line 13), so that the simulations can pass through the same child
node multiple times. If the action node $\boldsymbol{a}^{*}$ does not have
enough child nodes, a new observation $\boldsymbol{o}$ is generated using the
observation model $O$, and the posterior belief $B^{\prime}$ is computed with
the Bayes filter (Line 16-17). The $\omega$ is calculated based on the
observation model Eq. 1 (Line 18). On the other hand, if the action node
$\boldsymbol{a}^{*}$ already has enough child nodes, an existing belief child
node is sampled proportional to $\omega$ (Line 26). If the node
$(B^{\prime},\boldsymbol{p}^{\prime})$ has already been created, the Simulate
function is called recursively to visit and create deeper nodes (Line 23, 27).
Otherwise, the Rollout function is performed to estimate $J$ (Line 21, 34-44).
We introduce a heuristic method in Rollout to efficiently generate informative
paths. Specifically in RolloutPolicy function (Line 37), instead of random
action selection, the optimal action is chosen probabilistically based on the
Traversal algorithm with a preset planning horizon $D^{\prime}$. Finally, the
parameters of visited nodes are updated (Line 29-31).

##  V Simulation Results and Analyses

The proposed path planning approach is evaluated on various parking lot
scenarios. First, three manually designed, representative scenarios are
created to serve as a sanity check, qualitatively validating the proposed
approach in generating paths that align with intuitive expectations. Then,
$200$ random scenarios are generated to quantitatively evaluate the generality
of the proposed algorithm.

###  V-A Simulation Setup

We utilize the Matlab Automated Parking Valet toolbox to construct parking
lots. We construct three distinct parking lot models, the parameters of which
are listed in Tab. I. The locations of pose points are distributed in the
corridors, inside the zones, and around the periphery of the zones, and then
the pose points are interconnected based on the motion model of the vehicle.
Each parking space measures $w_{p}=3m$ and $h_{p}=6m$, and the parking space
state transition model Eq. 8 takes parameters $\mu=0.000378$ and
$\lambda=0.000624$, which is derived from the actual parking lot data [20].
The dimensions of FOV are $w_{f}=20m$ and $h_{f}=10m$, and the observation
model takes parameters $p_{1}=p_{2}=0.95$. The parameters in MCBFT are set as
$I=100$ and $D^{\prime}=5$. The simulations are conducted on a laptop with
Intel Core i7-1065G7<EMAIL_ADDRESS>and 16GB RAM.

###  V-B Qualitative Evaluation Using Three Representative Scenarios

We design three representative scenarios based on Parking Lot Model i@ to
evaluate the ability of the proposed algorithm in planning informative paths.
For each scenario, we conduct $10$ simulations, where the parking space state
transitions and the sensor measurements are simulated stochastically. We
denote the $30$ parking spaces in the bottom left, bottom right, mid left, mid
right, top left and top right area of the parking lot as the zone i@ to zone
vi@, respectively. If the $i$th parking space is not observed a priori, then
$b_{0}^{i}$ is set to be $0.5$. Otherwise, if the $i$th parking space is
initially observed as occupied or unoccupied, then $b_{0}^{i}$ is initialized
as $0.95$ or $0.05$, respectively. We evaluate the performance of the MCBFT
and the Traversal algorithm, both with a planning horizon of $20$ steps,
namely MCBFT-$20$ and Traversal-$20$, respectively. The baseline algorithms
include the Greedy algorithm, which essentially corresponds to the Traversal
algorithm with planning horizon $D=1$, and the Random Walk algorithm, where
the vehicle randomly selects an action at each pose. Fig. 2 visualizes the
three scenarios and the generated paths under different path planning
approaches.

TABLE I: Parking Lot Model Parameters

Model | $w_{l}(m)$ | $h_{l}(m)$ | $L_{c}(m)$ | $r$ | $c$ | $n_{z}$
---|---|---|---|---|---|---
Model i@ | $144$ | $56$ | $18$ | $3$ | $2$ | $30$
Model ii@ | $108$ | $130$ | $18$ | $7$ | $2$ | $18$
Model iii@ | $153$ | $74.5$ | $18$ | $4$ | $3$ | $18$

Figure 2: Path comparison in three representative scenarios. The red, blue,
cyan, and purple lines represent paths planned by MCBFT-$20$, Traversal-$20$,
the Greedy algorithm and the Random Walk algorithm, respectively. The colored
rectangle outlining each parking space represents the initial entropy, which
ranges from $0$ to $1$ as illustrated in the colorbar in (a). (a) Scenario<EMAIL_ADDRESS>(b) Scenario<EMAIL_ADDRESS>(c) Scenario<EMAIL_ADDRESS>
####  V-B1 Scenario i@

The initial pose of the autonomous vehicle is set as $(19.5,9.0,0.0)$, and the
right half of the parking lot is observed a priori (Fig. 2(a)). The paths
generated by MCBFT-$20$, Traversal-$20$, and the Greedy algorithm can all
accomplish sequential observations of the unobserved zones, thus rapidly
decreasing the information uncertainty. In contrast, in $90\%$ of the
simulations, the Random Walk algorithm produces unfavorable paths that
repeatedly detect zones that have already been observed, leading to small
reduction in entropy.

####  V-B2 Scenario ii@

The autonomous vehicle’s pose is now initialized at $(124.5,28.0,0.0)$, and
only the four parking spaces in zone i@ are previously observed (Fig. 2(b)).
To maximize the decrease of uncertainty, it is more reasonable for the
autonomous vehicle to make a left turn and observe zone vi@ and zone v@
sequentially, which corresponds to the blue Traversal planning path shown in
Fig. 2(b). In $8$ out of the $10$ simulation outcomes, MCBFT-$20$ generates
the same path as Traversal-$20$. In contrast, when utilizing the Greedy
algorithm or the Random Walk algorithm, only in $4$ out of $10$ simulations
does the vehicle first observe zone vi@ and<EMAIL_ADDRESS>
####  V-B3 Scenario iii@

The autonomous vehicle starts at pose $(61.5,28.0,0.0)$, and is most uncertain
about the parking occupancy status of zone iii@ (Fig. 2(c)). Under this
circumstance, both Traversal-$20$ and MCBFT-$20$ choose to turn around and
observe zone iii@ first. In contrast, in the Greedy and the Random Walk
algorithm, the number of simulations where the vehicle performs such behavior
is only $0$ and $1$ out of $10$, respectively, and in the remaining
simulations, the vehicle tends to stay in the already observed zones,
resulting in little information gain.

###  V-C Quantitative Evaluation via Randomly Generated Scenarios

To further evaluate the generality of the proposed approach, we consider two
different parking lot models, Parking Lot Model ii@ and iii@, and randomly
generate $100$ scenarios for each model. For each scenario, we randomly select
a set of zones to be pre-observed. Observed zones are initially assigned
random belief values between $0.3$ and $0.95$ for “occupied” spaces, and
between $0.05$ and $0.7$ for “unoccupied” spaces. Beliefs of unobserved zones
are set to be $0.5$. The pose of the autonomous vehicle is also randomly
initialized.

We evaluate the performance of five algorithms: the MCBFT with a horizon of
$10$ steps, namely MCBFT-$10$, the Greedy algorithm, the Random Walk
algorithm, and the Traversal algorithm with planning horizons of 10 and 5
steps, namely Traversal-$10$ and Traversal-$5$, respectively. All tested
algorithms share the same ground-truth occupancy status time series
$(\boldsymbol{x}_{0},\cdots,\boldsymbol{x}_{K})$, where $K$ denotes the
simulation length and is set to be three-fourths of the total number of
discrete locations reachable by the vehicle. We record the entropy $H_{k}$ of
the belief and the correctness rate $\alpha_{k}$ of the estimation, defined as
the ratio of the parking spaces where the estimation is consistent with the
actual state. Specifically, the vehicle estimates spaces with a belief
exceeding $0.6$ as “occupied”, and spaces with a belief below $0.4$ as
“unoccupied”. When the belief lies between $0.4$ and $0.6$, the vehicle cannot
confidently estimate the state of parking space, resulting in an “unsure”
estimation that is inconsistent with any actual state.

After simulation, we evaluate the performance of the aforementioned five
algorithms based on the three metrics: the improvement of the correctness rate
$\Delta\alpha\triangleq\alpha_{K}-\alpha_{0}$, the average computational time
$\overline{\tau}$ ($\text{s}/\text{step}$), and the entropy reduction
percentage $\frac{\Delta H}{H_{0}}$, where $\Delta H$ is defined as
$H_{0}-H_{K}$ that represents the decrease in entropy. Specifically, we denote
$N_{\alpha},N_{\Delta H}$ as the number of scenarios in which MCBFT-$10$
outperforms another algorithm in terms of $\Delta\alpha$ and $\frac{\Delta
H}{H_{0}}$, respectively. To further compare the path decisions of MCBFT and
the Traversal algorithm, we define the pose points where multiple actions are
available as the decision points, and introduce the concept of consistency
rate $\beta$, which is defined as the ratio of decision points where the
MCBFT-$10$ and Traversal-$10$ take the same action.

Figure 3: (a) Illustration of Parking Lot Model<EMAIL_ADDRESS>(b) and (c) compare the
simulation results of five algorithms in terms of entropy reduction percentage
$\frac{\Delta H}{H_{0}}$ and improvement of correctness rate $\Delta\alpha$,
respectively. TABLE II: Simulation Results in Parking Lot Model ii@

Alg | MCBFT-10 | Traversal-10 | Traversal-5 | Greedy | Random Walk
---|---|---|---|---|---
$\overline{\tau}$ | $0.35$ | $0.97$ | $4.8\times 10^{-2}$ | $2.7\times 10^{-3}$ | $1.6\times 10^{-5}$
$N_{\alpha}$ | $/$ | $45$ | $73$ | $84$ | $100$
$N_{\Delta H}$ | $/$ | $45$ | $79$ | $81$ | $100$

####  V-C1 Parking Lot Model ii@

We consider a parking lot model (Fig. 3(a)) that has four more rows than
Parking Lot Model<EMAIL_ADDRESS>The number of locations accessible to the vehicle is
$125$, and the simulation runs for $94$ time steps. As Fig. 3 shows,
MCBFT-$10$ and Traversal-$10$ outperform the other algorithms on both metrics:
$\Delta\alpha$ and $\frac{\Delta H}{H_{0}}$. As shown in Tab. II, MCBFT-$10$
demonstrates superior performance compared to Traversal-$5$ and the Greedy
algorithm in most of the scenarios. In nearly half of the scenarios, the
performance of MCBFT-$10$ surpasses that of Traversal-$10$. Moreover,
MCBFT-$10$ demonstrates an average running time that is approximately $36\%$
of that required by Traversal-$10$. Out of the $2634$ decision points,
MCBFT-$10$ and Traversal-$10$ make the same decisions in $2465$ instances,
resulting in $\beta=93.58\%$. These results demonstrate that MCBFT can obtain
optimal actions in the majority of cases while achieving high computational
efficiency.

####  V-C2 Parking Lot Model iii@

We design a parking lot model (Fig. 4(a)) that has one more column and row
than the Parking Lot Model<EMAIL_ADDRESS>The number of locations the vehicle can travel
to is $104$, and the simulation runs for $78$ time steps. The results, as
depicted in Fig. 4 and Tab. III, show similar trends as Section V-C1. The
$\beta$ reaches $92.84\%$, as MCBFT-$10$ and Traversal-$10$ select the same
option in $2413$ out of $2599$ decision points.

Figure 4: (a) Illustration of Parking Lot Model<EMAIL_ADDRESS>(b) and (c) compare the
simulation results of five algorithms in terms of entropy reduction percentage
$\frac{\Delta H}{H_{0}}$ and improvement of correctness rate $\Delta\alpha$,
respectively. TABLE III: Simulation Results in Parking Lot Model iii@

Alg | MCBFT-10 | Traversal-10 | Traversal-5 | Greedy | Random Walk
---|---|---|---|---|---
$\overline{\tau}$ | $0.20$ | $0.67$ | $4.0\times 10^{-2}$ | $2.3\times 10^{-3}$ | $1.2\times 10^{-5}$
$N_{\alpha}$ | $/$ | $55$ | $75$ | $84$ | $100$
$N_{\Delta H}$ | $/$ | $53$ | $77$ | $82$ | $100$

##  VI Conclusion

We present the problem of informative path planning for autonomous vehicles to
estimate parking occupancy and formulate it under the POMDP framework. We
develop MCBFT that utilizes progressive widening to mitigate the high
computational cost and enable online path planning. Simulation results show
that MCBFT achieves a favorable trade-off between optimality and computational
efficiency. Future work includes evaluating the proposed approach in realistic
parking lot environments and extension to multi-vehicle parking occupancy
estimation.

## References

  * [1] S. Yang and Z. S. Qian, “Turning meter transactions data into occupancy and payment behavioral information for on-street parking,” Transportation Research Part C: Emerging Technologies, vol. 78, pp. 165–182, 2017.
  * [2] A. Varghese and G. Sreelekha, “An efficient algorithm for detection of vacant spaces in delimited and non-delimited parking lots,” IEEE Transactions on Intelligent Transportation Systems, vol. 21, no. 10, pp. 4052–4062, 2019\.
  * [3] B. Assemi, A. Paz, and D. Baker, “On-street parking occupancy inference based on payment transactions,” IEEE Transactions on Intelligent Transportation Systems, vol. 23, no. 8, pp. 10680–10691, 2021.
  * [4] R. M. Nieto, A. Garcia-Martin, A. G. Hauptmann, and J. M. Martinez, “Automatic vacant parking places management system using multicamera vehicle detection,” IEEE Transactions on Intelligent Transportation Systems, vol. 20, no. 3, pp. 1069–1080, 2018.
  * [5] J. K. Suhr and H. G. Jung, “Automatic parking space detection and tracking for underground and indoor environments,” IEEE Transactions on Industrial Electronics, vol. 63, no. 9, pp. 5687–5698, 2016.
  * [6] R. Ke, Y. Zhuang, Z. Pu, and Y. Wang, “A smart, efficient, and reliable parking surveillance system with edge artificial intelligence on iot devices,” IEEE Transactions on Intelligent Transportation Systems, vol. 22, no. 8, pp. 4962–4974, 2021.
  * [7] Y. Wang and B. Ren, “Quadrotor-enabled autonomous parking occupancy detection,” in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2020.
  * [8] H. Zhou, L. Wei, M. Fielding, D. Creighton, S. Deshpande, and S. Nahavandi, “Car park occupancy analysis using uav images,” in 2017 IEEE International Conference on Systems, Man, and Cybernetics (SMC), pp. 3261–3265, 2017.
  * [9] S. Thrun, W. Burgard, and D. Fox, Probabilistic Robotics.  MIT Press, 2005.
  * [10] C. Liu and J. K. Hedrick, “Model predictive control-based target search and tracking using autonomous mobile robot with limited sensing domain,” 2017 American Control Conference (ACC), pp. 2937–2942, 2017.
  * [11] Y. Kantaros, B. Schlotfeldt, N. Atanasov, and G. J. Pappas, “Sampling-based planning for non-myopic multi-robot information gathering,” Autonomous Robots, vol. 45, no. 7, pp. 1029–1046, 2021.
  * [12] B. Du, K. Qian, C. Claudel, and D. Sun, “Parallelized active information gathering using multisensor network for environment monitoring,” IEEE Transactions on Control Systems Technology, vol. 30, no. 2, pp. 625–638, 2021\.
  * [13] S. C. Pinto, S. B. Andersson, J. M. Hendrickx, and C. G. Cassandras, “Multiagent persistent monitoring of targets with uncertain states,” IEEE Transactions on Automatic Control, vol. 67, no. 8, pp. 3997–4012, 2022\.
  * [14] C. H. Papadimitriou and J. N. Tsitsiklis, “The complexity of markov decision processes,” Mathematics of Operations Research, vol. 12, no. 3, pp. 441–450, 1987.
  * [15] C. Browne, E. J. Powley, D. Whitehouse, S. M. M. Lucas, P. I. Cowling, P. Rohlfshagen, S. Tavener, D. P. Liebana, S. Samothrakis, and S. Colton, “A survey of monte carlo tree search methods,” IEEE Transactions on Computational Intelligence and AI in Games, vol. 4, pp. 1–43, 2012.
  * [16] D. Silver and J. Veness, “Monte-carlo planning in large pomdps,” Advances in Neural Information Processing Systems, vol. 23, 2010.
  * [17] Z. Sunberg and M. J. Kochenderfer, “Online algorithms for pomdps with continuous state, action, and observation spaces,” in International Conference on Automated Planning and Scheduling, 2017.
  * [18] J. Fischer and Ö. S. Tas, “Information particle filter tree: An online algorithm for pomdps with belief-based rewards on continuous domains,” in International Conference on Machine Learning, 2020.
  * [19] M. Lauri, D. Hsu, and J. Pajarinen, “Partially observable markov decision processes in robotics: A survey,” IEEE Transactions on Robotics, vol. 39, pp. 21–40, 2022.
  * [20] J. Xiao, Y. Lou, and J. Frisby, “How likely am i to find parking? – a practical model-based framework for predicting parking availability,” Transportation Research Part B: Methodological, 2018.
  * [21] L. S. Mihaylova, T. Lefebvre, H. Bruyninckx, K. Gadeyne, and J. D. Schutter, “A comparison of decision making criteria and optimization methods for active robotic sensing,” in Numerical Methods and Application, 2002.


# The gap-size influence on the excitation of magnetorotational instability in
cylindrical Couette flows

G. Rüdiger1,2<EMAIL_ADDRESS>M. Schultz2  1University of Potsdam, Institute
of Physics and Astronomy, Karl-Liebknecht-Str. 24-25, 14476 Potsdam, Germany
2Leibniz-Institut für Astrophysik Potsdam, An der Sternwarte 16, D-14482
Potsdam, Germany

###### Abstract

The excitation conditions of the magnetorotational instability are studied for
axially unbounded Taylor-Couette flows of various gap widths between the
cylinders. The cylinders are considered as made from both perfect-conducting
or insulating material and the conducting fluid with a finite but small
magnetic Prandtl number rotates with a quasi-Keplerian velocity profile. The
solutions are optimized with respect to the wave number and the Reynolds
number of the rotation of the inner cylinder. For the axisymmetric modes we
find the critical Lundquist number of the applied axial magnetic field the
lower the wider the gap between the cylinders. A similar result is obtained
for the induced cell structure: the wider the gap the more spherical the cells
are. The marginal rotation rate of the inner cylinder – for fixed size of the
outer cylinder – always possesses a minimum for not too wide and not too
narrow gap widths. For perfect-conducting walls the minimum lies at $r_{\rm
in}\simeq 0.4$ while it is at $r_{\rm in}\simeq 0.5$ for insulating walls
where $r_{\rm in}$ is the normalized radius of the inner cylinder. The lowest
magnetic field amplitudes to excite the instability are required for Taylor-
Couette flows between perfect-conducting cylinders with gaps corresponding to
$r_{\rm in}\simeq 0.2$. For even wider and also for very thin gaps the needed
magnetic fields and rotation frequencies are shown to become rather huge.

Also the nonaxisymmetric modes with $|m|=1$ have been considered. Their
excitation generally requires stronger magnetic fields and higher magnetic
Reynolds numbers in comparison to those for the axisymmetric modes which is
primarily true for wide-gap containers with $r_{\rm in}\lower
3.01385pt\hbox{$\;\stackrel{{\scriptstyle\textstyle<}}{{\sim}}\;$}0.3$.

###### keywords:

Taylor-Couette flow – magnetohydrodynamics – magnetorotational instability

##  1 Introduction and motivation

Couette flows with conducting fluids between two rotating cylinders are
favorable for an experimental realization of the various versions of the
magnetorotational instability (MRI) for which the applied magnetic field is
always current-free (Velikhov, 1959; Rüdiger & Zhang, 2001; Ji et al., 2001;
Schartman et al., 2009, 2012). If the field is spiral rather than axial the
necessary Reynolds numbers and Hartmann numbers are surprisingly small
(Hollerbach & Rüdiger, 2005; Rüdiger et al., 2006) which made it possible to
investigate the corresponding helical version of MRI (HMRI) in much detail in
the Promise experiment (Stefani et al., 2006). The resulting instability modes
are axisymmetric and they are migrating along the rotation axis. The closely
related azimuthal MRI (AMRI) appears when working with current-free toroidal
fields for which the unstable modes are nonaxisymmetric (Ogilvie & Pringle,
1996; Hollerbach et al., 2010; Seilmayer et al., 2014).

The experiments mentioned above have used one and the same container
construction where the inner radius was 50% of the outer radius, i.e. $r_{\rm
in}=R_{\rm in}/R_{\rm out}=0.5$. For this geometry the rotation shear
$\mu={\it\Omega_{\rm out}}/{\it\Omega_{\rm in}}$ for quasi-Keplerian rotation
is $\mu=r_{\rm in}^{1.5}=0.35$. In a recent paper the Princeton group
presented experimental results related to the standard version of MRI (SMRI),
with a purely axial field being applied, using a container of $r_{\rm
in}=0.35$ and the aspect ratio $\Gamma=H/(R_{\rm out}-R_{\rm in})=2.1$ (Wang
et al., 2022b).

After previous results the absolute minimum of the magnetic Reynolds number
(optimized with respect to axial wave number) necessary for the excitation of
marginal instability is 21 for perfectly conducting walls and 14 for
insulating walls for $r_{\rm in}=0.5$ and for a shear flow with $\mu=0.33$
(Rüdiger et al., 2018a). For their experiment Wang et al. (2022b) report a
rotation rate ratio of $\mu=0.19$ slightly smaller than the value 0.21 for
quasi-Keplerian rotation.

The model with a homogeneous fluid contained between two vertically-infinite
rotating cylinders is used with an applied uniform magnetic field parallel to
the rotation axis. For viscous flows the most general form of the rotation law
${\it\Omega}(R)$ of the fluid is

${{\it\Omega}}(R)=a+{b\over R^{2}},$ (1)

where $a$ and $b$ are two constants related to the angular velocities
${\it\Omega}_{\rm{in}}$ and ${\it\Omega}_{\rm{out}}$ with which the inner and
the outer cylinders are rotating and $R$ is the distance from the rotation
axis. If $R_{\rm{in}}$ and $R_{\rm{out}}$ ($R_{\rm{out}}>R_{\rm{in}}$) are the
radii of the two cylinders then

$a={\mu-{r_{\rm in}}^{2}\over 1-{r_{\rm in}}^{2}}{{\it\Omega}}_{\rm in}\ \ \ \
\ {\textrm{and}}\ \ \ \ b=R_{\rm{in}}^{2}{1-\mu\over 1-{r_{\rm
in}}^{2}}{{\it\Omega}}_{\rm in}$ (2)

with the geometry parameters

$\mu={{{\it\Omega}}_{\rm{out}}\over{\it\Omega}_{\rm{in}}}\qquad{\textrm{and}}\qquad
r_{\rm in}={R_{\rm{in}}\over R_{\rm{out}}}.$ (3)

Following the Rayleigh stability criterion,

${{\rm d}(R^{2}{{\it\Omega}})^{2}\over{\rm d}R}>0,$ (4)

rotation laws are hydrodynamically stable for $a>0$, i.e. $\mu>r_{\rm
in}^{2}$. They should in particular be stable for resting inner cylinder, i.e.
$\mu\to\infty$.

The present paper has two motivations. The first one concerns the question of
the dependence of the eigenvalues on the gap width of the container. The
aspect ratio $r_{\rm in}$ of the cylinder radii is the only free parameter
describing the geometry of the axially unbounded container. It is so far
unknown how the gap width determines the critical rotation rate for given
magnetic field and also the wave number of the excited instability pattern.
The latter result will have consequences for necessary vertical extension of a
possible experimental setup.

We work with the magnetic Prandtl number

${\rm Pm}={\nu\over\eta},$ (5)

with $\nu$ as the kinematic viscosity and $\eta$ as the magnetic diffusivity.
The equations are solved here for the small magnetic Prandtl number Pm
$=10^{-5}$ close to the value for liquid sodium. As known, for small magnetic
Prandtl numbers the excitation conditions for the standard magnetorotational
instability only depend on the microscopic magnetic diffusivity rather than
the molecular viscosity, hence they do not depend on the actual value of
$\rm{Pm}$. The ratio of the container wall radii is varied from $r_{\rm
in}=0.1$ to $r_{\rm in}=0.9$. We shall see that between these values the
critical magnetic Reynolds number of rotation possesses a minimum while the
critical Lundquist number of the applied magnetic field linearly grows with
the aspect ratio $r_{\rm in}$.

The second question to be attacked is the difference of the excitation
conditions for axisymmetric and nonaxisymmetric modes. Though it is widely
known that the mode with the easiest excitation is the axisymmetric one it is
still important to know how much more difficult the excitation of a
nonaxisymmetric mode is.

##  2 Basic equations

The MHD equations which have to be solved are

${\partial{\mbox{\boldmath$u$}}\over\partial
t}+({\mbox{\boldmath$u$}}\nabla){\mbox{\boldmath$u$}}=-{1\over\rho}\nabla
p+\nu\Delta{\mbox{\boldmath$u$}}+{1\over\rho}{\mbox{\boldmath$J$}}\times{\mbox{\boldmath$B$}}$
(6)

and

${\partial{\mbox{\boldmath$B$}}\over\partial t}={\rm
curl}({\mbox{\boldmath$u$}}\times{\mbox{\boldmath$B$}})+\eta\Delta{\mbox{\boldmath$B$}},$
(7)

with the electric current

${\mbox{\boldmath$J$}}={1\over\mu_{0}}{\rm curl}{\mbox{\boldmath$B$}}$ (8)

and ${\rm div}\ {\mbox{\boldmath$u$}}={\rm div}{\mbox{\boldmath$B$}}=0$. They
are considered in cylindrical geometry with $R$, $\phi$, and $z$ as the
coordinates. A viscous electric-conducting incompressible fluid between two
rotating infinite cylinders in the presence of a uniform magnetic field
parallel to the rotation axis leads to the basic solution
$U_{R}=U_{z}=B_{R}=B_{\phi}=0,B_{z}=B_{0}={\rm const}\ {\rm and}\
U_{\phi}=aR+b/R$, with $U$ as the flow and $B$ as the magnetic field. We are
interested in the stability of this solution. The perturbed state of the flow
may be described by
$u^{\prime}_{R},\;u^{\prime}_{\phi},\;u^{\prime}_{z},\;p^{\prime},B^{\prime}_{R},\;B^{\prime}_{\phi},\;B^{\prime}_{z}$
with $p^{\prime}$ as the pressure perturbation.

In the following only the linear stability problem will be considered. By
expansion of the disturbances into normal modes the solutions of the
linearized equations are of the form

$\displaystyle{\mbox{\boldmath$u$}}^{\prime}={\mbox{\boldmath$u$}}(R)\ {\rm
e}^{{\rm i}(m\phi+kz+\omega t)},\ \ \ \ \ \ \ \ \ \ \ \
{\mbox{\boldmath$B$}}^{\prime}={\mbox{\boldmath$B$}}(R)\ {\rm e}^{{\rm
i}(m\phi+kz+\omega t)}.$ (9)

From here on all dashes are being omitted from the symbols of fluctuating
quantities. The marginal stability line is defined where the imaginary part
$\Im({\omega})$ vanishes. We shall always use the geometrical average

$R_{0}=\sqrt{(R_{\rm out}-R_{\rm in})R_{\rm in}}$ (10)

as the unit of length, $\eta/R_{0}$ as the unit of velocity and $B_{0}$ as the
unit of the magnetic field. We note the rather weak dependence of $R_{0}$ on
the value of $r_{\rm in}$ as long as $0.2\leq r_{\rm in}\leq 0.8$. The $R_{0}$
only becomes small for $r_{\rm in}\to 0$ or for $r_{\rm in}\to 1$, i.e. for
very wide or for very narrow gaps between the cylinders. Its maximum is
reached for $r_{\rm in}=0.5$. In order to denormalize quantities, $R_{0}^{-1}$
is used as the unit of wave numbers and $\nu/R_{0}^{2}$ as the unit of
frequencies.

Using the same symbols for normalized quantities, the equations can be written
as a system of 10 equations of first order, i.e.

$\displaystyle\frac{{\rm d}u_{R}}{{\rm d}R}=-\frac{u_{R}}{R}-{\rm
i}\frac{m}{R}u_{\phi}-{\rm{i}}ku_{z},$ (11) $\displaystyle\frac{{\rm
d}u_{\phi}}{{\rm d}R}=X_{2}-\frac{u_{\phi}}{R},$ (12) $\displaystyle\frac{{\rm
d}u_{z}}{{\rm d}R}=X_{3},$ (13) $\displaystyle\frac{{\rm d}X_{1}}{{\rm
d}R}=\left(\frac{m^{2}}{R^{2}}+k^{2}\right)u_{R}+{\rm i}(\omega+m{{\rm
Re}}\,\Omega)u_{R}+$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad+2{\rm{i}}\frac{m}{R^{2}}u_{\phi}-2{\rm{Re}}\,\Omega
u_{\phi}-{\rm{i}}k{\rm{Ha}}^{2}B_{R},$ (14) $\displaystyle\frac{{\rm
d}X_{2}}{{\rm
d}R}=\left(\frac{m^{2}}{R^{2}}+k^{2}\right)u_{\phi}+{\rm{i}}(\omega+m{\rm{Re}}\,\Omega)u_{\phi}-$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad-2{\rm{i}}\frac{m}{R^{2}}u_{R}+2a{\rm{Re}}u_{R}-{\rm{i}}k{\rm
Ha}^{2}B_{\phi}+$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad+\frac{m^{2}}{R^{2}}u_{\phi}+k\frac{m}{R}u_{z}-{\rm{i}}\frac{m}{R}X_{1},$
(15) $\displaystyle\frac{{\rm d}X_{3}}{{\rm
d}R}=\left(\frac{m^{2}}{R^{2}}+k^{2}\right)u_{z}+{\rm{i}}(\omega+m{\rm
Re}\,\Omega)u_{z}-$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad-\frac{X_{3}}{R}-{\rm
i}k{\rm{Ha}}^{2}B_{z}+k\frac{m}{R}u_{\phi}+k^{2}u_{z}-{\rm{i}}kX_{1},$ (16)
$\displaystyle\frac{{\rm d}B_{R}}{{\rm d}R}=-\frac{B_{R}}{R}-{\rm
i}\frac{m}{R}B_{\phi}-{\rm i}kB_{z},$ (17) $\displaystyle\frac{{\rm
d}B_{\phi}}{{\rm d}R}=X_{4}-\frac{B_{\phi}}{R},$ (18) $\displaystyle\frac{{\rm
d}B_{z}}{{\rm d}R}={\rm{i}}\left(\frac{m^{2}}{kR^{2}}+k\right)B_{R}-\frac{{\rm
Pm}}{k}(\omega+m{\rm Re}\,\Omega)B_{R}+u_{R}-\frac{m}{kR}X_{4},$ (19)
$\displaystyle\frac{{\rm d}X_{4}}{{\rm
d}R}=\left(\frac{m^{2}}{R^{2}}+k^{2}\right)B_{\phi}+{\rm i}{\rm
Pm}(\omega+m{\rm Re}\,\Omega)B_{\phi}-$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad-2{\rm{i}}\frac{m}{R^{2}}B_{R}-{\rm
i}ku_{\phi}+2{\rm PmRe}\frac{b}{R^{2}}B_{R}.,$ (20)

where $X_{1}$ is given by

$X_{1}=\frac{{\rm d}u_{R}}{{\rm d}R}+\frac{u_{R}}{R}-P-{\rm Ha}^{2}B_{z}$ (21)

with $P$ as the pressure fluctuation.

Here the dimensionless Reynolds number $\rm{Re}$ and the Hartmann number
$\rm{Ha}$ are defined as

${\rm Re}={R_{0}^{2}{{\it\Omega}}_{\rm in}\over\nu}\ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ {\rm Ha}={R_{0}B_{0}\over\sqrt{\mu_{0}\rho\nu\eta}}.$ (22)

For given Hartmann number and magnetic Prandtl number we shall compute with a
linear theory the critical Reynolds number of the rotation of the inner
cylinder, also for various azimuthal mode numbers $m$. We shall see that the
excitation conditions for SMRI can easily be expressed by the magnetic
Reynolds number $\rm{Rm}$ and the Lundquist number $\rm{S}$, with the
definitions

${\rm{Rm}}={R_{0}^{2}{{\it\Omega}}_{\rm in}\over\eta}\ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ {\rm{S}}={R_{0}B_{0}\over\sqrt{\mu_{0}\rho}\eta}$ (23)

without influence of the molecular viscosity. The ratio of both quantities
forms the magnetic Mach number of rotation,

$\rm{Mm}={\rm{Rm}\over\rm{S}},$ (24)

which describes the strength of the rotation normalized with the applied
magnetic field. The majority of cosmic objects is characterized by magnetic
Mach numbers larger than unity (except the magnetars). We shall discuss the
run of the magnetic Mach number on the normalized gap width between the
cylinders only for the characteristic constellation where the Reynolds number
is minimal for the excitation of the instability.

If for given size of the container and for given magnetic diffusivity the
rotation of the inner cylinder is needed for which the instability is excited
then the alternative dimensionless numbers

${\rm{Rm}}_{\rm out}={R_{\rm out}^{2}{{\it\Omega}}_{\rm in}\over\eta}\ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ {\rm{S}}_{\rm out}={R_{\rm
out}B_{0}\over\sqrt{\mu_{0}\rho}\eta}$ (25)

can be used for the normalized rotation rate and magnetic field strength,
respectively.

The below calculations were specifically carried out for the small magnetic
Prandtl number $\rm{Pm}=10^{-5}$ but as we know that the obtained eigenvalues
$\rm{Rm}$ and $\rm{S}$ are also correct for even smaller $\rm{Pm}$ (Rüdiger &
Shalybkov, 2002). The reason is that for small $\rm{Pm}$ the critical Reynolds
number runs with $1/\rm{Pm}$ so that the magnetic Reynolds number
$\rm{Rm}\simeq$ const. It is clear, therefore, that for $\rm{Pm}=0$ the
magnetorotational instability does not exist.

##  3 Boundary conditions

For the solution of the differential equations of 10th order a set of 10
boundary conditions is needed. Always no-slip conditions for the velocity on
the walls are used, i.e. $u_{R}=u_{\phi}={\rm d}u_{R}/{\rm d}R=0.$ The
magnetic boundary conditions depend on the electrical properties of the walls.
For perfectly conducting walls the tangential currents and the radial
component of the magnetic field vanish hence ${\rm d}B_{\phi}/{\rm
d}R+B_{\phi}/R=B_{R}=0$. These boundary conditions hold for both $R=R_{\rm
in}$ and $R=R_{\rm out}$.

For insulating walls the magnetic boundary conditions are different at $R_{\rm
in}$ and $R_{\rm out}$, i.e. for $R_{\rm in}$

$B_{R}+{\rm
i}\frac{B_{z}}{I_{m}(kR)}\left(\frac{m}{kR}I_{m}(kR)+I_{m+1}(kR)\right)=0,$
(26)

and for $R=R_{\rm out}$

$B_{R}+{\rm
i}\frac{B_{z}}{K_{m}(kR)}\left(\frac{m}{kR}K_{m}(kR)-K_{m+1}(kR)\right)=0,$
(27)

where $I_{m}$ and $K_{m}$ are the modified Bessel functions. The condition for
the toroidal field is $kRB_{\phi}=mB_{z}$ (Rüdiger et al., 2003). Neutral
stability of the solutions is reached for vanishing $\Im(\omega)$.

The homogeneous set of equations (2.6) – (20) with the boundary conditions
included determine the eigenvalue problem of the form ${\cal L}(k,m,{\rm
Re},{\rm Ha},\omega)=0$ for given Pm. $\cal L$ is a complex quantity, both its
real part and its imaginary part must vanish for the critical Reynolds number.
For nonaxisymmetric modes the real part, ${\Re}(\omega)$, of $\omega$
describes a drift of the pattern along the azimuth. It is the second quantity
fixed by the complex eigenequation. For a fixed Hartmann number, a fixed
Prandtl number and a given vertical wave number we also find the critical
$\rm{Re}$ of the system. It is minimal for a certain wave number defining a
marginally unstable mode. The corresponding value $\rm{Re}_{\rm min}$ is the
minimal Reynolds number if $\rm{Ha}$ is varied and $\rm{Ha}_{\rm min}$ is the
Hartmann number for which the $\rm{Re}_{\rm min}$ occurs. For nonaxisymmetric
modes the real part of the frequency $\omega$ is the second eigenvalue fixed
by the eigenequation.

##  4 General results

For perfect-conducting boundary conditions the resulting curves of marginal
stability are given in Fig. 1 and for vacuum boundary conditions in Fig. 2. In
both plots the left graph presents the magnetic Reynolds numbers (optimized by
the wave number) and the right one gives the resulting wave numbers, in both
cases as functions of the Lundquist number. The dashed lines belong to the
nonaxisymmetric modes with $|m|=1$. Generally, the latter require higher
values $\rm{Rm}_{\rm min}$ than the axisymmetric modes. The curves for $m=0$
exhibit their typical shape: They are rather steep for $\rm{S}<\rm{S}_{\rm
min}$ and they are much flatter for $\rm{S}>\rm{S}_{\rm min}$ where
$\rm{S}_{\rm min}$ is the minimum of the function $\rm{Rm}=\rm{Rm}(\rm{S})$.
For $\rm{S}=$O(1) the influence of the microscopic diffusion stops the
existence of the SMRI while for much larger $\rm{S}$ this existence is limited
by too strong external magnetic fields.

The curves for the axisymmetric modes only possess a lower limit of the
Reynolds number while the curves for the nonaxisymmetric modes always possess
a lower and an upper limit of the (normalized) rotation frequencies. The
nonaxisymmetric modes are stabilized by too fast differential rotation
(Rädler, 1986; Rüdiger et al., 2018a).

###  4.1 Medium gaps

The minimal $\rm{Rm}_{\rm min}$ in Fig. 1 for medium gap widths approximately
behave according to $(1-r_{\rm in})\rm{Rm}_{\rm min}\simeq 13$ (see Table 1).
This relation implies that the minimal rotation rates of the inner cylinder
for the considered gaps behave like ${\it\Omega}_{\rm in}\propto 1/(1-r_{\rm
in})$, i.e. in containers with wider gaps the instability is easier to excite.
For very wide gaps, however, both the magnetic Reynolds number as well as the
needed rotation rate of the inner cylinder grow to very large values. Minimal
rotation rates are only possible for experiments with medium $r_{\rm in}$.

In Table 1 the characteristic eigenvalues $\rm{Rm}_{\rm min}$ and $\rm{S}_{\rm
min}$ have been collected for the minima of the curves in Figs. 1 and 2 which
we shall call the characteristic values. Obviously, the influences of the
boundary conditions are only small for the containers with large gaps. One
finds the $\rm{Rm}_{\rm min}$ with vacuum boundary condition as always smaller
than for perfect-conductor conditions (except for very thin gaps) but the
$\rm{S}_{\rm min}$ are always larger (see also Fig. 4). Consequently, the
corresponding magnetic Mach numbers are much larger for cylinders made from
perfect conductors. For such containers the axial wave numbers are always
larger than for insulating walls.

With our normalizations the vertical extent $\delta z$ of one cell of the
instability pattern, normalized by the gap width $D=R_{\rm out}-R_{\rm in}$
between the cylinders, is given by

$\zeta={\delta z\over D}={\pi\over k}\sqrt{{r_{\rm in}\over 1-r_{\rm in}}}.$
(28)

Flat cells are described by $\zeta<1$ while axially elongated cells possess
$\zeta>1$. For $r_{\rm in}=0.5$ it is simply $\zeta={\pi}/{k}$ so that for
$k\simeq\pi$ the cell is almost circular. The examples given in Table 1,
however, show that the wave numbers do not reach the value of $\pi$ for medium
$r_{\rm in}$ hence the cells are always elongated in the axial direction. This
is true for models with both sorts of boundary conditions. For perfect-
conducting cylinder walls the $\zeta>1$ hardly varies with the $r_{\rm in}$ –
the cells are always oblong. Evidently, such cells are not suitable to provide
the angular momentum transport in accretion disks or galaxies.

For perfect-conducting cylinders it follows $\zeta\simeq 1.8$ for almost all
$r_{\rm in}$ while it can become significantly larger for insulating material.
The minimum axial extend of a container probing MRI pattern is thus $H\simeq
1.8D$. The aspect ratio $H/D$ of the Princeton experiment is 2.1. It happens
that the numerical $\zeta$-values differ by almost a factor of two for models
with the same geometry but with different boundary conditions

Figure 1: Left: Stability maps for $m=0$ (solid lines): $r_{\rm in}=0.1$,
$r_{\rm in}=0.2$, $r_{\rm in}=0.3$ (green), $r_{\rm in}=0.4$, $r_{\rm in}=0.5$
(red), $r_{\rm in}=0.6$, $r_{\rm in}=0.7$ $r_{\rm in}=0.8$. For $m=1$ (dashed
lines): $r_{\rm in}=0.3$ (green), $r_{\rm in}=0.4$, $r_{\rm in}=0.5$ (red),
$r_{\rm in}=0.6$. Right: The corresponding axial wave numbers. Quasi-Keplerian
differential rotation, $\rm{Pm}=10^{-5}$, perfect-conducting cylinder
material.

Figure 2: The unstable modes for containers with insulating cylinders. Left:
stability maps for $m=0$ (solid lines): $r_{\rm in}=0.1$, $r_{\rm in}=0.3$
(green), $r_{\rm in}=0.5$ (red), $r_{\rm in}=0.7$, $r_{\rm in}=0.9$ (blue).
For $m=1$ (dashed lines): $r_{\rm in}=0.3$ (green), $r_{\rm in}=0.4$, $r_{\rm
in}=0.5$ (red), $r_{\rm in}=0.7$. Right: the corresponding axial wave numbers.
Quasi-Keplerian differential rotation, $\rm{Pm}=10^{-5}$. Table 1: The
coordinates of the minima of the profiles in Figs. 1 and 2 (both left) for
several radii the inner cylinder for different boundary conditions (left:
perfect conduction; right: vacuum). $m=0$, $\rm{Pm}=10^{-5}$. All models are
for quasi-Keplerian rotation law.

$r_{\rm in}$ $\mu$ $\rm{Rm}_{\rm min}$ $\rm{S}_{\rm min}$ $k$ $\zeta$
$\rm{Mm}$ 0.1 0.031 20.0 1.58 0.7 1.4 12.7 0.2 0.089 17.8 2.05 1.0 1.6 8.81
0.3 0.16 18.5 2.76 1.2 1.7 6.93 0.4 0.25 20.8 3.38 1.4 1.8 6.15 0.5 0.35 24.7
4.11 1.7 1.8 6.00 0.6 0.46 31.1 5.00 2.1 1.8 6.22 0.7 0.59 42.0 6.32 2.6 1.8
6.64 0.8 0.72 64.1 8.22 3.4 1.8 7.79 0.9 0.85 131 12.6 5.1 1.8 10.4 0.95 0.93
265 17.6 7.3 1.9 15.0 $r_{\rm in}$ $\mu$ $\rm{Rm}_{\rm min}$ $\rm{S}_{\rm
min}$ $k$ $\zeta$ $\rm{Mm}$ 0.1 0.031 30.3 2.78 1.4 0.7 10.9 0.2 0.089 20.8
3.04 1.1 1.4 6.84 0.3 0.16 16.9 3.62 1.1 1.9 4.67 0.4 0.25 16.0 4.82 1.1 2.3
3.32 0.5 0.35 16.5 5.21 1.2 2.6 3.17 0.6 0.46 18.3 6.47 1.2 3.2 2.83 0.7 0.59
22.3 8.54 1.2 4.0 2.61 0.8 0.72 30.7 12.0 1.2 5.2 2.55 0.9 0.85 56.6 22.1 1.0
9.4 2.56 0.95 0.93 109 39.8 0.92 14.9 2.74

###  4.2 Extremal gaps, thin-shell approximation

We note that for very wide gaps with $r_{\rm in}\simeq 0.1$ the characteristic
Lundquist numbers $\rm{S}_{\rm min}$ – for which the associated Reynolds
number is minimal – are reduced to values of order unity while for very narrow
gaps with $r_{\rm in}\simeq 0.95$ the Lundquist numbers are maximal. In both
limits the Reynolds numbers $\rm{Rm}_{\rm min}$ possess enlarged values. As
also the linear dimension $R_{0}$ becomes small for small and/or large $r_{\rm
in}$ the necessary inner rotation rates become very large excluding the
applicability of containers with very wide and/or very narrow gaps between the
cylinders for experiments. One finds that containers with $0.3\lower
3.01385pt\hbox{$\;\stackrel{{\scriptstyle\textstyle<}}{{\sim}}\;$}r_{\rm
in}\lower
3.01385pt\hbox{$\;\stackrel{{\scriptstyle\textstyle<}}{{\sim}}\;$}0.6$ require
the least rotation rates for excitation of standard magnetorotational
instability. On the other hand, for both sorts of boundary conditions the
models with $r_{\rm in}\simeq 0.2$ require the weakest magnetic fields.

Table 1 also gives the the results for very thin gaps between the cylinders up
to $r_{\rm in}=0.95$ (Donnelly & Ozima, 1960, 1962). The most striking
difference due to the choice of the boundary conditions is here the numerical
value of the calculated magnetic Mach number (last column). For $r_{\rm in}\to
1$ the characteristic Reynolds number grows to larger and larger values. We
did not find a maximum of $\rm{Rm}_{\rm min}$ for $r_{\rm in}\to 1$. As the
$R_{0}^{2}$ runs with $d=1-r_{\rm in}$ for $d\to 0$ we find ${\it\Omega_{\rm
in}}\to\infty$ in this limit. It should thus not be possible to work with a
thin-shell approximation (Edmonds, 1958) for numerical or experimental
realizations of the standard MRI in TC flows.

##  5 The nonaxisymmetric modes

The excitation of nonaxisymmetric modes requires faster rotation and stronger
magnetic fields than the excitation of the axisymmetric modes (Fig. 1, Fig. 2
, dashed lines). For medium $r_{\rm in}$ of about 0.5 (red) the lines of
marginal stability do hardly depend on the radius of the inner cylinder. As
also the geometric radius $R_{0}$ is almost constant for different $r_{\rm
in}$ the rotation frequencies and magnetic field amplitudes needed to excite
nonaxisymmetric modes are almost the same for such values of $r_{\rm in}$.
However, for weak magnetic fields with $\rm{S}<\rm{S}_{\rm min}$ the Reynolds
numbers for excitation of the $m=1$ modes are much higher than those for the
excitation of the axisymmetric modes with $m=0$. We note that the curves for
the weak-field branch with $\rm{S}<\rm{S}_{\rm min}$ become very steep so that
the excitation of nonaxisymmetric modes requires very rapid rotation.

For wide gaps between the cylinders ($r_{\rm in}=0.3$, green) one finds that
the lowest Reynolds number belongs to much higher Lundquist numbers than for
$r_{\rm in}\approx 0.5$. In addition, the curves for weak Lundquist numbers
are much steeper than the curve for the corresponding axisymmetric mode.
Hence, for wide gaps with $r_{\rm in}\simeq 0.3$ it is almost impossible to
excite the axisymmetric and the nonaxisymmetric mode simultaneously by
experiments with Lundquist numbers not much higher than unity (see Wang et al.
(2022a)).

Figure 3: Drift rates $\omega_{\rm dr}=\Re(\omega)/{\it\Omega_{\rm in}}$ of
the nonaxisymmetric modes $m=1$ (dashed lines): $r_{\rm in}=0.3$ (green),
$r_{\rm in}=0.4$, $r_{\rm in}=0.5$ (red), $r_{\rm in}=0.6$. Left: perfect-
conducting cylinder material, Right: insulating cylinder material. Quasi-
Keplerian differential rotation, $\rm{Pm}=10^{-5}$.

The nonaxisymmetric modes are drifting in azimuthal direction. The drift
rates, $\omega_{\rm dr}=\Re(\omega)/{\it\Omega_{\rm in}}$, are given in Fig.
3. According to the relation

${{\partial\phi/\partial t}\over{\it\Omega_{\rm in}}}=-{\omega_{\rm dr}\over
m},$ (29)

the negative $\omega_{\rm dr}$ plotted in the figures indicate a migration of
the patterns in the direction of the global rotation - for both sorts of
boundary conditions. In all cases the azimuthal migration is faster then the
rotation of the outer cylinder ($\omega_{\rm dr}>\mu$).

##  6 Marginal stability for weak fields with $\rm{S}<\rm{S}_{\rm min}$

It is possible to operate with weak magnetic fields, i.e. with
$\rm{S}<\rm{S}_{\rm min}$. Then, however, the Reynolds numbers necessary for
excitation basically grow and the unstable wave numbers become smaller, i.e.
the cells become longer. Figures 4 demonstrate these weak-field solutions for
$r_{\rm in}=0.3$ and $r_{\rm in}=0.4$ where the vertical lines belong to the
magnetic fields ($\rm{S}=0.96$ for $B=2150$ G and $\rm{S}=1.22$ for $B=2750$
G) reported for gallium SMRI experiments by Wang et al. (2022b). We have also
to note the influence of the boundary conditions: for insulating cylinders the
curves are so steep that the actual Reynolds numbers exceed the $\rm{Rm}_{\rm
min}$ by more than an order of magnitude. The effect is reduced for perfect-
conducting cylinder material but still the enhancement of the critical
Reynolds number is by a factor of three relative to $\rm{Rm}_{\rm min}$. This
is even a minimum value as the real cylinders are by far not perfect-
conducting. For galinstan as the fluid and stainless steel as the (outer)
cylinder material the ratio

$\hat{\sigma}={\sigma_{\rm cyl}\over\sigma_{\rm fl}}$ (30)

of the electric conductivities of the cylinders and the fluid is 0.47 which
does neither well approach perfect-conduction nor vacuum boundary conditions.
For sodium experiments one finds $\hat{\sigma}\simeq 0.16$, hence the vacuum
boundary conditions might provide the appropriate description. Rüdiger et al.
(2018b) derived the form of the boundary conditions with finite values of (30)
and have shown that for $\hat{\sigma}$ of order unity the resulting
eigenvalues can approximately be interpolated between the values for
$\hat{\sigma}=0$ (insulating boundaries) and $\hat{\sigma}=\infty$ (perfect-
conducting boundaries).

Figure 4: Stability map of the axisymmetric modes for $r_{\rm in}=0.3$ (left)
and $r_{\rm in}=0.4$ (right) in relation to experimental results of (Wang et
al., 2022b). Containers with perfect-conducting cylinders (red) and with
insulating cylinders (black). The vertical lines mark the Lundquist numbers
after (23) for the Princeton MRI-experiments represented by the lowest and the
highest solid circles in their Fig. 2a while the asterisks belong the used
magnetic Reynolds numbers. The vertical lines correspond to applied magnetic
fields of 2150 G and 2750 G, resp. With the planned Dresdyn sodium experiment
the same magnetic field amplitudes belong to Lundquist numbers exceeding $15$.

Detailed consequences of the application of weak magnetic fields for the
excitation of the axisymmetric mode are shown in Fig. 4. The Lundquist number
(in our definition, see (23)) is marked by vertical lines for applied magnetic
fields of 2150 G and 2750 G for a container with $R_{\rm out}=20.3$ cm,
$r_{\rm in}=0.35$ filled with galinstan ($\rho=6.4$ g cm-3 and $\eta=2428$ cm2
s-1). The magnetic Prandtl number of galinstan is $1.4\times 10^{-6}$. These
numbers correspond to the gallium experiment by Wang et al. (2022b). The
maximally possible uniform magnetic field in this experiment given as 4800 G
corresponds to a Lundquist number after Eq. (23) of $\rm{S}=2.1$. One finds
with Fig. 4 that for such fields the minimum magnetic Reynolds number for
marginal stability must exceed 20, corresponding to a minimum rotation rate
with $Rm=10.8$ in the notation of Wang et al. (2022b). The left vertical lines
in Fig. 4 represent the alleged SMRI realization for $B_{0}=0.2$ and $Rm=3.4$
(lower asterisk) provided by their Fig. 2a. The experiment with the fastest
rotation ($Rm=4.5$, upper asterisk) belongs to a Lundquist number of
$\rm{S}=1.22$ (right vertical lines).

If the typical parameters of the experiment with almost Keplerian flow and
with maximal $Rm=4.5$ (in their notation) are transformed to our definitions
(23), one obtains maximal Reynolds numbers of $\rm{Rm}=8.4$ (right asterisks
in the plots) which does not reach, however, the values required for marginal
instability at any Lundquist number. After Fig. 4 this deficit is particularly
drastic for vacuum boundary conditions. It is thus risky for experiments to
work with Lundquist numbers which are basically smaller than the $\rm{S}_{\rm
min}$ given in the Table 1.

Our results comply with the finding of Goodman & Ji (2002) in their Figs. 1
and 2, that containers with insulating walls of $r_{\rm in}=0.33$ including a
quasi-Keplerian flow do not allow the excitation of SMRI with an applied field
of less than 2750 G. Even with perfect-conducting cylinders the curves are so
steep for $\rm{S}<\rm{S}_{\rm min}$ that the needed rotation rates of the
cylinders are unrealistically high.

##  7 Flat cells

As already demonstrated in Table 1 the cell structure for the axisymmetric
modes is nearly spherical for wide gaps and rather elongated in axial
direction for narrow gaps. It is thus basically unclear whether the standard
MRI is also able to produce flat cells with $\zeta\ll 1$ which are necessary
to appear in experiments with a flat container ($H<D$) and/or in flat cosmical
objects such as accretion disks and galaxies. For the latter, however, the
magnetic Mach number of rotation does not exceed values of 10 as it is the
case for the solutions given in Table 1.

We have to probe, therefore, whether eigensolutions exist for finite $\rm{Rm}$
and $\rm{S}$ when, for example, $\zeta=0.1$ is required. From (28) one obtains
for (say) $r_{\rm in}=0.1$ that solutions with $k=10.5$ are matched for
$\zeta=0.1$ and $k=21$ for $\zeta=0.05$.

Figure 5: Lines of marginal stability for $r_{\rm in}=0.1$ and a fixed axial
wave number of $k=10.5$ for $\zeta=0.1$ and $k=21$ for $\zeta=0.05$. Perfect-
conducting cylinders (red), insulating cylinders (green). The magnetic Mach
numbers of the solutions grow for growing flatness. The influence of the
boundary conditions is rather weak. $m=0$.

Figure 5 shows the results. Indeed, the envisaged flat cells exist in the
axially unbounded container and even for similar Lundquist numbers as for the
solutions with the lowest Reynolds numbers and the elongated cells. The actual
Reynolds numbers for flat cell structures, however, are much higher than
before. Hence, the magnetic Mach numbers for flat cells are also higher than
for the elongated cells summarized in Table 1. They are larger than 10 and
grow for growing flatness. This result complies with that of a global model of
a flat galaxy for magnetic Prandtl number $\rm{Pm}\geq 1$ (Kitchatinov &
Rüdiger, 2004). The more flat the cells, the stronger the dissipation and the
harder the differential rotation must work to excite the instability while the
magnetic field needed for the rotation minimum remains unchanged. We also
learn from Fig. 5 that the influence of the actual boundary condition is
remarkably weak.

##  8 Discussion

The influence of the position $r_{\rm in}$ of the inner cylinder of Taylor-
Couette flows on the excitation of the magnetorotational instability has been
studied. To demonstrate the results we shall switch to the representations
(25) of the magnetic Reynolds number and the Lundquist number which for given
outer cylinder size $R_{\rm out}$ form minimal normalized inner rotation rates
and magnetic field amplitudes needed for excitation of the instability. These
quantities are plotted as function of $r_{\rm in}$ by Fig. 6 for both sorts of
boundary conditions. The main result is that too narrow or too wide gaps would
require very high rotation rates or very strong magnetic fields. For
$0.3\lower
3.01385pt\hbox{$\;\stackrel{{\scriptstyle\textstyle<}}{{\sim}}\;$}r_{\rm
in}\lower
3.01385pt\hbox{$\;\stackrel{{\scriptstyle\textstyle<}}{{\sim}}\;$}0.6$ the
dependence of the eigenvalues on $r_{\rm in}$ is rather weak. The minima for
conducting cylinders are at $r_{\rm in}\simeq 0.4$ for the rotation rate and
at $r_{\rm in}\simeq 0.2$ for the magnetic field. We also note that the
influence of the boundary conditions is opposite for rotation and field. For
vacuum conditions the needed rotation rates are mostly lower than for perfect-
conduction conditions but the needed magnetic fields are higher for insulating
cylinders.

Figure 6: The magnetic Reynolds number (left) and the Lundquist number (right)
after the definitions (25) representing the normalized inner rotation rate and
the magnetic field amplitude needed for excitation vs. the inner cylinder
position. Quasi-Keplerian differential rotation, $\rm{Pm}=10^{-5}$, perfect-
conducting cylinder material (red), insulating cylinders (black). The numbers
are taken from Table 1.

The experiment with the fastest rotation ($Rm=4.5$) by Wang et al. (2022b)
corresponds to $\rm{Rm}_{\rm out}=37$ which is certainly below the curves in
Fig. 6 (left). A similar situation holds with respect to the magnetic field:
the given maximal possible field of 4800 G provides a Lundquist number of
$\rm{S}_{\rm out}=4.4$ which again does not reach the minimum value in the
right plot of Fig. 6. Compared with our calculations the reported experiments
are subcritical. There is no value of $r_{\rm in}$ for which the considered
rotation rates and/or magnetic fields are supercritical by a large margin. For
a correct interpretation of the obtained results the modification of the
angular velocity profiles and the specific role of the Ekman-Hartmann layers
close to the copper lids must be understood in much more details.

###### Acknowledgements.

Frank Stefani (Dresden-Rossendorf) is acknowledged for discussions of the
presented problem and a critical reading of the manuscript.

## References

  * Donnelly & Ozima (1960) Donnelly, R. J. & Ozima, M. 1960 Hydromagnetic Stability of Flow Between Rotating Cylinders. Physical Review Letters 4, 497–498.
  * Donnelly & Ozima (1962) Donnelly, R. J. & Ozima, M. 1962 Experiments on the Stability of Flow Between Rotating Cylinders in the Presence of a Magnetic Field. Proceedings of the Royal Society of London Series A 266, 272–286.
  * Edmonds (1958) Edmonds, Jr., F. N. 1958 Hydromagnetic Stability of a Conducting Fluid in a Circular Magnetic Field. Physics of Fluids 1, 30–41.
  * Goodman & Ji (2002) Goodman, J. & Ji, H. 2002 Magnetorotational instability of dissipative Couette flow. Journal of Fluid Mechanics 462, 365–382.
  * Hollerbach & Rüdiger (2005) Hollerbach, R. & Rüdiger, G. 2005 New Type of Magnetorotational Instability in Cylindrical Taylor-Couette Flow. Physical Review Letters 95 (12), 124501.
  * Hollerbach et al. (2010) Hollerbach, R., Teeluck, V. & Rüdiger, G. 2010 Nonaxisymmetric Magnetorotational Instabilities in Cylindrical Taylor-Couette Flow. Physical Review Letters 104 (4), 044502.
  * Ji et al. (2001) Ji, H., Goodman, J. & Kageyama, A. 2001 Magnetorotational instability in a rotating liquid metal annulus. Month. Not. Roy. Astr. Soc.  325, L1–L5.
  * Kitchatinov & Rüdiger (2004) Kitchatinov, L. L. & Rüdiger, G. 2004 Seed fields for galactic dynamos by the magnetorotational instability.  Astronomy & Astrophysics  424, 565–570.
  * Ogilvie & Pringle (1996) Ogilvie, G. I. & Pringle, J. E. 1996 The non-axisymmetric instability of a cylindrical shear flow containing an azimuthal magnetic field. Month. Not. Roy. Astr. Soc.  279, 152–164.
  * Rädler (1986) Rädler, K. H. 1986 Effect of Differential Rotation on Magnetic Fields of Cosmical Bodies. In Plasma Astrophysics (ed. T. Duc Guyenne & L. M. Zeleny), ESA Special Publication, vol. 251, p. 569.
  * Rüdiger et al. (2018a) Rüdiger, G., Gellert, M., Hollerbach, R., Schultz, M. & Stefani, F. 2018a Stability and instability of hydromagnetic Taylor-Couette flows. Physics Reports  741, 1–89, arXiv: 1703.09919.
  * Rüdiger et al. (2006) Rüdiger, G., Hollerbach, R., Stefani, F., Gundrum, T., Gerbeth, G. & Rosner, R. 2006 The Traveling-Wave MRI in Cylindrical Taylor-Couette Flow: Comparing Wavelengths and Speeds in Theory and Experiment.  The Astrophysical Journal  649, L145–L147.
  * Rüdiger et al. (2003) Rüdiger, G., Schultz, M. & Shalybkov, D. 2003 Linear magnetohydrodynamic Taylor-Couette instability for liquid sodium.  Physical Review E  67 (4), 046312.
  * Rüdiger et al. (2018b) Rüdiger, G., Schultz, M., Stefani, F. & Hollerbach, R. 2018b Magnetorotational instability in Taylor-Couette flows between cylinders with finite electrical conductivity. Geophysical and Astrophysical Fluid Dynamics 112 (4), 301–320, arXiv: 1804.01501.
  * Rüdiger & Shalybkov (2002) Rüdiger, G. & Shalybkov, D. 2002 Stability of axisymmetric Taylor-Couette flow in hydromagnetics.  Physical Review E  66 (1), 016307.
  * Rüdiger & Zhang (2001) Rüdiger, G. & Zhang, Y. 2001 MHD instability in differentially-rotating cylindric flows.  Astronomy & Astrophysics  378, 302–308.
  * Schartman et al. (2009) Schartman, E., Ji, H. & Burin, M. J. 2009 Development of a Couette-Taylor flow device with active minimization of secondary circulation. Review of Scientific Instruments 80 (2), 024501–024501.
  * Schartman et al. (2012) Schartman, E., Ji, H., Burin, M. J. & Goodman, J. 2012 Stability of quasi-Keplerian shear flow in a laboratory experiment.  Astronomy & Astrophysics  543, A94.
  * Seilmayer et al. (2014) Seilmayer, M., Galindo, V., Gerbeth, G., Gundrum, T., 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. Physical Review Letters 113 (2), 024505.
  * Stefani et al. (2006) Stefani, F., Gundrum, T., 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. Physical Review Letters 97 (18), 184502\.
  * Velikhov (1959) Velikhov, E. 1959 Stability of an ideally conducting liquid flowing between cylinders rotating in a magnetic field. Soviet. Phys. JETP 36, 1389–1404.
  * Wang et al. (2022a) Wang, Y., Gilson, E. P., Ebrahimi, F., Goodman, J., Caspary, K. J., Winarto, H. W. & Ji, H. 2022a Identification of a non-axisymmetric mode in laboratory experiments searching for standard magnetorotational instability. Nature Communications 13, 4679, arXiv: 2209.08410.
  * Wang et al. (2022b) Wang, Y., Gilson, E. P., Ebrahimi, F., Goodman, J. & Ji, H. 2022b Observation of Axisymmetric Standard Magnetorotational Instability in the Laboratory. Physical Review Letters  129 (11), 115001, arXiv: 2209.08457.