# Grid-Based Simulations of Polar Circumbinary Disks: Polar Alignment and
Vortex Formation

Ian Rabago1,2 , Zhaohuan Zhu1,2 , Rebecca G. Martin1,2 and Stephen H. Lubow3
1Department of Physics and Astronomy, University of Nevada, Las Vegas 4505 S.
Maryland Parkway Las Vegas, NV 89154, USA
2Nevada Center for Astrophysics, University of Nevada, Las Vegas 4505 S.
Maryland Parkway Las Vegas, NV 89154, USA
3Space Telescope Science Institute, Baltimore, MD 21218
0000-0001-5008-2794 0000-0003-3616-6822 0000-0003-2401-7168
0000-0002-4636-7348

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

###### Abstract

We describe the first grid-based simulations of the polar alignment of a
circumbinary disk. We simulate the evolution of an inclined disk around an
eccentric binary using the grid-based code ATHENA++. The use of a grid-based
numerical code allows us to explore lower disk viscosities than have been
examined in previous studies. We find that the disk aligns to a polar
orientation when the $\alpha$ viscosity is high, while disks with lower
viscosity nodally precess with little alignment over 1000 binary orbital
periods. The timescales for polar alignment and disk precession are compared
as a function of disk viscosity, and are found to be in agreement with
previous studies. At very low disk viscosities (e.g. $\alpha=10^{-5}$),
anticyclonic vortices are observed along the inner edge of the disk. These
vortices can persist for thousands of binary orbits, creating azimuthally
localized overdensities as well as multiple pairs of spiral arms. The vortex
is formed at $\sim 3-4$ times the binary semi-major axis, close to the inner
edge of the disk, and orbits at roughly the local Keplerian speed. The
presence of a vortex in the disk may play an important role in the evolution
of circumbinary systems, such as driving episodic accretion and accelerating
the formation of polar circumbinary planets.

###### keywords:

accretion, accretion discs – protoplanetary discs – hydrodynamics – binaries:
general – planets and satellites: formation – methods: numerical

††pubyear: 20XX††pagerange: Grid-Based Simulations of Polar Circumbinary
Disks: Polar Alignment and Vortex Formation–11

##  1 Introduction

Protoplanetary disks are produced during the star formation process. Stars are
frequently formed in pairs or clusters during the collapse of a molecular
cloud (Duquennoy & Mayor, 1991; Kouwenhoven et al., 2007), and so disks around
multiple stars are expected to be common (Monin et al., 2007). Analogous to
how circumstellar protoplanetary disks can give rise to planets, the
circumstellar and circumbinary disks in a binary star system may act as the
birthsites of the S-type and P-type circumbinary planets, respectively
(Dvorak, 1986). Observations of binary star systems at different evolutionary
stages have revealed different types of protoplanetary disks and planetary
systems, including circumbinary disks (Simon & Guilloteau, 1992; Czekala et
al., 2017; Kennedy et al., 2019; Bi et al., 2020; Kraus et al., 2020),
individual circumstellar disks (Cruz-Sáenz de Miera et al., 2019; Keppler et
al., 2020), and exoplanet systems in both S- and P-type orbits (e.g. Doyle et
al. 2011; Martin 2018).

Binary and higher-order star systems exhibit rich dynamics that are not seen
in single-star systems. The gravitational torque from the inner binary creates
a region of dynamical instability, limiting the range of stable orbits for S-
and P-type planets (Holman & Wiegert, 1999; Quarles et al., 2018; Chen et al.,
2020). Companion stars on faraway orbits can perturb the system over long
timescales, causing orbiting particles to undergo von Zeipel-Lidov-Kozai
oscillations (von Zeipel, 1910; Kozai, 1962; Lidov, 1962). For simplicity, we
restrict most of the discussion in this paper to the binary case.

For cicumbinary or P-type orbits around a circular orbit binary, misalignment
of the orbit to the binary orbital plane causes the orbit to nodally precess
about the binary angular momentum vector. If the inner binary is eccentric, a
second stable configuration exists in which the orbiting body precesses about
the binary eccentricity vector in a polar orientation (Verrier & Evans, 2009;
Farago & Laskar, 2010; Doolin & Blundell, 2011). For circumbinary disks around
the eccentric binary, inclination damping due to viscous forces causes the
disk to settle at either $0\degree$ or $90\degree$, depending on the initial
misalignment. This can lead to the creation of polar disks, where the disk is
aligned perpendicular to the orbits of the binary stars. The stability of
these two orientations has previously been confirmed using analytic theory and
SPH numerical simulations (Martin & Lubow, 2017, 2018; Lubow & Martin, 2018;
Zanazzi & Lai, 2018; Smallwood et al., 2019).

A handful of such “polar disks” have been observed. The HD 98800 system is a
hierarchical double binary system with a circumbinary disk observed around the
B binary component (Kennedy et al., 2019). The A and B binary systems orbit
each other on a $67$ AU orbit with moderate eccentricity $(e\sim 0.5)$, while
the disk-hosting BaBb binary components orbit on a $1$ AU orbit with high
eccentricity $(e\sim 0.8)$ (Zúñiga-Fernández et al., 2021). The HD 98800 B
disk is truncated on both edges by the inner and outer binary, with the dust
disk extending between 2.5 and 4.6 AU (Kennedy et al., 2019). More recently, a
polar circumbinary disk was observed around V773 Tau B (Kenworthy et al.,
2022). Observations of the IRS 43 system (Brinch et al., 2016) have found an
edge-on circumbinary disk with the binary star components orbiting outside of
the circumbinary disk plane. Both components of the binary are also surrounded
by circumstellar disks, each at a different angle to the circumbinary disk,
suggesting complex accretion from the circumbinary disk onto the central
stars.

To date, all simulations of polar disks have been performed using smoothed-
particle hydrodynamics (SPH) (e.g., Martin & Lubow 2017; Kennedy et al. 2019;
Smallwood et al. 2019; Cuello & Giuppone 2019; Martin et al. 2022). These
simulations typically make use of over $10^{6}$ particles and are able to
capture the global dynamical evolution of the disk. However, the range of
viscosities that can be included in SPH simulations is limited by the number
of particles in a given region. The minimum disk effective viscosity (Shakura
& Sunyaev, 1973) is roughly $\alpha=0.01$ for most simulations (Price et al.,
2018). Protoplanetary disks around single stars have been found with
viscosities of $\alpha\lesssim 10^{-4}$ (Pinte et al., 2016; Villenave et al.,
2022), and so the study of how polar disks behave at lower viscosities is
important but has been unexplored.

In this paper we work to expand on previous results using a grid-based
hydrodynamic code to simulate the circumbinary disk. This allows us to extend
the results of previous works to lower viscosities than have previously been
examined using SPH simulations. Our paper is organized as follows. In Section
2, we outline our computational setup. Our results are presented in Section 3.
We discuss the implications to circumbinary disk evolution in Section 4, and
conclude in Section 5.

##  2 Methods

Figure 1:  Visualization of the $\alpha=10^{-2}$ simulation setup at $t=1000$
binary orbits. The disk starts inclined at $t=0$ by rotating the disk along
the x-axis, keeping $\overset{\rightharpoonup}{L}_{\rm disk}$ in the yz-plane,
and undergoes nodal precession about the z-axis due to the binary torque.
Colors show a slice of density along the xz-plane, with red contour showing a
3D isosurface. The orbits of the binary stars are visible in the center. The
simulation grid is shown at half resolution on the right half of the
simulation domain. Regions close to the binary and the z-axis are not in the
computational domain. An animation of the disk evolution is available in the
online material associated with this paper and can be downloaded at
https://youtu.be/Ny-ggFHALMA.  Figure 2: Schematic of the central binary and
circumbinary disk, as well as the angles used to measure the disk orientation.
The binary orbits in the $xz$-plane, with the angular momentum vector
$\overset{\rightharpoonup}{L}_{\rm b}$ pointing along the positive y-axis and
the eccentricity vector $\overset{\rightharpoonup}{e}_{\rm b}$ pointing along
the z-axis. The disk is oriented in 3D space by its angular momentum vector
$\overset{\rightharpoonup}{L}_{\rm disk}$, where it forms an angle $i$ with
the binary angular momentum vector $\overset{\rightharpoonup}{L}_{\rm b}$. The
gray dashed line denotes the disk’s “line of nodes”, where the disk crosses
the $xz$-plane at an angle $\mathbf{\Omega}$ with the eccentricity vector
$\overset{\rightharpoonup}{e}_{\rm b}$.

We simulate a circumbinary disk using the grid-based hydrodynamical code
ATHENA++ (Stone et al., 2020), using spherical-polar coordinates
$(r,\theta,\phi)$ for the simulations. Figure 1 shows a visualization of our
simulation setup. A disk that begins with an inclination in the range for
polar alignment evolution nodally precesses about the binary eccentricity
vector. Therefore, we place the binary in the $xz$-plane of the simulation, so
that the binary angular momentum vector $\overset{\rightharpoonup}{L}_{\rm b}$
points along the positive y-axis and the eccentricity vector
$\overset{\rightharpoonup}{e}_{\rm b}$ points along the positive z-axis. With
this placement of the binary, a disk aligning towards a “polar” orientation
settles towards the $xy$-plane of the simulation domain $(\theta=\pi/2)$. This
binary orientation keeps the disk from precessing outside of the specified
domain during the simulation. This orientation has the additional benefit of
orienting the major axes of the binary towards the poles of the spherical-
polar simulation domain, which prevents strong gravitational interactions from
occurring close to the inner radial boundary.

The primary simulation domain spans a range of $r=[1.07a_{b},10.7a_{b}]$,
$\theta=[5\degree,175\degree]$, and $\phi=[0,2\pi]$. We use a grid of 144
logarithmically spaced cells in $r$, 176 uniformly spaced cells in $\theta$,
and 384 uniformly spaced cells in $\phi$. Our disk setup is based on the setup
used in Martin & Lubow (2017). We initialize the midplane density profile of
the disk according to the density power-law profile

$\rho(R,z=0)=\rho_{0}\left(\frac{R}{R_{0}}\right)^{-p}$ (1)

and set the vertical profile by numerically integrating the density at each
grid cell to establish hydrostatic equilibrium according to the disk scale
height $H=c_{s}/\Omega_{K}$. Here, $R$ is the cylindrical radius to the disk
axis of rotation, $R_{0}=a_{b}$, $\rho_{0}=1$ at $R_{0}$, $c_{s}$ is the local
isothermal sound speed $\sqrt{P/\rho}$ at $R$, and $\Omega_{K}$ is the
Keplerian orbital frequency at $R$. The disk temperature is initialized using
the power-law profile

$T(r)=T_{0}\left(\frac{r}{r_{0}}\right)^{-q},$ (2)

where $r_{0}=a_{b}$ and $T_{0}=1$ at $R_{0}$. We use power-law exponents of
$p=2.25$ and $q=1.5$, corresponding to an initial disk surface density profile
$\Sigma\propto r^{-3/2}$ and a constant kinematic viscosity $\nu=\alpha c_{\rm
s}^{2}/\Omega_{K}$. The density profile is truncated outside of the range
$[2a_{b},5a_{b}]$ using an exponential cutoff of the form
$\exp{\left(-|r-\mu_{r}|/\sigma_{r}\right)}$, where $\mu_{r}$ are the disk
edges of $2a_{b}$ and $5a_{b}$ and the relative scale length $\sigma_{r}$ is
$\sigma_{r,\rm in}=0.35a_{b}$ for the inner edge and $\sigma_{r,\rm
out}=1.8a_{b}$ for the outer edge.

We initialize the disk with a spherically isothermal temperature profile, as
well as an initial scale height of $H/r=0.1$ at $R=3.5a_{b}$. We use the
orbital cooling scheme outlined in Equation 5 of Zhu et al. (2015), with a
dimensionless cooling time of $t_{\rm cool}=0.01\Omega_{K}$. At the simulation
boundaries, we instate a one-way outflow boundary condition in the radial
direction and a reflecting boundary condition in the polar direction. We use a
spherically symmetric density floor with a value of $\rho_{\rm
floor}=10^{-4}\rho_{0}$ at $r=1$ and a power-law slope of $d=2.25$, identical
to that of the density profile. The disk velocity is initialized with a
Keplerian profile using the total mass of the binary $M_{\rm
tot}=m_{1}+m_{2}$. Since the axis of rotation is now within the simulation
domain, regions near the poles of rotation are initialized with a velocity

$\mathrm{v}=\sqrt{\frac{GM_{\rm tot}}{r}}\frac{R}{r},$ (3)

which prevents excessively high velocities occurring near the poles of
rotation.

The binary components are simulated as gravitational bodies, with equal masses
$m_{1}=m_{2}=0.5$, which we place in an eccentric orbit with semi-major axis
$a_{b}$ and eccentricity $e_{b}$. We choose an eccentricity of $e_{b}=0.5$.
This sets the initial tilt required for disk polar alignment of a low mass
disc to be $i_{\rm crit}<i<180\degree-i_{\rm crit}$, where $i_{\rm
crit}=\arcsin{\sqrt{3/8}}\simeq 38\degree$, based on the behavior of test
particles (Farago & Laskar, 2010). For massive disks, the critical angle
changes due to angular momentum exchange between the disk and the binary
(Martin & Lubow, 2019). Since the binary does not feel the gravitational force
from the disk and disk self-gravity is not included from our simulations, the
disk can be considered to be in the low mass regime. We use a 2nd-order
leapfrog integrator to solve the equations of motion for the binary. We use
the CFL fluid timestep of ATHENA++ as the timestep for the integrator.

Figure 2 shows the geometric setup of the binary-disk system, with inclination
angle $i$ and ascending node $\mathbf{\Omega}$ labeled. The angles $i$ and
$\mathbf{\Omega}$ are measured relative to the binary plane using the
equations

$i=\arccos{\left(\frac{L_{y}}{L}\right)}$ (4)

and

$\mathbf{\Omega}=\arctan{\left(\frac{L_{z}}{L_{x}}\right)},$ (5)

where $L_{x}$, $L_{y}$ and $L_{z}$ are the components of the disk angular
momentum $L_{\rm disk}$ in the simulation $x$, $y$ and $z$ axes.

We initialize the disk at an inclination of $i=120\degree$ and ascending node
$\mathbf{\Omega}=90^{\circ}$ with respect to the plane of the binary by
converting to tilted disk coordinates (Zhu 2019, Eq. 44, also see Appendix A).
The high initial inclination ensures the disk starts within the librating
region that will evolve towards a polar alignment. For test particles around
an equal-mass binary, there is no difference between the choice of prograde
$(60\degree)$ versus retrograde $(120\degree)$ orbits, as the evolution of the
test particle only depends on whether the angle between the binary and disk
planes is beyond the critical inclination angle $i_{\rm crit}$. When full
three-body systems are considered, angular momentum can be exchanged between
the binary system and the circumbinary particle (Farago & Laskar, 2010; Martin
& Lubow, 2019). This creates a difference between the prograde and retrograde
cases, and large angular momentum ratios $j=L_{\rm disk}/L_{\rm b}$ can lead
to the formation of “crescent” librating orbits (Chen et al., 2019; Abod et
al., 2022). Since we do not model the exchange of angular momentum between the
disk and the binary, we expect the initially prograde and retrograde cases to
be similar.

We simulate the disk evolution for 1000 binary orbits, using the vl2 second-
order van Leer predictor-corrector time integrator (Stone & Gardiner, 2009)
and third-order PPM reconstruction. We perform a set of five simulations,
varying the $\alpha$-viscosity parameter between $\alpha=10^{-1},5\times
10^{-2},10^{-2},10^{-3}$, and $10^{-5}$. The disk viscosity is constant
throughout most of the simulation domain. Close to the inner radial boundary,
where $R<R_{\rm{min}}=1.07a_{b}$, we reduce the value of the kinematic
viscosity $\nu$ using an exponential cutoff

$\nu=\alpha\frac{c_{\rm s}^{2}}{\Omega_{K}}\exp{\left(\frac{R-R_{\rm
min}}{0.268}\right)}.$ (6)

This prevents the diffusive timestep from becoming too small and restricting
the simulation speed. For simulations with $\alpha>0.01$, the viscous
timescale $t_{\rm visc}\sim R^{2}/\nu$ is short enough that material spreading
of the disk requires that we expand our simulation region. For these
simulations, we expand the outer radial edge of our simulation domain to
$R_{\rm out}=35.7a_{b}$ using 216 logarithmically spaced radial cells in order
to match the increased viscous spreading of the disk.

##  3 Results

Figure 3: Time evolution of the disk for different values of $\alpha$, showing
large changes in orientation. _Top:_ Disk inclination $i$. _Top Middle:_
Longitude of ascending node $\mathbf{\Omega}$. _Middle:_ Angle between the
disk and the binary eccentricity vector $\theta_{\rm polar}$. _Bottom Middle:_
Longitude of ascending node in the xy-plane $\mathbf{\Omega}_{\rm polar}$.
_Bottom:_ Disk mass as a fraction of the initial mass. All angular quantities
are measured in reference to the binary orbital plane at a distance of
$R=3a_{b}$.  Figure 4: Trajectories of the simulations in
$i\mathrm{\cos\mathbf{\Omega}}-i\mathrm{sin\mathbf{\Omega}}$ phase space. Each
trajectory starts at the black ‘$\times$’ and spirals inwards counterclockwise
towards the point $(0.0,0.5)$, with the final state of each simulation at
$1000$ binary orbits marked by the colored dots. Angular quantities are
measured in reference to the binary orbital plane at a distance of $R=3a_{b}$,
as in Figure 3. Gray paths show the trajectories of test particle orbits,
spaced by $10\degree$. An animation of the test particle precession, as well
as the projection onto the phase space, is available with the online version
of this paper and can be downloaded at https://youtu.be/ZBkhG6pkOAQ.

###  3.1 Polar Alignment

We observe polar libration of the disk in all of our simulations. The top two
panels of Figure 3 show the disk inclination and ascending node oscillating
about $i=\mathbf{\Omega}=90\degree$, with a precession time $t_{\rm prec}$
that increases with disk viscosity. The third and fourth panels of Figure 3
show the polar libration of the disk in reference to the simulation xy-plane,
using the angles $\theta_{\rm polar}$ and $\mathbf{\Omega}_{\rm polar}$, the
angle between $\overset{\rightharpoonup}{L}_{\rm disk}$ and
$\overset{\rightharpoonup}{e}_{\rm b}$ and longitude of ascending node in the
xy-plane , respectively. These angles are calculated as

$i_{\rm polar}=\arccos{\left(\frac{L_{z}}{L}\right)}$ (7)

and

$\mathbf{\Omega}_{\rm polar}=\arctan{\left(\frac{L_{y}}{L_{x}}\right)}.$ (8)

Precession of the disk causes $\mathbf{\Omega}_{\rm polar}$ to circulate. As
the disk oscillates, inclination damping causes the oscillation amplitude to
shrink over time, allowing the disk to settle towards a polar orientation and
$\theta_{\rm polar}$ to settle towards 0. For fixed disk density and
temperature distributions, the timescale for polar alignment $t_{\rm damp}$ is
inversely proportional to $\alpha$ (King et al., 2013; Lubow & Martin, 2018),
and so disks with a higher $\alpha$-viscosity settle towards a polar alignment
quicker, performing less oscillations before reaching a polar configuration.

From Figure 3, the polar alignment timescale can be estimated to be roughly
$300$ binary orbits for the $\alpha=0.1$ simulation, and roughly $500$ binary
orbits for the $\alpha=0.05$ simulation. This is roughly consistent with the
analytical predictions given by Equations 29-30 of Lubow & Martin (2018),
which give $t_{\rm damp}=220\ T_{b}$ and $455\ T_{b}$, respectively. The
difference in density evolution of these two cases likely accounts in part for
how $t_{\rm damp}$ departs from a pure $1/\alpha$ dependence as $\alpha$
changes. For simulations with viscosities $\alpha<0.05$, the oscillation
amplitude is not reduced significantly during the course of our simulation
time, suggesting that the timescale for polar alignment is at least on the
order of several thousand binary orbits. The polar alignment timescales for
the $\alpha=0.01$ simulation are also similar to the SPH simulations with high
resolution ($10^{6}$ particles). SPH simulations at lower resolution show a
quicker dampening of the tilt oscillations and thus a lower alignment
timescale.

We measure the precession periods in our $\alpha=0.1$ and $\alpha=0.01$
simulations to be roughly $300-400$ binary orbits. The precession period
increases over time as the disk material is redistributed outwards. These
values are in rough agreement with the SPH simulations of the same Shakura-
Sunyaev viscosity111The equivalent SPH artificial viscosity for these
simulations is $\alpha_{\rm{SPH}}=4$ and $0.4$ respectively (Lodato & Price,
2010). in Martin & Lubow (2017) and Martin & Lubow (2018) (see their Eq. 2),
which measure precession periods of roughly $300$ binary orbits for disks that
start with a tilt of $80^{\circ}$.

The settling behavior towards polar alignment is shown clearly in Figure 4,
which plots the trajectory of the disk in the
$i\,\mathrm{\cos\mathbf{\Omega}}-i\,\mathrm{sin\mathbf{\Omega}}$ phase space,
along with trajectories of test particles in gray. After an initial damping
phase, present in all simulations, the disk trajectories spiral inwards
towards the point $(0.0,0.5)$, corresponding to a disk oriented at exactly
$90\degree$ inclination. High-viscosity disks $(\alpha\gtrsim 0.05)$ spiral
inwards and align towards polar orientation quickly. Disks with lower
viscosities $(\alpha<0.05)$ follow nearly closed oscillating trajectories
close to the $70\degree$ test particle trajectory, with lower values of
$\alpha$ settling towards trajectories slightly closer to a polar orientation.
As noted above, the SPH simulations that start at an inclination of
$80^{\circ}$ with $\alpha=0.01$ show a significant decay of tilt over 1000
binary orbital periods.

Figure 5: _Top:_ Spherically integrated radial surface density profiles at
$t=1000$ binary orbits, plotted on a logarithmic scale. _Bottom:_ Radially
averaged disk eccentricity at $t=1000$ binary orbits, calculated with Equation
16 of Shi et al. (2012).  Figure 6: Face-on density profiles of the disk at
$t=675$ binary orbits. Figure 7: Closeup of the vortex in the $\alpha=10^{-5}$
simulation at $t=692$ binary orbits. _Top left:_ Gas midplane density. The
local overdensity and spiral arm pairs are clearly visible. _Top right:_ Local
vorticity. Strong anticyclonic regions are indicated in blue. Spiral arms are
visible as red diagonal lines in the vicinity of the vortex. In both plots,
black curves show contours of density. _Middle:_ Azimuthal density across the
vortex center. _Bottom:_ Radial density across the vortex center.

Figure 5 shows the disk surface density profiles at $t=1000$ binary orbits.
The inner edge of the disk drops off quickly in all simulations at roughly
$2a$, with higher viscosity disks exhibiting a smaller inner cavity. This
behavior is in agreement with previous numerical simulations (Franchini et
al., 2019), as well as analytical predictions (Miranda & Lai, 2015; Lubow &
Martin, 2018) suggesting that disks in a polar orientation should be truncated
somewhat closer to the binary than in the coplanar case. The model of Miranda
& Lai (2015) suggests that the disks in these simulations should be should be
truncated at the 1:3 commensurability for an outer Lindblad resonance. Streams
of gas that flow into the central gap may carry much or all of the gas
required for a steady-state accretion disk, as is known in the coplanar case
(e.g., Artymowicz & Lubow, 1996; Shi et al., 2012; Muñoz & Lai, 2016). The
profiles generally follow a wider, thinner distribution with increasing
$\alpha$-viscosity. Since the diffusion coefficient $\nu$ is proportional to
$\alpha$, the disks with a larger $\alpha$ are able to spread out further due
to their shorter viscous timescale $t_{\rm visc}\sim R^{2}/\nu$. As the disk
material spreads radially, the change in mass distribution increases the
disk’s moment of inertia, increasing the precession timescale as seen in
Figure 3. Notably, the case of $\alpha=10^{-5}$ does not follow this trend,
having a slightly wider surface density profile than the $\alpha=10^{-3}$
simulation. This may be due to the formation of vortices at low values of
$\alpha$, which generates spiral wakes in the disk. These spirals can increase
the outward transport of material and widen the surface density profile. We
discuss the details of vortex formation and their impacts on the disk in
Section 3.2.

The bottom panel of Figure 5 shows the shell-averaged eccentricity of the disk
at $t=1000$ binary orbits. We calculate the eccentricity of the disk using
Equation 16 of Shi et al. (2012). The disk eccentricity increases and
eventually saturates as time increases, so these values can be taken as a
maximum value for each simulation. We find that simulations with
$\alpha\gtrsim 0.01$ are able to maintain low eccentricities, while disks with
lower viscosities are excited to moderate $(e\sim 0.1)$ eccentricity.

Figure 8: Snapshots of the $\alpha=10^{-5}$ simulation, showing the local disk
vorticity. In the left four panels, blue regions denote areas of large
anticylonic vorticity. The black contours highlight high-density regions
$(\rho=0.6\rho_{\rm max})$ within the disk midplane. _Far Left:_ Early on, the
RWI triggers with an $m=2$ mode, creating two vortices on opposite sides of
the disk. _Left:_ Differences in orbital velocities cause the vortices to move
together and merge. _Middle:_ The result is a single vortex which orbits at
the inner edge of the disk. _Right:_ The remaining vortex is long-lived, and
persists until the end of the simulation. _Far Right:_ Midplane disk density,
showing the vortex as well as the multiple spiral arms. The density scale is
shifted compared to Figure 7 in order to highlight the spiral arms.

###  3.2 Vortex Formation

Figure 9: Orbital motion of the vortex in the $\alpha=10^{-5}$ simulation.
_Top:_ Orbital radius of the vortex. _Middle:_ Orbital frequency, in units of
the binary orbital frequency $\Omega_{b}$. The dotted line shows the Keplerian
orbital speed at the vortex’s orbital radius. _Bottom:_ Vortex aspect ratio
$\chi=r\Delta\phi/\Delta r$.

A gallery of disk midplane density profiles is shown in Figure 6. Each panel
shows a face-on view of the disk at $t=675$ binary orbits. Similar to Figure
5, material in disks with higher viscosity is more spread out with a lower
surface density. Some light spiral features are present towards the inner edge
of the disk in all simulations, which vary in strength as the disk precesses
and changes its orientation relative to the binary, but most of the disks are
azimuthally symmetric and largely featureless.

The notable exception to this behavior is the nearly inviscid $\alpha=10^{-5}$
simulation, which shows a prominent density enhancement along the inner edge
of the disk, as well as an associated one-armed spiral. Horseshoe-shaped
overdensity features have previously been observed in ALMA dust continuum
images (Casassus et al., 2013; van der Marel et al., 2013; Casassus, 2016),
and have been explained by disk material moving on eccentric orbits (Ragusa et
al., 2017), or vortices generated by sharp surface density gradients via the
Rossby Wave Instability (RWI, Lovelace et al. 1999; Li et al. 2000). Our
simulations are the first instances where vortices have been observed within
polar disks. The steep surface density profile along the inner edge of the
disk creates a minimum in the local vortensity, an important characteristic in
generating RWI vortices (Bae et al., 2015). Vortex formation can be suppressed
by small amounts of viscosity $(\alpha\lesssim 10^{-3})$ (de Val-Borro et al.,
2007; Fu et al., 2014; Zhu & Stone, 2014; Owen & Kollmeier, 2017), so they
were not observed in previous SPH simulations of polar disks.

A close-up of the overdensity feature is shown in Figure 7. We calculate the
midplane density and vorticity by transforming the disk data as outlined in
Appendix A. The left panel shows the disk density, where the overdense region
is clearly shown. The right panel shows the local disk vorticity
$(\omega-\omega_{0})/\omega_{0}$, where
$\omega=-(\nabla\times\mathrm{v}_{\theta})$ and
$\omega_{0}=\frac{1}{2}\Omega_{K}$ is the initial Keplerian vorticity. At the
disk midplane after the coordinate transformation, the unit vector in $\theta$
points downwards, opposite in direction to the unit vector in $z$, so we add a
negative sign to our calculation of vorticity to match the sign convention of
vorticity in other coordinate systems, i.e.
$\omega=\nabla\times\mathrm{v}_{z}$. Large areas of anticyclonic motion are
coincident with the overdensity region, and confirm that the overdense region
is associated with the creation of a vortex. The density cuts across the
vortex show that the peak gas density $\rho_{max}$ is enhanced by a factor of
$\gtrsim 8$ times the background, which is consistent with previous
simulations of RWI vortices (e.g. Lyra et al. 2008; Bae et al. 2015).

Two spiral arms are visible on the leading side of the overdensity, extending
along slanted density contours and trailing inwards towards the star. A
similar pair of outward pointing spirals is present along the trailing edge of
the vortex, but are less defined, causing the two spirals to partially combine
and form a single, wide spiral arm behind the vortex. In the vorticity plot,
the spirals are visible in red as diagonal lines of high positive vorticity.
Together, these create two distinct pairs of spiral arms originating from the
vortex: a “central” pair of spirals, which begin along the center line of the
vortex (roughly $3\ a_{\rm b}$ in the figure) and extend directly away from
the vortex, and a “tangential” pair of arms, which are created along the inner
and outer radial edges of the vortex and run tangent to the edge of the
vortex.

The initial triggering of the RWI and creation of the vortex follows the
evolution outlined in Bae et al. (2015). We show panels of the vortex growth
and evolution in Figure 8. The large vortensity minimum at the inner edge of
the disk generates several RWI vortices during the linear growth phase of the
instability (far left). The number of initial vortices can vary, but our
simulations show a consistent $m=2$ mode initially present within the first
tens of binary orbits. The initial appearance of this mode may be related to
the binary components, as previous single-star numerical calculations find the
$m=3-5$ modes to have the fastest growth rates (Li et al., 2000; Lin, 2012).
Once the vortices are formed, differences in their orbital speeds will cause
them to migrate towards each other and merge (left), eventually forming a
singe anticyclonic vortex. The merging process is relatively quick, and all
vortices are combined into a single vortex within a few hundred orbits of
their creation (middle). The remaining vortex persists until the end of the
simulation (right). A comparison with the disk density (far right) shows the
coincident overdensity and spiral arm features as in Figure 7.

The orbital motion of the vortex is shown in Figure 9. We locate cells of
density $\rho>0.5\rho_{max}$ within the vortex and calculate the average of
their locations as the vortex radius and azimuthal position in the local disk
coordinates. We also use the radial and azimuthal extent of these cells to
calculate the aspect ratio of the vortex $\chi=r\Delta\phi/\Delta r$. After
the initial vortices merge, the resulting single vortex orbits at roughly 3
times the binary separation at Keplerian speed. The vortex aspect ratio is
somewhat more variable even after the vortices have finished merging, but is
generally around a value of $\chi\sim 5$.

##  4 Discussion

###  4.1 The effect of vortices vs. lumps in the disk

The overdensity feature observed in the $\alpha=10^{-5}$ simulation is
noticeably different than similar features created in coplanar circumbinary
disks. Several coplanar disks have been found to exhibit horseshoe features,
but these features have been explained as a density “lump”, created by
eccentric gas at the inner edge of the disk interacting with inflow streams
created by the action of the binary torque (Shi et al., 2012; Miranda et al.,
2017; Ragusa et al., 2017, 2020). In contrast, the overdensity feature in the
polar disk corresponds to the vorticity minimum in the disk, and is closer in
nature to vortex-induced clumps caused at the edges of gaps opened by massive
planets (de Val-Borro et al., 2007; Zhu et al., 2014; Hammer et al., 2017).

Overdensities in circumbinary disks can drive periodic accretion onto the
central binary. Simulations of coplanar curcumbinary disks and the binary
components by Muñoz & Lai (2016) show the lump at the inner edge causes pulsed
accretion at the rate of the lump’s orbital period. When binary eccentricity
is included, the accretion flows become pulsed at the rate of the binary
orbital period, and can preferentially favor one star for several hundred
orbits. In Miranda et al. (2017), mass accretion is found to vary around
circular binaries on short timescales of $\frac{1}{2}T_{b}$, as well as on
longer timescales of roughly $5T_{b}$.

Accretion rates in polar disks have been less studied. Smallwood et al. (2022)
perform simulations of the HD 98800 system and find that binary accretion
rates show no periodicity and are roughly constant in time for a polar disk
with $\alpha=0.01$. The high viscosity in these simulations means they do not
form vortices, even during the passage of the outer AaAb binary.

Figure 10: Shell-averaged radial mass flux of the $\alpha=10^{-5}$ simulation.
The dashed line marks the approximate orbital distance of the vortex. Short-
period outflows are visible close to the binary $(2a_{b})$ as thin red
stripes. Long-period variability, associated with the orbital period of the
vortex, is visible from $3-4a_{b}$ as alternating red and blue bands.

The vortex in our $\alpha=10^{-5}$ simulation exhibits similar variable
behavior on both short and long timescales. Figure 10 shows the shell-averaged
radial mass flux over roughly 10 binary orbits. Long-term flows are present at
distances of $3-4a_{b}$, oscillating with a period of roughly $5T_{b}$. This
low frequency mass flux is due to the orbital motion of the vortex around the
binary, carrying the overdense region on an eccentric orbit at its local
orbital period (Figures 5, 9). Closer to the binary, at roughly $2a_{b}$,
short-term outflows are visible with a periodicity of roughly $T_{b}$. These
outflows are caused by small streams of material pulled from the inner edge of
the disk and flung outwards by the binary torque, similar in nature to those
seen by Shi et al. (2012). Since our inner binary is eccentric, these streams
do not pile up and form an overdense lump (Miranda et al., 2017). Overall the
vortex in our polar circumbinary disks is different from the lump in coplanar
circumbinary disks in two ways: 1) the vortex only forms when $\alpha$ is very
low, while the lump can form in disks with a much higher $\alpha$ (even in MHD
simulations with $\alpha\sim 0.1$, Shi et al. 2012); 2) the vortex exists in
polar circumbinary disks with an eccentric binary, while the lump is only
found when the binary is nearly circular ($e_{\rm b}\lesssim 0.05$, Miranda et
al. 2017).

The overdensity created by the vortex changes the local pressure gradient on
the neighboring gas, producing regions of super-Keplerian and sub-Keplerian
velocity just inside and outside of the vortex orbital radius, respectively
(Kuznetsova et al., 2022). These deviations from Keplerian rotation may be
observable as a kinematic signature (Boehler et al., 2021), which are inverted
but similar in strength to the signatures formed in the gaps of protoplanetary
disks thought to be formed by young planets (Teague et al., 2018; Teague et
al., 2019; Izquierdo et al., 2022). Future high-resolution studies with ALMA
may be able to distinguish vortex creation in polar disks.

Observations of HD 98800 B from Kennedy et al. (2019) show a uniform disk in
both gas and dust, with no evidence of large scale disk asymmetry. The absence
of vortex-like features may be a sign that the viscosity in HD 98800 B is high
enough to prevent vortices from forming. Another explanation could be that
vortices are relatively transient features within the disk. The lifetime of
vortices in the disk is on the order of thousands of orbital periods (Zhu &
Stone, 2014; Hammer et al., 2021). For HD 98800 B, this corresponds to a
lifetime of $\lesssim 1000$ years for a vortex formed at the inner edge of the
disk, which is much shorter than the disk lifetime of $\sim 10$ Myr (Barrado Y
Navascués, 2006). Therefore, it is possible that any vortices that initially
formed in the disk have long since dispersed, unless vortices are continuously
generated by the inner binary.

The rightmost panel of Figure 8 shows the vortex and spiral arms as they
appear in the disk. Notably, the appearance of the spiral arms resembles the
scattered light images of HD 142527 (Avenhaus et al., 2014; Avenhaus et al.,
2017), which is also known to have an overdensity feature (Casassus et al.,
2013). Spirals along the inner edge of the disk can be created, without the
creation of a vortex, due to the gravitational interaction of the binary with
the disk (Price et al., 2018), so it is unclear if vortex-generated spirals
are the cause of the spiral arms in HD 142527.

###  4.2 Prospects for Polar Circumbinary Planets

It is unknown if circumbinary planets are able to form within polar disks.
Currently, no planets have been observed in polar orbits around binaries, nor
planet candidates in polar-aligned disks. The lack of detections is at least
partially due to a strong observational bias; current methods for detecting
exoplanets rely on identifying data with multiple, periodic transits, which
favor single-star, coplanar systems. Planets on circumbinary orbits exhibit
large transit timing variations, on the order of hours or days (Armstrong et
al., 2013). Inclined orbits are more likely to produce single or irregular
transits due to the orbital precession induced by the binary torque (Martin &
Triaud, 2014; Chen et al., 2019; Chen et al., 2022), rendering most detection
methods used for single planets unusable and requiring the use of special
transit folding methods (e.g. Martin & Fabrycky 2021). If circumbinary planets
are equally likely to form in polar disks as they are in coplanar disks, they
may be as common as planets around single-star systems (Armstrong et al.,
2014; Martin & Triaud, 2014).

The formation of vortices within low-viscosity disks can have significant
effects on planet formation history. An outstanding problem in planet
formation involves the formation of large planetary embryos in a disk, as gas
drag from the disk can cause growing particles to rapidly spiral inwards onto
the star, before they are able to decouple from the gas (Weidenschilling,
1977). Anticyclonic vortices in a disk act as grain traps, and can intercept
over half of the grains that cross its orbit during their inward radial
migration (Fromang & Nelson, 2005). Once captured in the vortex, the local
vorticity works to focus particles towards the vortex “eye”, which can
facilitate accelerated growth of planetary material all the way up to Jupiter-
mass planets (Klahr & Bodenheimer, 2006; Lyra et al., 2008; Zhu et al., 2014;
Owen & Kollmeier, 2017). In polar disks, vortices may also lead to direct
planet formation or seed the disk with a large initial planetesimal
population, which can allow further planetary growth and formation terrestrial
planets in polar orbits (Childs & Martin, 2021, 2022).

Figures 7 and 9 show that the vortex orbits at a few times the binary
separation, close to the inner edge of the disk. Indeed, the sharp surface
density profile created at the inner cavity wall allows the RWI to produce
vortices in the absence of viscosity. This is close to the theoretical
stability limit for circumbinary planets (Dvorak, 1986). Curiously, many of
the circumbinary planets discovered by _Kepler_ and _TESS_ orbit just outside
the inner stability limit (Martin & Triaud, 2014; Welsh et al., 2014).
Together, this may suggest a potential in-situ formation scenario for close-in
circumbinary planets, where the sharp inner cavity of the circumbinary disk
generates long-lived vortices which allows for efficient trapping of solid
materials. Rapid planet formation at this distance would produce circumbinary
planets close to the circumbinary inner stability limit.

##  5 Conclusion

In this paper, we describe the first grid-based simulations on the evolution
of polar-aligning circumbinary disks. We find the disks align towards a polar
orientation, with timescales in rough agreement with previous SPH studies and
analytic estimates. In nearly inviscid disks, RWI vortices can be generated
close to the disk inner edge while polar alignment occurs. These vortices can
persist for thousands of binary orbits in a polar disk, and create
overdensities within the disk. Two pairs of spiral arms are seen originating
from the vortex, along the center and tangent to the vortex core. The
overdensities created by the vortices are fundamentally different from similar
features created through eccentricity excitation within the disk. The
combination of vortices and spiral arms may create a detectable signature in
polar disks, though none have been observed as of yet. The presence of
vortices and overdensities in polar circumbinary disks may also aid in the
formation of polar circumbinary planets, accelerating planet growth close to
the inner stability limit.

## Acknowledgements

This material is based upon work supported in part by UNLV’s Faculty Top Tier
Doctoral Graduate Research Assistantship Grant Program (TTDGRA) and the
National Aeronautics and Space Administration under Grant No. 80NSSC20M0043
issued through the Nevada NASA Space Grant Consortium. Z. Z. acknowledges
support from the National Science Foundation under CAREER Grant Number
AST-1753168. RGM and SHL acknowledge support from NASA through grants
80NSSC21K0395 and 80NSSC19K0443. Figures in this paper were made with the help
of Matplotlib (Hunter, 2007), NumPy (Harris et al., 2020), and VisIt (Childs
et al., 2012). The animations included in the online version of this paper
were made with the help of the NAS Visualization Team.

## Data Availability

The data used in this paper is available upon request to the corresponding
author.

## References

  * Abod et al. (2022) Abod C. P., Chen C., Smallwood J., Rabago I., Martin R. G., Lubow S. H., 2022, MNRAS, 517, 732
  * Armstrong et al. (2013) Armstrong D., et al., 2013, MNRAS, 434, 3047
  * Armstrong et al. (2014) Armstrong D. J., Osborn H. P., Brown D. J. A., Faedi F., Gómez Maqueo Chew Y., Martin D. V., Pollacco D., Udry S., 2014, MNRAS, 444, 1873
  * Artymowicz & Lubow (1996) Artymowicz P., Lubow S. H., 1996, ApJ, 467, L77
  * Avenhaus et al. (2014) Avenhaus H., Quanz S. P., Schmid H. M., Meyer M. R., Garufi A., Wolf S., Dominik C., 2014, ApJ, 781, 87
  * Avenhaus et al. (2017) Avenhaus H., et al., 2017, AJ, 154, 33
  * Bae et al. (2015) Bae J., Hartmann L., Zhu Z., 2015, ApJ, 805, 15
  * Barrado Y Navascués (2006) Barrado Y Navascués D., 2006, A&A, 459, 511
  * Bi et al. (2020) Bi J., et al., 2020, ApJ, 895, L18
  * Boehler et al. (2021) Boehler Y., et al., 2021, A&A, 650, A59
  * Brinch et al. (2016) Brinch C., Jørgensen J. K., Hogerheijde M. R., Nelson R. P., Gressel O., 2016, ApJ, 830, L16
  * Casassus (2016) Casassus S., 2016, Publ. Astron. Soc. Australia, 33, e013
  * Casassus et al. (2013) Casassus S., et al., 2013, Nature, 493, 191
  * Chen et al. (2019) Chen C., Franchini A., Lubow S. H., Martin R. G., 2019, MNRAS, 490, 5634
  * Chen et al. (2020) Chen C., Lubow S. H., Martin R. G., 2020, MNRAS, 494, 4645
  * Chen et al. (2022) Chen C., Lubow S. H., Martin R. G., 2022, MNRAS, 510, 351
  * Childs & Martin (2021) Childs A. C., Martin R. G., 2021, ApJ, 920, L8
  * Childs & Martin (2022) Childs A. C., Martin R. G., 2022, ApJ, 927, L7
  * Childs et al. (2012) Childs H., et al., 2012, in , High Performance Visualization–Enabling Extreme-Scale Scientific Insight. pp 357–372
  * Cruz-Sáenz de Miera et al. (2019) Cruz-Sáenz de Miera F., Kóspál Á., Ábrahám P., Liu H. B., Takami M., 2019, ApJ, 882, L4
  * Cuello & Giuppone (2019) Cuello N., Giuppone C. A., 2019, A&A, 628, A119
  * Czekala et al. (2017) Czekala I., et al., 2017, ApJ, 851, 132
  * Doolin & Blundell (2011) Doolin S., Blundell K. M., 2011, MNRAS, 418, 2656
  * Doyle et al. (2011) Doyle L. R., et al., 2011, Science, 333, 1602
  * Duquennoy & Mayor (1991) Duquennoy A., Mayor M., 1991, A&A, 500, 337
  * Dvorak (1986) Dvorak R., 1986, A&A, 167, 379
  * Farago & Laskar (2010) Farago F., Laskar J., 2010, MNRAS, 401, 1189
  * Franchini et al. (2019) Franchini A., Lubow S. H., Martin R. G., 2019, ApJ, 880, L18
  * Fromang & Nelson (2005) Fromang S., Nelson R. P., 2005, MNRAS, 364, L81
  * Fu et al. (2014) Fu W., Li H., Lubow S., Li S., 2014, ApJ, 788, L41
  * Hammer et al. (2017) Hammer M., Kratter K. M., Lin M.-K., 2017, MNRAS, 466, 3533
  * Hammer et al. (2021) Hammer M., Lin M.-K., Kratter K. M., Pinilla P., 2021, MNRAS, 504, 3963
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Holman & Wiegert (1999) Holman M. J., Wiegert P. A., 1999, AJ, 117, 621
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Izquierdo et al. (2022) Izquierdo A. F., Facchini S., Rosotti G. P., van Dishoeck E. F., Testi L., 2022, ApJ, 928, 2
  * Kennedy et al. (2019) Kennedy G. M., et al., 2019, Nature Astronomy, 3, 230
  * Kenworthy et al. (2022) Kenworthy M. A., et al., 2022, arXiv e-prints, p. arXiv:2207.05575
  * Keppler et al. (2020) Keppler M., et al., 2020, A&A, 639, A62
  * King et al. (2013) King A. R., Livio M., Lubow S. H., Pringle J. E., 2013, MNRAS, 431, 2655
  * Klahr & Bodenheimer (2006) Klahr H., Bodenheimer P., 2006, ApJ, 639, 432
  * Kouwenhoven et al. (2007) Kouwenhoven M. B. N., Brown A. G. A., Portegies Zwart S. F., Kaper L., 2007, A&A, 474, 77
  * Kozai (1962) Kozai Y., 1962, AJ, 67, 591
  * Kraus et al. (2020) Kraus S., et al., 2020, Science, 369, 1233
  * Kuznetsova et al. (2022) Kuznetsova A., Bae J., Hartmann L., Low M.-M. M., 2022, ApJ, 928, 92
  * Li et al. (2000) Li H., Finn J. M., Lovelace R. V. E., Colgate S. A., 2000, ApJ, 533, 1023
  * Lidov (1962) Lidov M. L., 1962, Planet. Space Sci., 9, 719
  * Lin (2012) Lin M.-K., 2012, ApJ, 754, 21
  * Lodato & Price (2010) Lodato G., Price D. J., 2010, MNRAS, 405, 1212
  * Lovelace et al. (1999) Lovelace R. V. E., Li H., Colgate S. A., Nelson A. F., 1999, ApJ, 513, 805
  * Lubow & Martin (2018) Lubow S. H., Martin R. G., 2018, MNRAS, 473, 3733
  * Lyra et al. (2008) Lyra W., Johansen A., Klahr H., Piskunov N., 2008, A&A, 491, L41
  * Martin (2018) Martin D. V., 2018, in Deeg H. J., Belmonte J. A., eds, , Handbook of Exoplanets. p. 156, doi:10.1007/978-3-319-55333-7_156
  * Martin & Fabrycky (2021) Martin D. V., Fabrycky D. C., 2021, AJ, 162, 84
  * Martin & Lubow (2017) Martin R. G., Lubow S. H., 2017, ApJ, 835, L28
  * Martin & Lubow (2018) Martin R. G., Lubow S. H., 2018, MNRAS, 479, 1297
  * Martin & Lubow (2019) Martin R. G., Lubow S. H., 2019, MNRAS, 490, 1332
  * Martin & Triaud (2014) Martin D. V., Triaud A. H. M. J., 2014, A&A, 570, A91
  * Martin et al. (2022) Martin R. G., Lepp S., Lubow S. H., Kenworthy M. A., Kennedy G. M., Vallet D., 2022, ApJ, 927, L26
  * Miranda & Lai (2015) Miranda R., Lai D., 2015, MNRAS, 452, 2396
  * Miranda et al. (2017) Miranda R., Muñoz D. J., Lai D., 2017, MNRAS, 466, 1170
  * Monin et al. (2007) Monin J. L., Clarke C. J., Prato L., McCabe C., 2007, in Reipurth B., Jewitt D., Keil K., eds, Protostars and Planets V. p. 395 (arXiv:astro-ph/0604031)
  * Muñoz & Lai (2016) Muñoz D. J., Lai D., 2016, ApJ, 827, 43
  * Owen & Kollmeier (2017) Owen J. E., Kollmeier J. A., 2017, MNRAS, 467, 3379
  * Pinte et al. (2016) Pinte C., Dent W. R. F., Ménard F., Hales A., Hill T., Cortes P., de Gregorio-Monsalvo I., 2016, ApJ, 816, 25
  * Price et al. (2018) Price D. J., et al., 2018, MNRAS, 477, 1270
  * Quarles et al. (2018) Quarles B., Satyal S., Kostov V., Kaib N., Haghighipour N., 2018, ApJ, 856, 150
  * Ragusa et al. (2017) Ragusa E., Dipierro G., Lodato G., Laibe G., Price D. J., 2017, MNRAS, 464, 1449
  * Ragusa et al. (2020) Ragusa E., Alexander R., Calcino J., Hirsh K., Price D. J., 2020, MNRAS, 499, 3362
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 500, 33
  * Shi et al. (2012) Shi J.-M., Krolik J. H., Lubow S. H., Hawley J. F., 2012, ApJ, 749, 118
  * Simon & Guilloteau (1992) Simon M., Guilloteau S., 1992, ApJ, 397, L47
  * Smallwood et al. (2019) Smallwood J. L., Lubow S. H., Franchini A., Martin R. G., 2019, MNRAS, 486, 2919
  * Smallwood et al. (2022) Smallwood J. L., Lubow S. H., Martin R. G., 2022, MNRAS, 514, 1249
  * Stone & Gardiner (2009) Stone J. M., Gardiner T., 2009, New Astron., 14, 139
  * Stone et al. (2020) Stone J. M., Tomida K., White C. J., Felker K. G., 2020, ApJS, 249, 4
  * Teague et al. (2018) Teague R., Bae J., Bergin E. A., Birnstiel T., Foreman-Mackey D., 2018, ApJ, 860, L12
  * Teague et al. (2019) Teague R., Bae J., Bergin E. A., 2019, Nature, 574, 378
  * Verrier & Evans (2009) Verrier P. E., Evans N. W., 2009, MNRAS, 394, 1721
  * Villenave et al. (2022) Villenave M., et al., 2022, ApJ, 930, 11
  * Weidenschilling (1977) Weidenschilling S. J., 1977, MNRAS, 180, 57
  * Welsh et al. (2014) Welsh W. F., Orosz J. A., Carter J. A., Fabrycky D. C., 2014, in Haghighipour N., ed., 1 Vol. 293, Formation, Detection, and Characterization of Extrasolar Habitable Planets. pp 125–132 (arXiv:1308.6328), doi:10.1017/S1743921313012684
  * Zanazzi & Lai (2018) Zanazzi J. J., Lai D., 2018, MNRAS, 473, 603
  * Zhu (2019) Zhu Z., 2019, MNRAS, 483, 4221
  * Zhu & Stone (2014) Zhu Z., Stone J. M., 2014, ApJ, 795, 53
  * Zhu et al. (2014) Zhu Z., Stone J. M., Rafikov R. R., Bai X.-n., 2014, ApJ, 785, 122
  * Zhu et al. (2015) Zhu Z., Dong R., Stone J. M., Rafikov R. R., 2015, ApJ, 813, 88
  * Zúñiga-Fernández et al. (2021) Zúñiga-Fernández S., et al., 2021, A&A, 655, A15
  * de Val-Borro et al. (2007) de Val-Borro M., Artymowicz P., D’Angelo G., Peplinski A., 2007, A&A, 471, 1043
  * van der Marel et al. (2013) van der Marel N., et al., 2013, Science, 340, 1199
  * von Zeipel (1910) von Zeipel H., 1910, Astronomische Nachrichten, 183, 345

##  Appendix A Coordinate Transform of Disk Data

Simulations of disks are usually oriented so that their local coordinate axes
are aligned with the simulation axes. This alignment simplifies the analysis
of important disk parameters such as azimuthal slicing along the midplane, or
calculation of vorticity along the disk’s vertical axis.

Disks that are inclined, warped, or precess over time are not always aligned
with the simulation axes, and so these quantities are harder to acquire
without some prior manipulation. We show here our method of manipulating the
simulation data of an inclined disk so that it may be analyzed in a similar
way as a flat, zero-inclination disk. As in the paper, we use spherical-polar
coordinates for this section, though the basic process can be used for any
coordinate system.

Consider a disk that is inclined to the simulation axes at an angle $i$. Any
point in 3D space can be described with _simulation coordinates_
$(r,\theta,\phi)$ or _disk coordinates_ $(r,\theta^{\prime},\phi^{\prime})$,
which are not necessarily the same. The transformation between these two
coordinate systems is given by the equations

$\sin^{2}{\theta^{\prime}}=\left(\sin{\theta}\cos{\phi}\right)^{2}+\left(\sin{\theta}\sin{\phi}\sin{i}+\cos{\theta}\sin{i}\right)^{2}$
(9)

and

$\tan{\phi^{\prime}}=\frac{\sin{\theta}\sin{\phi}\cos{i}+\cos{\theta}\sin{i}}{\sin{\theta}\cos{\phi}}.$
(10)

Equation 9 is equivalent to Equation 44 of Zhu (2019). Since this equation
only gives $\sin^{2}{\theta^{\prime}}$, a degeneracy exists when trying to map
coordinates on the top and bottom halves of the sphere. Therefore, we
calculate $\cos{\theta^{\prime}}=\sqrt{1-\sin^{2}{\theta^{\prime}}}$ and
multiply the value by $-1$ if the original polar angle $\theta$ lies below the
equator of the mapped coordinates. This corresponds to

$\theta=\frac{\pi}{2}-\tan^{-1}{\left(\tan{i}\sin{\phi}\right)}.$ (11)

To construct an untilted disk, we start by creating an evenly spaced grid
along the disk coordinates $\theta^{\prime}$ and $\phi^{\prime}$. For each
grid cell, we transform the grid coordinates to simulation coordinates, locate
the nearest data cell in the original simulation coordinates, and copy that
cell data from the simulation into the grid cell. The result is a data cube
oriented along the disk axes, which can be analyzed using traditional methods.

For velocity, additional transformations must be done to account for the
rotation of the velocity components. The azimuthal distance element in the
disk coordinates $d\phi^{\prime}$ can be expressed in terms of the simulation
coordinates as

$d\phi^{\prime 2}=(rd\theta)^{2}+(r\sin{\theta}d\phi)^{2}.$ (12)

This allows us to calculate the angle between the disk and simulation grids
$\psi$ with

$\tan{\psi}=\frac{d\theta}{\sin{\theta}d\phi}.$ (13)

In practice, we calculate $\psi$ by comparing adjacent cells in the disk
coordinates and converting back to the simulation coordinates, allowing the
distances to be calculated using simple differences. Once the local rotation
angle is found, the velocity in the disk coordinates
$\mathrm{v}_{\theta}^{\prime}$ and $\mathrm{v}_{\phi}^{\prime}$ can be
calculated with

$\mathrm{v}_{\theta}^{\prime}=\mathrm{v}_{\theta}\cos{\psi}-\mathrm{v}_{\phi}\sin{\psi}$
(14)

and

$\mathrm{v}_{\phi}^{\prime}=\mathrm{v}_{\theta}\sin{\psi}+\mathrm{v}_{\phi}\cos{\psi}.$
(15)

The radial velocity remains unchanged during the rotation, i.e.
$\mathrm{v}_{r}^{\prime}=\mathrm{v}_{r}$.

Figure 11 shows an example of the transformation applied to a 2D spherical
slice of our simulation data at $t=0$. The transformed density and velocity
are similar to that of a disk aligned with the coordinate system. Small
artifacts are visible in the areas close to the poles of the simulation
coordinates due to the lack of initial data in these regions. Small vertical
striations can also be seen along the transformed data,

For simulations in this paper, we transform the disk data using the total
angular momentum vector $\overset{\rightharpoonup}{\mathbf{L}}$ to calculate
$i$ and $\mathbf{\Omega}$ using Equations 4 and 5. Warped or broken disks can
be analyzed in a similar fashion by using the local angular momentum vector
$\overset{\rightharpoonup}{\mathbf{L}}(R)$ to transform each annulus
separately.

Figure 11: Example of the coordinate transform using a 2D slice of the
simulation data at $t=0$. _Top Row:_ Coordinate grid lines for the simulation
and disk coordinates. _Middle Row:_ Disk density in the simulation and disk
coordinates. _Bottom Row:_ Azimuthal velocity $\mathrm{v}_{\phi}$ in the
simulation and disk coordinates.


Combining inputs from multiple sensor modalities effectively in reinforcement learning (RL) is an open problem.
While many self-supervised representation learning approaches exist to improve performance and sample complexity for image-based RL, they usually neglect other available information, such as robot proprioception.
However, using this proprioception for representation learning can help algorithms to focus on relevant aspects and guide them toward finding better representations.
In this work, we systematically analyze representation learning for RL from multiple sensors by building on Recurrent State Space Models.
We propose a combination of reconstruction-based and contrastive losses, which allows us to choose the most appropriate method for each sensor modality.
We demonstrate the benefits of joint representations, particularly with distinct loss functions for each modality, for model-free and model-based RL on complex tasks.
Those include tasks where the images contain distractions or occlusions and a new locomotion suite.
We show that combining reconstruction-based and contrastive losses for joint representation learning improves performance significantly compared to a post hoc combination of image representations and proprioception and can also improve the quality of learned models for model-based RL.

[Agarwal et al., 2021]
Rishabh Agarwal, Max Schwarzer, Pablo Samuel Castro, Aaron C Courville, and
Marc Bellemare.
Deep reinforcement learning at the edge of the statistical precipice.
Advances in neural information processing systems,
34:0 29304–29320, 2021.

[Ba et al., 2016]
Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E Hinton.
Layer normalization.
arXiv preprint arXiv:1607.06450, 2016.

[Banijamali et al., 2018]
Ershad Banijamali, Rui Shu, Hung Bui, Ali Ghodsi, et al.
Robust locally-linear controllable embedding.
In International Conference on Artificial Intelligence and
Statistics, pages 1751–1759. PMLR, 2018.

[Becker and Neumann, 2022]
Philipp Becker and Gerhard Neumann.
On uncertainty in deep state space models for model-based
reinforcement learning.
Transactions on Machine Learning Research, 2022.

[Brohan et al., 2022]
Anthony Brohan, Yevgen Chebotar, Chelsea Finn, Karol Hausman, Alexander Herzog,
Daniel Ho, Julian Ibarz, Alex Irpan, Eric Jang, Ryan Julian, et al.
Do as i can, not as i say: Grounding language in robotic affordances.
In 6th Annual Conference on Robot Learning, 2022.

[Cho et al., 2014]
Kyunghyun Cho, Bart Van Merriënboer, Dzmitry Bahdanau, and Yoshua Bengio.
On the properties of neural machine translation: Encoder-decoder
arXiv preprint arXiv:1409.1259, 2014.

[Deng et al., 2022]
Fei Deng, Ingook Jang, and Sungjin Ahn.
Dreamerpro: Reconstruction-free model-based reinforcement learning
with prototypical representations.
In International Conference on Machine Learning, pages
4956–4975. PMLR, 2022.

[Driess et al., 2023]
Danny Driess, Fei Xia, Mehdi SM Sajjadi, Corey Lynch, Aakanksha Chowdhery,
Brian Ichter, Ayzaan Wahid, Jonathan Tompson, Quan Vuong, Tianhe Yu, et al.
Palm-e: An embodied multimodal language model.
arXiv preprint arXiv:2303.03378, 2023.

[Finn et al., 2016]
Chelsea Finn, Ian Goodfellow, and Sergey Levine.
Unsupervised learning for physical interaction through video
Advances in neural information processing systems, 29, 2016.

[Fu et al., 2021]
Xiang Fu, Ge Yang, Pulkit Agrawal, and Tommi Jaakkola.
Learning task informed abstractions.
In International Conference on Machine Learning, pages
3480–3491. PMLR, 2021.

[Fu et al., 2022]
Zipeng Fu, Ashish Kumar, Ananye Agarwal, Haozhi Qi, Jitendra Malik, and Deepak
Coupling vision and proprioception for navigation of legged robots.
In Proceedings of the IEEE/CVF Conference on Computer Vision
and Pattern Recognition (CVPR), pages 17273–17283, June 2022.

[Ha and Schmidhuber, 2018]
David Ha and Jürgen Schmidhuber.
World models.
arXiv preprint arXiv:1803.10122, 2018.

[Haarnoja et al., 2018]
Tuomas Haarnoja, Aurick Zhou, Kristian Hartikainen, George Tucker, Sehoon Ha,
Jie Tan, Vikash Kumar, Henry Zhu, Abhishek Gupta, Pieter Abbeel, et al.
Soft actor-critic algorithms and applications.
arXiv preprint arXiv:1812.05905, 2018.

[Hafner et al., 2019]
Danijar Hafner, Timothy Lillicrap, Ian Fischer, Ruben Villegas, David Ha,
Honglak Lee, and James Davidson.
Learning latent dynamics for planning from pixels.
In International conference on machine learning, pages
2555–2565. PMLR, 2019.

[Hafner et al., 2020]
Danijar Hafner, Timothy Lillicrap, Jimmy Ba, and Mohammad Norouzi.
Dream to control: Learning behaviors by latent imagination.
In International Conference on Learning Representations, 2020.

[Hafner et al., 2021]
Danijar Hafner, Timothy P Lillicrap, Mohammad Norouzi, and Jimmy Ba.
Mastering atari with discrete world models.
In International Conference on Learning Representations, 2021.

[Hafner et al., 2022]
Danijar Hafner, Kuang-Huei Lee, Ian Fischer, and Pieter Abbeel.
Deep hierarchical planning from pixels.
In Advances in Neural Information Processing Systems, 2022.

[Hafner et al., 2023]
Danijar Hafner, Jurgis Pasukonis, Jimmy Ba, and Timothy Lillicrap.
Mastering diverse domains through world models.
arXiv preprint arXiv:2301.04104, 2023.

[Kalashnikov et al., 2018]
Dmitry Kalashnikov, Alex Irpan, Peter Pastor, Julian Ibarz, Alexander Herzog,
Eric Jang, Deirdre Quillen, Ethan Holly, Mrinal Kalakrishnan, Vincent
Vanhoucke, et al.
Scalable deep reinforcement learning for vision-based robotic
In Conference on Robot Learning, pages 651–673. PMLR, 2018.

[Kay et al., 2017]
Will Kay, Joao Carreira, Karen Simonyan, Brian Zhang, Chloe Hillier, Sudheendra
Vijayanarasimhan, Fabio Viola, Tim Green, Trevor Back, Paul Natsev, et al.
The kinetics human action video dataset.
arXiv preprint arXiv:1705.06950, 2017.

[Kingma and Ba, 2015]
Diederik P Kingma and Jimmy Ba.
Adam: A method for stochastic optimization.
In ICLR, 2015.

[Kingma and Welling, 2013]
Diederik P Kingma and Max Welling.
Auto-encoding variational bayes.
arXiv preprint arXiv:1312.6114, 2013.

[Lee et al., 2020]
Alex X Lee, Anusha Nagabandi, Pieter Abbeel, and Sergey Levine.
Stochastic latent actor-critic: Deep reinforcement learning with a
latent variable model.
Advances in Neural Information Processing Systems,
33:0 741–752, 2020.

[Levine et al., 2016]
Sergey Levine, Chelsea Finn, Trevor Darrell, and Pieter Abbeel.
End-to-end training of deep visuomotor policies.
The Journal of Machine Learning Research, 170
(1):0 1334–1373, 2016.

[Ma et al., 2020]
Xiao Ma, Siwei Chen, David Hsu, and Wee Sun Lee.
Contrastive variational reinforcement learning for complex
In Conference on Robot Learning, pages 959–972. PMLR, 2020.

[Majumdar et al., 2023]
Arjun Majumdar, Karmesh Yadav, Sergio Arnaud, Yecheng Jason Ma, Claire Chen,
Sneha Silwal, Aryan Jain, Vincent-Pierre Berges, Pieter Abbeel, Jitendra
Malik, et al.
Where are we in the search for an artificial visual cortex for
embodied intelligence?
arXiv preprint arXiv:2303.18240, 2023.

[Mees et al., 2022]
Oier Mees, Lukas Hermann, and Wolfram Burgard.
What matters in language conditioned robotic imitation learning over
unstructured data.
IEEE Robotics and Automation Letters, 70 (4):0
11205–11212, 2022.

[Murphy, 2012]
Kevin P Murphy.
Machine learning: a probabilistic perspective.
MIT press, 2012.

[Nair et al., 2022]
Suraj Nair, Aravind Rajeswaran, Vikash Kumar, Chelsea Finn, and Abhinav Gupta.
R3m: A universal visual representation for robot manipulation.
In 6th Annual Conference on Robot Learning, 2022.

[Nguyen et al., 2021]
Tung D Nguyen, Rui Shu, Tuan Pham, Hung Bui, and Stefano Ermon.
Temporal predictive coding for model-based planning in latent space.
In International Conference on Machine Learning, pages
8130–8139. PMLR, 2021.

[Okada and Taniguchi, 2021]
Masashi Okada and Tadahiro Taniguchi.
Dreaming: Model-based reinforcement learning by latent imagination
without reconstruction.
In 2021 IEEE International Conference on Robotics and
Automation (ICRA), pages 4209–4215. IEEE, 2021.

[Okada and Taniguchi, 2022]
Masashi Okada and Tadahiro Taniguchi.
Dreamingv2: Reinforcement learning with discrete world models without
In 2022 IEEE/RSJ International Conference on Intelligent Robots
and Systems (IROS), pages 985–991. IEEE, 2022.

[Oord et al., 2018]
Aaron van den Oord, Yazhe Li, and Oriol Vinyals.
Representation learning with contrastive predictive coding.
arXiv preprint arXiv:1807.03748, 2018.

[Parisi et al., 2022]
Simone Parisi, Aravind Rajeswaran, Senthil Purushwalkam, and Abhinav Gupta.
The unsurprising effectiveness of pre-trained vision models for
In International Conference on Machine Learning, pages
17359–17371. PMLR, 2022.

[Poole et al., 2019]
Ben Poole, Sherjil Ozair, Aaron Van Den Oord, Alex Alemi, and George Tucker.
On variational bounds of mutual information.
In International Conference on Machine Learning, pages
5171–5180. PMLR, 2019.

[Radford et al., 2021]
Alec Radford, Jong Wook Kim, Chris Hallacy, Aditya Ramesh, Gabriel Goh,
Sandhini Agarwal, Girish Sastry, Amanda Askell, Pamela Mishkin, Jack Clark,
et al.
Learning transferable visual models from natural language
In International Conference on Machine Learning, pages
8748–8763. PMLR, 2021.

[Ramesh et al., 2022]
Aditya Ramesh, Prafulla Dhariwal, Alex Nichol, Casey Chu, and Mark Chen.
Hierarchical text-conditional image generation with clip latents.
arXiv preprint arXiv:2204.06125, 2022.

[Rezende et al., 2014]
Danilo Jimenez Rezende, Shakir Mohamed, and Daan Wierstra.
Stochastic backpropagation and approximate inference in deep
generative models.
In International conference on machine learning, pages
1278–1286. PMLR, 2014.

[Seo et al., 2022]
Younggyo Seo, Danijar Hafner, Hao Liu, Fangchen Liu, Stephen James, Kimin Lee,
and Pieter Abbeel.
Masked world models for visual control.
In 6th Annual Conference on Robot Learning, 2022.

[Shridhar et al., 2022]
Mohit Shridhar, Lucas Manuelli, and Dieter Fox.
Cliport: What and where pathways for robotic manipulation.
In Conference on Robot Learning, pages 894–906. PMLR, 2022.

[Shridhar et al., 2023]
Mohit Shridhar, Lucas Manuelli, and Dieter Fox.
Perceiver-actor: A multi-task transformer for robotic manipulation.
In Conference on Robot Learning, pages 785–799. PMLR, 2023.

[Shu et al., 2020]
Rui Shu, Tung Nguyen, Yinlam Chow, Tuan Pham, Khoat Than, Mohammad Ghavamzadeh,
Stefano Ermon, and Hung Bui.
Predictive coding for locally-linear control.
In International Conference on Machine Learning, pages
8862–8871. PMLR, 2020.

[Sohn et al., 2015]
Kihyuk Sohn, Honglak Lee, and Xinchen Yan.
Learning structured output representation using deep conditional
generative models.
Advances in neural information processing systems,
28:0 3483–3491, 2015.

[Srinivas et al., 2020]
Aravind Srinivas, Michael Laskin, and Pieter Abbeel.
Curl: Contrastive unsupervised representations for reinforcement
arXiv preprint arXiv:2004.04136, 2020.

[Srivastava et al., 2021]
Nitish Srivastava, Walter Talbott, Martin Bertran Lopez, Shuangfei Zhai, and
Joshua M. Susskind.
Robust robotic control from pixels using contrastive recurrent
state-space models.
In Deep RL Workshop NeurIPS 2021, 2021.

[Stooke et al., 2021]
Adam Stooke, Kimin Lee, Pieter Abbeel, and Michael Laskin.
Decoupling representation learning from reinforcement learning.
In International Conference on Machine Learning, pages
9870–9879. PMLR, 2021.

[Tassa et al., 2018]
Yuval Tassa, Yotam Doron, Alistair Muldal, Tom Erez, Yazhe Li, Diego de Las
Casas, David Budden, Abbas Abdolmaleki, Josh Merel, Andrew Lefrancq, et al.
Deepmind control suite.
arXiv preprint arXiv:1801.00690, 2018.

[Tassa et al., 2020]
Yuval Tassa, Saran Tunyasuvunakool, Alistair Muldal, Yotam Doron, Siqi Liu,
Steven Bohez, Josh Merel, Tom Erez, Timothy Lillicrap, and Nicolas Heess.
dm_control: Software and tasks for continuous control.
Software Impacts, 6:0 100022, 2020.

[Tomar et al., 2023]
Manan Tomar, Utkarsh Aashu Mishra, Amy Zhang, and Matthew E Taylor.
Learning representations for pixel-based control: What matters and
Transactions on Machine Learning Research, 2023.

[Wahlström et al., 2015]
Niklas Wahlström, Thomas B Schön, and Marc Peter Deisenroth.
From pixels to torques: Policy learning with deep dynamical models.
arXiv preprint arXiv:1502.02251, 2015.

[Wang et al., 2022]
Tongzhou Wang, Simon Du, Antonio Torralba, Phillip Isola, Amy Zhang, and
Yuandong Tian.
Denoised mdps: Learning world models better than the world itself.
In International Conference on Machine Learning, pages
22591–22612. PMLR, 2022.

[Watter et al., 2015]
Manuel Watter, Jost Springenberg, Joschka Boedecker, and Martin Riedmiller.
Embed to control: A locally linear latent dynamics model for control
from raw images.
In Advances in neural information processing systems, pages
2746–2754, 2015.

[Wu et al., 2022]
Philipp Wu, Alejandro Escontrela, Danijar Hafner, Pieter Abbeel, and Ken
Daydreamer: World models for physical robot learning.
In 6th Annual Conference on Robot Learning, 2022.

[Xiao et al., 2022]
Tete Xiao, Ilija Radosavovic, Trevor Darrell, and Jitendra Malik.
Masked visual pre-training for motor control.
arXiv preprint arXiv:2203.06173, 2022.

[Yarats et al., 2021]
Denis Yarats, Rob Fergus, Alessandro Lazaric, and Lerrel Pinto.
Reinforcement learning with prototypical representations.
In International Conference on Machine Learning, pages
11920–11931. PMLR, 2021a.

[Yarats et al., 2021]
Denis Yarats, Amy Zhang, Ilya Kostrikov, Brandon Amos, Joelle Pineau, and Rob
Improving sample efficiency in model-free reinforcement learning from
In Proceedings of the AAAI Conference on Artificial
Intelligence, volume 35, pages 10674–10681, 2021b.

[Yarats et al., 2022]
Denis Yarats, Rob Fergus, Alessandro Lazaric, and Lerrel Pinto.
Mastering visual continuous control: Improved data-augmented
reinforcement learning.
In International Conference on Learning Representations, 2022.

[You et al., 2022]
Bang You, Oleg Arenz, Youping Chen, and Jan Peters.
Integrating contrastive learning with dynamic models for
reinforcement learning from images.
Neurocomputing, 476:0 102–114, 2022.

[Zhang et al., 2020]
Amy Zhang, Rowan Thomas McAllister, Roberto Calandra, Yarin Gal, and Sergey
Learning invariant representations for reinforcement learning without
In International Conference on Learning Representations, 2020.