# Badminton Birdie-Like Aerodynamic Alignment of Drifting Dust Grains by
Subsonic Gaseous Flows in Protoplanetary Disks

Zhe-Yu Daniel Lin,1,2 Zhi-Yun Li,1 Haifeng Yang,3 Leslie W. Looney,4,5 Ian W.
Stephens,6 Manuel Fernández-López,7 Rachel E. Harrison8
1Department of Astronomy, University of Virginia, 530 McCormick Rd.,
Charlottesville 22904, Virginia, USA
2Earth and Planets Laboratory, Carnegie Institution for Science, 5241 Broad
Branch Rd. NW, Washington, DC 20015, USA
3Kavli Institute for Astronomy and Astrophysics, Peking University, Yi He Yuan
Lu 5, Haidian Qu, Beijing 100871, People’s Republic of China
4Department of Astronomy, University of Illinois, 1002 W Green St., Urbana,
Illinois 61801, USA
5National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville,
Virginia 22903, USA
6Department of Earth, Environment, and Physics, Worcester State University,
Worcester, Massachusetts 01602, USA
7Instituto Argentino de Radioastronomía, CCT-La Plata (CONICET), C.C.5, 1894,
Villa Elisa, Argentina
8School of Physics and Astronomy, Monash University, Clayton VIC 3800,
Australia
E-mail<EMAIL_ADDRESS>or<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Recent (sub)millimeter polarization observations of protoplanetary disks
reveal toroidally aligned, effectively prolate dust grains large enough (at
least $\sim 100$ $\mu$m) to efficiently scatter millimeter light. The
alignment mechanism for these grains remains unclear. We explore the
possibility that gas drag aligns grains through gas-dust relative motion when
the grain’s center of mass is offset from its geometric center, analogous to a
badminton birdie’s alignment in flight. A simple grain model of two non-
identical spheres illustrates how a grain undergoes damped oscillations from
flow-induced restoring torques which align its geometric center in the flow
direction relative to its center of mass. Assuming specular reflection and
subsonic flow, we derive an analytical equation of motion for spheroids where
the center of mass can be shifted away from the spheroid’s geometric center.
We show that a prolate or an oblate grain can be aligned with the long axis
parallel to the gas flow when the center of mass is shifted along that axis.
Both scenarios can explain the required effectively prolate grains inferred
from observations. Application to a simple disk model shows that the alignment
timescales are shorter than or comparable to the orbital time. The grain
alignment direction in a disk depends on the disk (sub-)structure and grain
Stokes number (St) with azimuthal alignment for large St grains in sub-
Keplerian smooth gas disks and for small St grains near the gas pressure
extrema, such as rings and gaps.

###### keywords:

polarization – protoplanetary disks

††pubyear: 2015††pagerange: Badminton Birdie-Like Aerodynamic Alignment of
Drifting Dust Grains by Subsonic Gaseous Flows in Protoplanetary Disks–C

##  1 Introduction

Dust grains within gaseous protoplanetary disks are the raw material for
forming planets (e.g. Drążkowska et al., 2023). Their studies have been
revolutionized by the Atacama Large Millimeter/submillimeter Array (ALMA),
especially through disk-scale observations of linear continuum polarization
(e.g. Kataoka et al., 2016b; Stephens et al., 2017; Lee et al., 2018; Girart
et al., 2018; Ohashi et al., 2018; Alves et al., 2018; Sadavoy et al., 2018;
Bacciotti et al., 2018; Dent et al., 2019; Takahashi et al., 2019; Harrison et
al., 2019; Sadavoy et al., 2019; Stephens et al., 2020; Aso et al., 2021; Tang
et al., 2023; Stephens et al., 2023; Lin et al., 2024a; Liu et al., 2024;
Harrison et al., 2024). The origin of disk polarization is often attributed to
scattering when dust grains have grown large enough to efficiently scatter
millimeter/sub-millimeter light. Measuring polarization can constrain the
properties of dust grains, like the distribution, grain sizes, porosity, etc
(e.g. Kataoka et al., 2015; Yang et al., 2016a; Kataoka et al., 2016a; Yang et
al., 2017; Kirchschlager & Bertrang, 2020; Lin et al., 2020b; Yang & Li, 2020;
Zhang et al., 2023; Zamponi et al., 2024; Yang et al., 2024). Alternatively,
if grains are aligned, then the thermal emission is intrinsically polarized
(e.g. van de Hulst, 1957; Yang et al., 2016b; Kirchschlager et al., 2019; Lin
et al., 2020a, 2022).

Polarized thermal emission from aligned grains has been used to trace the
magnetic field in the envelopes around protostars. Polarization successfully
measured the expected hour-glass morphology of the magnetic field (e.g. Girart
et al., 2006; Stephens et al., 2013; Maury et al., 2018; Kwon et al., 2019;
Huang et al., 2024). The widely accepted explanation is through radiative
alignment torques, i.e., RATs (Dolginov, 1972; Dolginov & Mytrophanov, 1976;
Draine & Weingartner, 1996, 1997; Lazarian & Hoang, 2007a) which allow grains
to be aligned to the magnetic field or radiation field through internal and
external alignment (see, e.g., Lazarian 2007, Lazarian et al. 2015 for a
review). Internal alignment refers to the alignment of the angular momentum,
$\bm{J}$, to the axis of maximum moment of inertia due to energy dissipation
(Purcell, 1979; Lazarian & Efroimsky, 1999; Hoang & Lazarian, 2009; Hoang et
al., 2022). As a result, the long axis of the grain should be perpendicular to
$\bm{J}$ and the polarization of an ensemble of grains should appear
effectively oblate. External alignment occurs when $\bm{J}$ becomes aligned
with an external field, like the magnetic field $\bm{B}$ or the radiation
field (Lazarian & Hoang, 2007a; Lazarian et al., 2015). When producing
polarization after internal and external alignments, the grains as an ensemble
should be effectively oblate with their short axes aligned along the external
field.

In the disk-scales, however, recent observations have revealed compelling
evidence for the presence of grains that appear effectively prolate with their
long axes aligned toroidally around the disk. For example, multiwavelength
observations of HL Tau have unveiled transitions in the polarization
morphology from the Very Large Array (VLA) Q-band ($\lambda=7.1$ mm) to ALMA
Bands 3, 4, 5, 6, and 7 ranging from $\lambda=$3.1 mm to 0.87 mm (Stephens et
al., 2014, 2017; Lin et al., 2024a). The change in pattern cannot be solely
attributed to scattering. In particular, Band 3 ($\lambda=$3.1 mm) shows an
azimuthally oriented polarization pattern that was shown to require toroidally
aligned, effectively prolate grains based on the azimuthal variation in the
polarization (Kataoka et al., 2017; Mori & Kataoka, 2021; Yang et al., 2019).
Since the optical depth decreases with increasing wavelength, scattering can
become less of a factor at longer wavelengths. Multiwavelength models have
shown that scattering of toroidally aligned, prolate grains can produce the
observed polarization transition (Lin et al., 2022, 2024a). In addition,
recent high angular resolution observations ($5$ au resolution) at ALMA Band 7
($\lambda=0.87$ mm) revealed azimuthally oriented polarization in the first
gap, but signatures of scattering polarization in the rings. Stephens et al.
(2023) demonstrated that by incorporating optical depth changes between rings
and gaps, scattering of toroidally aligned, prolate grains can also explain
the observed polarization substructure.

Signatures of toroidally aligned, prolate grains akin to those observed from
HL Tau also exist in disks around other sources. Haro 6-13, V892 Tau, DG Tau,
and GG Tau show polarization that is predominantly parallel to the disk minor
axis (a signature of scattering) at the shorter wavelength, $\lambda\sim 0.9$
mm, and azimuthally oriented polarization at the longer wavelength,
$\lambda\sim 3$ mm (Bacciotti et al., 2018; Harrison et al., 2019; Tang et
al., 2023; Ohashi et al., 2023; Harrison et al., 2024). For AS 209, the
polarization along the disk major axis has polarization angles parallel to the
disk minor axis, but the polarization in the outer regions appears azimuthally
oriented (at Band 7 from Mori et al. 2019 and at Band 6 from Harrison et al.
2021). The polarization pattern is reminiscent of HL Tau where the optical
depth in the inner region is large while the optical depth in the outer region
is low leading to a scattering signature in the inner region and a signature
of toroidally aligned, prolate grains in the outer region (Stephens et al.,
2017; Lin et al., 2024a).

An alternative to radiative alignment is mechanical alignment where
interactions with the gas can align the dust grains. One type of mechanical
alignment is the Gold mechanism that allows the alignment of grains along
their long axes under the presence of gas-dust relative motion (Gold, 1952a,
b; Purcell, 1969). Although the mechanism can explain the necessary prolate
grains, the required relative motion should be supersonic which does not apply
to protoplanetary disks (e.g. Nakagawa et al., 1986; Lesur et al., 2023).

Another type of mechanical alignment is through the helicity of grains, where
grains spin up under the presence of gas-dust relative motion (Lazarian &
Hoang, 2007b). This mechanism also requires internal alignment, which, along
with external alignment, ensures that the short axes of grains become aligned
to the direction of the gas flow or the magnetic field $\bm{B}$. While the
proposed mechanism no longer requires supersonic drift, it is expected to
produce effectively oblate grains because of internal alignment. The
observation of prolate grains suggests a lack of internal alignment which may
not be too surprising since internal alignment timescales are much longer than
the gas damping timescale for grains larger than $\sim 10$ $\mu$m in
protoplanetary disk environments (Hoang et al., 2022).

A potential solution is through a modification to the radiative alignment
paradigm using the so-called “wrong” internal alignment (Hoang & Lazarian,
2009). For large grains without internal alignment, the short axes can become
perpendicular to $\bm{J}$ (in contrast to the “right” alignment mentioned
above where the short axes are parallel to $\bm{J}$). When $\bm{J}$ becomes
externally aligned to the expected toroidal $\bm{B}$-field of a protoplanetary
disk, it can match the desired toroidally aligned prolate grains if grains
contain large amounts of iron clusters (Nguyen Tat et al., 2023). However,
alignment with the $\bm{B}$-field appears difficult in disk environments
because the Larmor precession timescale can be much longer than the gas
damping timescale unless the former is greatly shortened by, e.g.,
superparamagnetic inclusions (Tazaki et al., 2017; Yang, 2021).

The various examples from observations call for an answer to why grains are
toroidally aligned and effectively prolate, but a robust explanation remains
inconclusive. Understanding the grain alignment mechanism is not only
interesting in its own right, but it will also permit the measurement of the
underlying grain-aligning vector field in the disk. In this paper, we
investigate a novel type of mechanical alignment, where the offset between the
center of mass of a grain and its geometric center creates restoring torques
when the grain feels a systematic flow of gas with respect to its center of
mass, i.e., an aerodynamic flow (denoted by the vector “$\bm{A}$” hereafter).
The mechanism is motivated by how a badminton birdie flying through the air
oscillates and quickly aligns aerodynamically along the direction of flight
against the gaseous headwind in terrestrial environments 111Other examples
from human activities include the jiànzi or the hago (in the game of
hanetsuki).  (Cohen et al., 2015). To highlight the central role of the offset
between geometric and mass centers, we will term the mechanism “badminton
birdie-like alignment" (or “birdie-like alignment" for short) to distinguish
it from other forms of mechanical alignment (the expected degree of the offset
will be discussed below). We ignore internal relaxation since the effect is
slow for large grains and only consider the rotational motion of grains due to
torques provided by gas drag. In disks with subsonic gas-dust bulk relative
motion, we will show that the proposed mechanism can achieve fast alignment
within an orbital time for most regions in disks, which opens a new window on
mapping the field of aerodynamic flow acting on grains — the $\bm{A}$-field —
through dust continuum polarization observations.

The paper is organized as follows. Section 2 builds the physical intuition by
using a simple grain model composed of two spheres to illustrate the key
alignment behavior and alignment timescales under the presence of a gas flow.
Section 3 derives the torque on a grain of arbitrary shape and calculates the
alignment of spheroids (both prolate and oblate) whose solutions for
polarization are known enabling connections to observed disk polarization from
aligned grains. Section 4 implements the spheroid alignment model in a simple
disk model with gas and dust velocity fields to assess the alignment timescale
and alignment directions. Section 5 provides a discussion, and we summarize
our results in Section 6.

##  2 Double-Sphere Illustration

This section begins with the simplest model to illustrate the alignment under
gas drag using a double-sphere model (Cohen et al., 2015). Given the simple
spherical shape, the gas drag in the Epstein regime has a simple analytical
form (Epstein, 1924). We first describe the necessary formulation of the
problem, and then parameterize the model to focus on the physical quantities
that determine alignment. We will show that the asymmetry between the two
spheres induces oscillation due to a systematic gas-dust relative motion. The
spin of the grain through gas leads to damping of the oscillation and
eventually to alignment.

###  2.1 Problem Setup

Figure 1:  Schematic of the double sphere model. The dust grain is composed of
two spheres, labeled “1” and “2,” with radii $\varsigma_{1}$ and
$\varsigma_{2}$, respectively. Point $O$ is the center of mass of the grain.
The radius vectors from the center of mass are $\bm{r}_{1}$ and $\bm{r}_{2}$,
respectively. $\bm{A}$ is the gas velocity relative to the center of mass of
the grain, i.e., the aerodynamic flow from the perspective of the grain, and
is along $\bm{e}_{3}$. The spin of the grain is along $\bm{e}_{2}$. $\theta$
is the angle of $\bm{r}_{2}$ from $\bm{e}_{3}$ in the counter-clockwise
direction in this figure. Since alignment occurs at $\theta=0$ when
$m_{1}/m_{2}>\varsigma_{1}^{2}/\varsigma_{2}^{2}$, sphere 1 is colored darker
to denote a larger mass than sphere 2 and is drawn smaller than sphere 2.

Consider a grain composed of two spheres, sphere 1 and sphere 2, connected by
a massless rigid pole. We use point $O$ to denote the center of mass and let
$\bm{A}$ be the systematic velocity of the gas with respect to $O$. Fig. 1
shows a schematic of the double-sphere. The radius vectors from $O$ to the
center of each spheres are $\bm{r}_{1}$ and $\bm{r}_{2}$, respectively, with
corresponding lengths denoted as $r_{1}$ and $r_{2}$. By definition of the
center of mass, $r_{i}$ is related to the mass through

$\displaystyle m_{1}r_{1}=m_{2}r_{2}\text{ ,}$ (1)

where $m_{1}$ and $m_{2}$ are the masses of the spheres. The $i$th sphere is
characterized by radius $\varsigma_{i}$ and its radius and mass are related
through

$\displaystyle m_{i}=\frac{4\pi}{3}\varsigma_{i}^{3}\rho_{i}\text{ ,}$ (2)

where $\rho_{i}$ is the material density for the $i$th sphere, which can be
different for the spheres.

We use the Cartesian coordinate system centered on the center of mass $O$ with
unit directions $\bm{e}_{1}$, $\bm{e}_{2}$, and $\bm{e}_{3}$. Let $\bm{e}_{3}$
be along the direction of the gas flow $\bm{A}$. We define the angle between
$\bm{r}_{2}$ and $\bm{e}_{3}$ as $\theta$. Without loss of generality, let
$\bm{e}_{1}$ be in the plane formed by $\bm{r}_{2}$ and $\bm{e}_{3}$. The
aerodynamic drag force will be in this plane, producing a torque around the
$\bm{e}_{2}$ axis, which will be the axis of the flow-induced grain rotating
motion, denoted by $\bm{\omega}=\omega\bm{e}_{2}$.

The drag force in the Epstein regime is (Epstein, 1924):

$\displaystyle\bm{F}=-\frac{4\pi}{3}\rho_{g}\varsigma^{2}v_{\text{th}}\bm{u}$
(3)

for a spherical grain of radius $\varsigma$ embedded in a gas with mass
density $\rho_{g}$. For convenience, we will also use the gas number density
$n_{g}\equiv\rho_{g}/(\mu m_{p})$ where $m_{p}$ is the mass of a proton and
$\mu$ is the mean molecule weight which we assume as 2.3 for molecular
hydrogen-dominated gas. $v_{\text{th}}$ is the average speed of the gas
molecules with a Maxwell-Boltzmann distribution

$\displaystyle v_{\text{th}}\equiv\sqrt{\frac{8kT}{\pi\mu m_{p}}}$ (4)

where $k$ is the Boltzmann constant and $T$ is the gas temperature. $\bm{u}$
is the velocity of the sphere relative to the gas 222There are several
quantities related to velocity used throughout the paper that deserve
clarification. We use “$A$” to denote the bulk motion of gas relative to the
center of mass of a grain. We use “$u$” to denote the motion of a piece of a
grain relative to the gas. For Sec. 2, $u$ is the speed for one of the
spheres, while for Sec. 3, $u$ corresponds to the speed of an infinitesimal
surface. In Sec. 4, we will use “$v$” to describe the velocity field of the
disk. . The Epstein regime applies if the mean free path of the gas molecules
is much larger than the size of the particle. For disks, the mean free path is
generally of order $\sim$cm near $\sim 1$ au and larger at larger radii where
the gas density is lower (Armitage, 2015). Thus, the Epstein regime is
appropriate in the outer disk where resolved dust polarization is detected.

Each sphere provides a cross-section for gas drag to induce a torque
$\bm{\Gamma}_{i}\equiv\bm{r}_{i}\times\bm{F}_{i}$ on the grain. The velocity
$\bm{u}_{i}$ for the $i$th sphere traveling relative to the gas is

$\displaystyle\bm{u}_{i}\equiv\bm{\omega}\times\bm{r}_{i}-\bm{A}$ (5)

when incorporating the rotation of the grain (with angular velocity
$\bm{\omega}$). Applying the relative velocity from Eq. 5 to the drag force
from Eq. 3 gives

$\displaystyle\bm{\Gamma}_{i}=-\frac{4\pi}{3}\rho_{g}v_{\text{th}}\varsigma_{i}^{2}[\bm{r}_{i}\times(\bm{\omega}\times\bm{r}_{i})-\bm{r}_{i}\times\bm{A}].$
(6)

The total torque on the grain is $\bm{\Gamma}=\bm{\Gamma}_{1}+\bm{\Gamma}_{2}$
from both spheres.

We can easily identify that the torque $\bm{\Gamma}$ is only non-zero along
$\bm{e}_{2}$ meaning we only have to consider the component $\Gamma_{2}$
($\Gamma_{1}=\Gamma_{3}=0$; not to be confused with $\bm{\Gamma}_{2}$ which is
the torque from the second sphere). There is only one equation of motion:

$\displaystyle\Gamma_{2}=I\dot{\omega}$ (7)

where $I$ is the moment of inertia and $\dot{\omega}$ is the angular
acceleration. For the double-sphere model, we can easily derive the moment of
inertia as

$\displaystyle
I=\sum_{i}m_{i}\bigg{(}\frac{2}{5}\varsigma_{i}^{2}+r_{i}^{2}\bigg{)}.$ (8)

where the second term comes from the parallel axis theorem.

Using Eq. 6 and 7, we get a second-order differential equation:

$\displaystyle I\dot{\omega}+D\omega+P\sin\theta$ $\displaystyle=0$ (9)

where

$\displaystyle D$
$\displaystyle\equiv\frac{4\pi}{3}\rho_{g}v_{\text{th}}(\varsigma_{1}^{2}r_{1}^{2}+\varsigma_{2}^{2}r_{2}^{2})$
$\displaystyle P$
$\displaystyle\equiv\frac{4\pi}{3}\rho_{g}v_{\text{th}}A(-\varsigma_{1}^{2}r_{1}+\varsigma_{2}^{2}r_{2}).$

The second term is a torque that is proportional to $\omega$ and acts to
oppose it, serving as the damping torque. The third term is a
$\theta$-dependent torque induced by the gas flow through relative dust-gas
drift.

In the limit of small $\theta$, Eq. 9 simply describes a damped harmonic
oscillator. The motion of a damped harmonic oscillator is characterized by two
timescales: the undamped period of oscillation $t_{o}$ and the damping time
$t_{d}$. The undamped period of oscillation is simply

$\displaystyle t_{o}\equiv\frac{2\pi}{\omega_{o}}$ (10)

where $\omega_{o}$ is the undamped angular frequency

$\displaystyle\omega_{o}\equiv\sqrt{\frac{P}{I}}.$ (11)

The damping time describes the timescale for the oscillation amplitude to
decrease:

$\displaystyle t_{d}\equiv\frac{2I}{D}.$ (12)

Although these expressions are derived from analytical solutions to the damped
harmonic oscillator, it remains beneficial to use them as characteristic
timescales for Eq. 9 beyond the limit of small $\theta$. Evidently, with the
birdie-like alignment mechanism, the timescale to reach alignment is in fact
the damping time, which we discuss in more detail below.

###  2.2 Alignment Direction

When $P\neq 0$ (i.e., when the two grains are not identical and $A\neq 0$), we
can easily understand why this grain has to be aligned by defining the
potential energy from the $\theta$-dependent torque:

$\displaystyle
U\equiv-\int_{0}^{\theta}(-P\sin\theta^{\prime})d\theta^{\prime}=-P\cos\theta$
(13)

where $\theta^{\prime}$ is a dummy variable and the reference potential energy
at $\theta=0$ is $-P$. We first consider $P>0$, in which case $U$ is a minimum
at $\theta=0$. The $-\cos\theta$ dependence of $U$ is easy to understand
through Fig. 2 which shows a schematic of how the orientation of the grain
relates to $U$ described by Eq. 13. The direction of the $\theta$-dependent
torque depends on the orientation of the grain and acts against the
displacement of $\theta$ from $0$ attempting to trap the grain in alignment.
At $\theta=0$, the torque vanishes and the grain is kept stable with sphere 2
following the direction of the flow (the lighter and/or larger sphere follows
$\bm{A}$). The anti-alignment point, $\theta=\pi$, also does not produce a
torque, but is unstable. The $\theta$-dependent torque is thus a restoring
torque within the potential well, which we call the “flow-induced restoring
torque.”

The existence of the damping torque ($\omega$-dependent torque) diminishes the
spin, which eventually decreases the rotational energy until $\theta=0$, i.e.,
the grain becomes aligned, which corresponds to the minimum energy state.
Alternatively, if $P<0$ (for example, sphere 2 becomes the heavier sphere),
then $U$ reaches a minimum at $\theta=\pi$ though the lighter and/or larger
sphere still follows the direction of the flow (recall that $\theta$ is
specifically defined using $\bm{r}_{2}$; see Fig. 1).

Figure 2:  A schematic of how the orientation of the grain relates to the
potential well created by the existence of aerodynamic flow. Top panel: the
orientation of an example double-sphere grain model that satisfies
$m_{1}/m_{2}>\varsigma_{1}^{2}/\varsigma_{2}^{2}$ like in Fig. 1 The heavier
sphere is colored gray, while the lighter sphere is white. The green lines
with arrows denote the direction of the gas flow. The circular arrow depicts
the direction of the flow-induced torque, which always tries to align the
grain and resists any angular displacement. As such, the flow-induced torque
is a “restoring" torque. When aligned, the lighter sphere follows the
direction of the flow, analogous to how the heavier head of the badminton
birdie leads the lighter tail against a headwind. Bottom panel: the energy
potential $U$ as a function of $\theta$. The depth of the potential well is
$2P$. Alignment occurs at $\theta=0$ where $U$ is minimal.

From Eq. 13, the peak of the potential well is simply $2P$. If the total
energy (rotational plus potential energy) of a grain is greater than $2P$, the
grain will not be trapped in a potential well, but continue to spin (and not
oscillate). For an initially aligned grain, one can easily derive that the
escape angular speed is $\omega_{\text{esc}}=2\omega_{o}$, i.e., a grain
spinning faster than $\omega_{\text{esc}}$ cannot be trapped in oscillation.
As gas damping dissipates the rotational energy of the grain, the grain will
inevitably land in a potential well and begin oscillating until it reaches
alignment at $\theta=0$.

Observing the expression of $P$ from Eq. 9, we can realize that the quantity
$-\varsigma_{1}^{2}r_{1}+\varsigma_{2}^{2}r_{2}$ is the first moment of the
cross-sections of the spheres. We can define the “geometrical center” of the
double-sphere grain by

$\displaystyle
r_{g}\equiv\frac{-\varsigma_{1}^{2}r_{1}+\varsigma_{2}^{2}r_{2}}{\varsigma_{1}^{2}+\varsigma_{2}^{2}}$
(14)

and use $\bm{r}_{g}$ to denote the vector from the center of mass to the
geometrical center. When $r_{g}>0$, the geometrical center shifts towards
sphere 2, and the final alignment angle at $\theta=0$ is such that
$\bm{r}_{g}\parallel\bm{A}$. Conversely, when $r_{g}<0$, the geometrical
center moves closer to sphere 1, and the final alignment angle at $\theta=\pi$
also means $\bm{r}_{g}\parallel\bm{A}$. In other words, the asymmetry between
the two spheres causes an offset between the center of mass and the
geometrical center and the grain becomes susceptible to alignment with
$\bm{A}$. The final alignment orientation is such that the geometrical center
points along the direction of the flow relative to the center of mass. When
$r_{g}=0$, the geometrical center corresponds to the center of mass leading to
$P=0$.

Since the center of mass is related to the mass of the spheres, we can express
how alignment depends on the mass or density. To fulfill $P>0$, it requires

$\displaystyle\frac{m_{1}}{m_{2}}>\frac{\varsigma_{1}^{2}}{\varsigma_{2}^{2}}$
(15)

where we used the definition of $P$ and Eq. 1. If the spheres have equal mass,
$m_{1}=m_{2}$, then the first sphere should be smaller than the second sphere
($\varsigma_{2}>\varsigma_{1}$) to reach alignment at $\theta=0$.

Equivalently, we can re-express $P>0$ using the material density and obtain

$\displaystyle\frac{\rho_{1}}{\rho_{2}}\frac{\varsigma_{1}}{\varsigma_{2}}>1.$
(16)

This expression explains that both an asymmetry in the material density and/or
in the size can create alignment. If $\rho_{1}=\rho_{2}$, the first sphere
needs to be larger than the second sphere $\varsigma_{1}>\varsigma_{2}$ to
reach alignment at $\theta=0$ (note the opposite relation because the material
density is fixed as opposed to keeping the total mass fixed).

Lastly, we consider the behavior when the flow-induced torque is non-existent.
From Eq. 9, we can see that $P$ can equal 0 when the two spheres are identical
or when $A=0$, while $D$ is always greater than 0. When $P=0$, the equation
reduces to $I\dot{\omega}+D\omega=0$, which can be integrated once to yield
the time evolution of the grain misalignment angle relative to the flow:

$\displaystyle\theta(t)=\theta(0)+\frac{1}{2}\omega(0)t_{d}(1-e^{-2t/t_{d}})$
(17)

where $\theta(0)$ and $\omega(0)$ are the initial conditions of $\theta$ and
$\omega$. Eq. 17 means that the grain does not oscillate. The final
orientation of $\theta$($t\rightarrow\infty$) is $\theta(0)+\omega(0)t_{d}/2$.
In other words, without the restoring torque ($P=0$), the final alignment
angle depends on the initial conditions. If the initial conditions for the
grains are random and the direction of rotation is also random, then the
ensemble should not have any preferred alignment direction (see Section 5 for
a discussion on the resulting polarization).

###  2.3 Parameterization

While the above discussion allows a qualitative description of alignment, we
utilize a parameterization of the double-sphere model to facilitate the
quantitative description. In particular, we aim to quantify how $t_{o}$ and
$t_{d}$ depend on the characteristic properties of the grain, in addition to
the level of asymmetry between the spheres. We define the characteristic
length $l$ through

$\displaystyle l\equiv r_{1}+r_{2}$ (18)

which is simply the length between the center of the two spheres. We define
the ratio of the sizes through

$\displaystyle\epsilon\equiv\frac{\varsigma_{1}}{\varsigma_{2}}$ (19)

and the sum of the sizes is constrained by

$\displaystyle\varsigma_{1}+\varsigma_{2}\equiv\lambda l$ (20)

where $\lambda$ is the fraction of $l$ that the radius of each sphere occupies
with $\lambda\leq 1$. $\lambda=1$ means the two spheres are in contact. The
sizes are thus

$\displaystyle\varsigma_{1}$
$\displaystyle=\frac{\epsilon\lambda}{1+\epsilon}l,$
$\displaystyle\varsigma_{2}$ $\displaystyle=\frac{\lambda}{1+\epsilon}l.$ (21)

In addition to the size, we define the ratio of the densities through

$\displaystyle\kappa\equiv\frac{\rho_{1}}{\rho_{2}}.$ (22)

We denote the total mass of the grain through $m\equiv m_{1}+m_{2}$ and along
with definitions of $\epsilon$ and $\kappa$, one can find that

$\displaystyle m_{1}$
$\displaystyle=\frac{\epsilon^{3}\kappa}{1+\epsilon^{3}\kappa}m,$
$\displaystyle m_{2}$ $\displaystyle=\frac{1}{1+\epsilon^{3}\kappa}m.$ (23)

Since $r_{i}$ for each sphere is related to the center of mass through Eq. 1,
we can re-express $r_{i}$ by

$\displaystyle r_{1}$ $\displaystyle=\frac{1}{1+\epsilon^{3}\kappa}l,$
$\displaystyle r_{2}$
$\displaystyle=\frac{\epsilon^{3}\kappa}{1+\epsilon^{3}\kappa}l.$ (24)

It is also convenient to consider the overall material density $\rho_{s}$
which is the total mass encompassed by the total volume of the two spheres:

$\displaystyle\rho_{s}=\frac{m}{l^{3}}\frac{3}{4\pi\lambda^{3}}\frac{(1+\epsilon)^{3}}{1+\epsilon^{3}}$
(25)

One can find that

$\displaystyle\rho_{1}$
$\displaystyle=\frac{1+\epsilon^{3}}{1+\epsilon^{3}\kappa}\kappa\rho_{s}$
$\displaystyle\rho_{2}$
$\displaystyle=\frac{1+\epsilon^{3}}{1+\epsilon^{3}\kappa}\rho_{s}$ (26)

We can now determine the coefficients to Eq. 9. The moment of inertia is

$\displaystyle
I=\frac{4\pi}{3}\rho_{s}l^{5}\frac{\lambda^{3}(1+\epsilon^{3})}{(1+\epsilon)^{3}(1+\epsilon^{3}\kappa)}\bigg{[}\frac{\epsilon^{3}\kappa}{1+\epsilon^{3}\kappa}+\frac{2}{5}\lambda^{2}(1+\epsilon^{5}\kappa)\bigg{]}$
(27)

The coefficient to the drag term is

$\displaystyle
D=\frac{4\pi}{3}\rho_{g}v_{\text{th}}l^{4}\frac{\lambda^{2}\epsilon^{2}(1+\epsilon^{4}\kappa^{2})}{(1+\epsilon)^{2}(1+\epsilon^{3}\kappa)^{2}}$
(28)

Lastly, we find

$\displaystyle
P=\frac{4\pi}{3}\rho_{g}v_{\text{th}}Al^{3}\frac{\lambda^{2}\epsilon^{2}(\epsilon\kappa-1)}{(1+\epsilon)^{2}(1+\epsilon^{3}\kappa)}.$
(29)

Note that the factor $\epsilon\kappa-1$ is equivalent to Eq. 16 and determines
the sign of $P$, which, in turn, determines the alignment direction.

Figure 3:  Dimensionless factors for the double-sphere model as a function of
$\epsilon$ (ratio of radii between the two spheres) and $\kappa$ (ratio of
material density). Panel a: the damping time dimensionless factor,
$\hat{t}_{d,d}$. Panel b: the oscillation time dimensionless factor,
$\hat{t}_{o,d}$. As $\epsilon\rightarrow 1$ and $\kappa\rightarrow 1$,
$\hat{t}_{o,d}\rightarrow\infty$ meaning there is no oscillation when the two
spheres are equivalent.

With the coefficient determined, we can assess the damping time $t_{d}$ and
oscillation time $t_{o}$. We will use $t_{d,d}$ and $t_{o,d}$ to denote
$t_{d}$ and $t_{o}$ of the double-sphere model, respectively. The damping time
for the double-sphere is

$\displaystyle t_{d,d}$
$\displaystyle=\frac{\rho_{s}l}{\rho_{g}v_{\text{th}}}\hat{t}_{d,d}\text{ ,}$
$\displaystyle\hat{t}_{d,d}$
$\displaystyle\equiv\frac{2\lambda(1+\epsilon^{3})(1+\epsilon^{3}\kappa)}{\epsilon^{2}(1+\epsilon)(1+\epsilon^{4}\kappa^{2})}\bigg{[}\frac{\epsilon^{3}\kappa}{1+\epsilon^{3}\kappa}+\frac{2}{5}\lambda^{2}(1+\epsilon^{5}\kappa)\bigg{]}$
(30)

where $\hat{t}_{d,d}$ is a dimensionless factor that encapsulates the
asymmetry. We will use the “ $\hat{}$ ” symbol to denote the dimensionless
factors that only depend on the geometry of the grain. The physical quantity
$\rho_{s}l/(\rho_{g}v_{\text{th}})$ is equivalent to the stopping time of a
spherical grain with radius $l$ and material density $\rho_{s}$. The stopping
time is the characteristic timescale for the aerodynamic drag force to stop
the motion of a moving grain (Armitage, 2015). The relation between the
damping time and the stopping time of a single sphere is not too surprising.
As the double-sphere grain spins, each sphere moves relative to the gas.
Requiring each sphere to stop its relative motion with respect to the gas is
equivalent to stopping the spin of the system, i.e., rotational damping. Thus,
the damping time should be on the order of the stopping time. However, this
effect cannot be captured for a single spherical grain as we will see in
Section 3. Plugging in some typical numbers for a protoplanetary disk for the
stopping time , we get

$\displaystyle t_{s}$
$\displaystyle\equiv\frac{\rho_{s}l}{\rho_{g}v_{\text{th}}}$
$\displaystyle\sim 25\text{ years }\bigg{(}\frac{\rho_{s}}{3\text{ g
cm${}^{-3}$ }}\bigg{)}\bigg{(}\frac{l}{1\text{ mm
}}\bigg{)}\bigg{(}\frac{10^{9}\text{ cm${}^{-3}$
}}{n_{g}}\bigg{)}\bigg{(}\frac{1\text{ km/s }}{v_{\text{th}}}\bigg{)}.$ (31)

Since increasing gas density dampens the grain spin more quickly
($t_{d}\propto 1/\rho_{g}$), we expect that grains are better aligned in
higher gas densities. Also, since $t_{d}\propto l$, larger grains are more
difficult to be damped and should be less aligned.

Fig. 3a shows the asymmetry factor $\hat{t}_{d,d}$ for different values of
size ratio $\epsilon$ and density ratio $\kappa$ while adopting $\lambda=1$
(i.e., contacting spheres). We find that the value is of order unity across
the parameter space. When $\epsilon=1$ and $\kappa=1$ (the symmetric case),
then $\hat{t}_{d,d}=1.4$. Note that the flow speed $A$ does not contribute to
the damping time, but only determines the oscillation potential (and hence
frequency) around the dust-gas drift direction.

The oscillation time for the double-sphere is

$\displaystyle t_{o,d}$
$\displaystyle=\sqrt{\frac{\rho_{s}l^{2}}{\rho_{g}v_{\text{th}}A}}\hat{t}_{o,d}\text{
,}$ $\displaystyle\hat{t}_{o,d}$ $\displaystyle\equiv
2\pi\sqrt{\frac{\lambda(1+\epsilon^{3})}{\epsilon^{2}(1+\epsilon)(\epsilon\kappa-1)}\bigg{[}\frac{\epsilon^{3}\kappa}{1+\epsilon^{3}\kappa}+\frac{2}{5}\lambda^{2}(1+\epsilon^{5}\kappa)\bigg{]}}$
(32)

where $\hat{t}_{o,d}$ is a dimensionless factor for the oscillation time and
the physical quantity in front of it is the characteristic oscillation time.
Plugging in the same numbers for a protoplanetary disk as Eq. 2.3 and adopting
$A=10$ m/s (see Sec. 4 below), the characteristic oscillation time is $\sim 5$
minutes.

Fig. 3b shows $\hat{t}_{o,d}$ as a function of size ratio $\epsilon$ for
different values of density ratio $\kappa$. While most of the parameter space
shows $\hat{t}_{o,d}\sim 10$, its value increases rapidly as $\epsilon$ and
$\kappa$ decreases to 1. Nevertheless, even if $\kappa=1.01$ and $\epsilon=1$,
or $\kappa=1$ and $\epsilon=1.01$, $\hat{t}_{o,d}\sim 50$. Thus, $t_{o,d}$ is
much shorter than the damping time and the Keplerian orbital periods of the
outer disk even for $1\%$ asymmetry.

Since $t_{o,d}\ll t_{d,d}$, the grain undergoes many oscillations before the
amplitude of oscillation diminishes. The oscillation behavior is more akin to
an underdamped harmonic oscillator which enters its equilibrium state (i.e.,
becomes aligned) gradually. We can understand why through the damping ratio
defined by $\zeta\equiv D/(2\sqrt{IP})$ where $\zeta<1$, $\zeta=1$, and
$\zeta>1$ correspond to underdamped, critically damped, and overdamped
oscillators, respectively. Applying the constants from Eq. 9, we have

$\displaystyle\zeta_{d}$
$\displaystyle=\sqrt{\frac{\rho_{g}v_{\text{th}}}{\rho_{s}A}}\hat{\zeta}_{d}\text{
, }$ $\displaystyle\hat{\zeta}_{d}$
$\displaystyle\equiv\frac{\epsilon(1+\epsilon^{4}\kappa^{2})}{2(1+\epsilon^{3}\kappa)}\sqrt{\frac{(1+\epsilon)^{3}}{\lambda(1+\epsilon^{3})(\epsilon\kappa-1)}}\bigg{[}\frac{\epsilon^{3}\kappa}{1+\epsilon^{3}\kappa}+\frac{2}{5}\lambda^{2}(1+\epsilon^{5}\kappa)\bigg{]}^{-\frac{1}{2}}$
(33)

where $\hat{\zeta}_{d}$ encapsulates the dimensionless, asymmetry factors.
Using typical conditions of disks used in Eq. 2.3, we have $\zeta_{d}\sim
4\times 10^{-7}\hat{\zeta}_{d}$ and $\hat{\zeta}_{d}\sim 6$ for
$\epsilon=1.01$ and $\kappa=1$. One can see that $\zeta_{d}\ll 1$ because the
gas mass density $\rho_{g}$ is drastically smaller than grain material density
$\rho_{s}$ (by 15 orders of magnitude with $n_{g}=10^{9}$ cm-3), which is more
than enough to compensate for the difference between the dust-gas relative
speed $A$ and the thermal speed $v_{\text{th}}$.

Fig. 4 shows the phase portrait (or trajectory maps) of the damped
oscillations given by Eq. 9. The phase portrait shows the evolution of
$\theta$ and $\omega$ for the dynamical system333We use the Python package
phaseportrait available at https://phaseportrait.github.io. We used the same
values from Eq. 2.3 along with $\epsilon=1.01$, $\kappa=1$, and $A=10$ m/s.
Fig. 4 also shows the separatrix described by
$(\omega/\omega_{0})^{2}=2(1+\cos\theta)$ which separates the bounded and
unbounded region. One can derive this relation by equating the total energy
(kinetic and potential energy) of the system to the energy at the unstable
point. Grains in the unbounded regions continually spin without regard to the
alignment direction. However, as the grain loses energy through damping, it
will cross the separatrix at some point and become bound within a potential
well (see Appendix A for an example phase portrait when the grain is more
damped). We can see that $\theta$ of a grain oscillates around the alignment
point at $\theta=0$ and is repelled at $\theta=-\pi$ and $\pi$.

Figure 4:  The phase portrait (trajectory map) of the double-sphere grain
under the presence of gas flow $\bm{A}$. The horizontal axis is $\theta$,
while the vertical axis is the angular velocity $\omega$ normalized by the
oscillation frequency $\omega_{o}$. The black circles correspond to stable
equilibrium points that serve as the attractors for alignment. The orange plus
signs are the unstable equilibrium points. The blue curve is the separatrix
that separates the bounded and unbounded regions. Regions inside the
separatrix are bounded and oscillate around the attractors.

To conclude the results from the double-sphere model, the asymmetry in the
structure (a combination of the cross-section and mass) and the existence of
gas flow lead to restoring torques that make the grain oscillate instead of a
spin-up. As the grain oscillates, the torque from gas damping opposes any
grain rotating motion which dampens the magnitude of oscillation. Even a $1\%$
level of asymmetry allows the grain to undergo severely underdamped
oscillations. The damped oscillations eventually reach an orientation where
the grain geometric center aligns along the gas flow with respect to its
center of mass. The timescale to align the grain is related to the damping
timescale which is comparable to the stopping time of the entire grain. In
other words, the grain asymmetry and gas flow determine the direction of
alignment and the damping time determines how quickly the grain reaches
alignment.

##  3 Spheroid Model

The double-sphere model of a grain allowed estimates of the oscillation
behavior depending on the level of asymmetry and the timescales for alignment.
In addition, the model will always be aligned along the long axis connecting
the two spheres, which hints at a potential solution that can match the
requirement of effectively aligned prolate grains constrained from
observations. However, the shape of the model is too idealized to assess its
polarization from thermal emission. In this section, we utilize an
axisymmetric ellipsoid, i.e., spheroid (prolate and oblate), and assess the
alignment behavior. The polarization from spheroids is well established, both
in the dipole limit (van de Hulst, 1957) and for arbitrary sizes using the
T-matrix technique (Waterman, 1971; Mishchenko & Travis, 1994). The spheroidal
alignment model can be directly compared to observed polarization.

###  3.1 Problem Setup

Similar to the previous section, we aim to solve the total torque
$\bm{\Gamma}$ of a spinning grain embedded in a flowing background of gas.
Obtaining $\bm{\Gamma}$ will allow us to derive the equation of motion. While
the double-sphere model only adds torque from two discrete points (the centers
of the two spheres), a smooth object, like a spheroid, requires us to
integrate the torque over its surface $S$:

$\displaystyle\bm{\Gamma}=\oint_{S}\bm{r}\times d\bm{F}$ (34)

where $\bm{r}$ is the radius vector of a surface element $d\sigma$ from the
center of mass $O$ and $d\bm{F}$ is the infinitesimal force from the gas
acting on the surface element.

Assuming specular reflection and subsonic relative motion, Appendix B shows
that $d\bm{F}$ at each point of the surface follows

$\displaystyle
d\bm{F}=-\rho_{g}\bigg{[}c_{s}^{2}\bm{n}+v_{\text{th}}\textbf{{N}}(\bm{u})\bigg{]}d\sigma$
(35)

where $\bm{n}$ is the unit normal direction of $d\sigma$ and $\bm{u}$ is the
velocity of the surface relative to the gas. $\textbf{{N}}\equiv\bm{n}\bm{n}$
is a projection tensor that simply returns a vector from the part of $\bm{u}$
that is normal to the surface $d\sigma$. The negative sign means the force is
opposite to the direction of $\bm{n}$. $c_{s}$ is the isothermal sound speed
$c_{s}\equiv\sqrt{kT/(\mu m_{p})}$.

Assume the grain spins as a solid body with angular velocity $\bm{\omega}$
around the center of mass. Each point of the surface will travel relative to
the gas as

$\displaystyle\bm{u}=\bm{\omega}\times\bm{r}-\bm{A}.$ (36)

The expression is directly related to Eq. 5, but in this case, the $\bm{u}$ is
the relative velocity of a surface element and not a complete sphere that was
assumed for the double-sphere model. Substituting Eq. 36 and Eq. 35 into Eq.
34 yields three terms to the total torque on the spheroid:

$\displaystyle\bm{\Gamma}=-\rho_{g}c_{s}^{2}\bm{K}+\rho_{g}v_{\text{th}}\textbf{{L}}\bm{\omega}+\rho_{g}v_{\text{th}}\textbf{{M}}\bm{A}$
(37)

with

$\displaystyle\bm{K}$ $\displaystyle\equiv\oint_{S}\bm{r}\times\bm{n}d\sigma,$
L $\displaystyle\equiv\oint_{S}\bm{r}\times\textbf{{N}}(\bm{r}\times)d\sigma,$
M $\displaystyle\equiv\oint_{S}\bm{r}\times\textbf{{N}}d\sigma$

where $\bm{K}$ is a vector and L and M are tensors. We have rearranged the
expression such that only quantities that depend on the location of the
surface are kept inside each integral. Note that the cross product,
$\bm{r}\times$, is a 2-rank tensor that can be implemented as a matrix
calculation. The first term is a driving torque that exists regardless of any
$\bm{A}$ or $\bm{\omega}$. The second term is the damping torque since it is
related to $\bm{\omega}$. The third term depends on $\bm{A}$ and, as we will
see later, produces the restoring torque.

Eq. 37 is particularly useful since each integral only depends on the
structure of the grain and does not depend on the environmental properties,
like $\bm{A}$ or $\rho_{g}$, or its dynamic state, like $\theta$ or
$\bm{\omega}$. This opens the possibility of precalculating these integrals
for any grain as its unique physical property before considering its dynamic
behavior.

###  3.2 Strictly Axisymmteric Spheroid

Thus far, we have not utilized any assumption about the geometry of the
spheroid. We now evaluate each term in Eq. 37 using a dedicated coordinate
system shown in Fig. 5. The setup is equivalent to Sec. 2. Let $\bm{e}_{1}$,
$\bm{e}_{2}$, and $\bm{e}_{3}$ represent the Cartesian, unit directions of the
lab frame that is centered at the origin $O$. We define $\bm{e}_{3}$ to be
along the direction of gas flow $\bm{A}$ and let the grain rotate around
$\bm{e}_{2}$.

Figure 5:  Schematic of a prolate and an oblate where the center of mass is
shifted along the axis of symmetry. We use ellipses to represent the cross-
sections of the two spheroids in the $e_{1}e_{3}$-plane. For the prolate, the
axis of symmetry is the long axis, while for the oblate, the axis of symmetry
is the short axis. $\bm{e}_{1}$, $\bm{e}_{2}$, and $\bm{e}_{3}$ are lab frame
unit vectors where $\bm{e}_{2}$ points out of the page. $\bm{b}_{1}$,
$\bm{b}_{2}$, and $\bm{b}_{3}$ are body frame unit vectors, where $\bm{b}_{2}$
also points out of the page and $\bm{b}_{3}$ is along the axis of symmetry.
The offset of the center of the spheroid from the center of mass $O$ is
$\bm{g}$. Since the center of mass is shifted along the axis of symmetry,
$\bm{g}\parallel\bm{b}_{3}$. The green arrow is the direction of the gas flow
$\bm{A}$.

Let the body frame of the spheroid also be a Cartesian coordinate system with
its center denoted by $G$. We use $\bm{b}_{1}$, $\bm{b}_{2}$, and $\bm{b}_{3}$
to denote the unit vectors. $\bm{b}_{3}$ is along the axis of asymmetry and
$\bm{b}_{2}$ is parallel to $\bm{e}_{2}$. The angle between $\bm{b}_{3}$ and
$\bm{e}_{3}$ is $\theta$. The relation between these two sets of unit vectors
is simply

$\displaystyle\bm{b}_{i}=W_{ij}\bm{e}_{j}$ (38)

where we adopt the Einstein summation convention and $i$ and $j$ are indices
from 1 to 3. $W_{ij}$ are the elements of the matrix

$\displaystyle\textbf{{W}}=\begin{pmatrix}\cos\theta&0&-\sin\theta\\\ 0&1&0\\\
\sin\theta&0&\cos\theta\end{pmatrix}.$ (39)

In the body frame, the surface of the spheroid follows

$\displaystyle\frac{x^{2}+y^{2}}{a^{2}}+\frac{z^{2}}{c^{2}}=1$ (40)

where $x$, $y$, and $z$ are coordinates in the body frame along $\bm{b}_{1}$,
$\bm{b}_{2}$, and $\bm{b}_{3}$, respectively. In addition, $c$ is the length
of the spheroid along the axis of symmetry (not to be confused with the
isothermal sound speed $c_{s}$ or the speed of light) and $a$ is the length
along the axis perpendicular to that. A prolate has $c>a$, while an oblate has
$c<a$.

Considering axisymmetric particles makes the surface integral easier to
compute via a surface of revolution. We parameterize the surface by
introducing the azimuthal angle $\phi$ around the $\bm{b}_{3}$ axis

$\displaystyle\begin{pmatrix}x\\\
y\end{pmatrix}=\begin{pmatrix}w(z)\cos\phi\\\ w(z)\sin\phi\end{pmatrix}$ (41)

where $\phi$ starts from $\bm{b}_{1}$ and increases following the right-hand
rule. Furthermore, $w(z)$ is the perpendicular distance from the surface to
the axis of symmetry and serves as the generatrix for the surface integration:

$\displaystyle w(z)\equiv a\sqrt{1-\bar{z}^{2}}$ (42)

where $\bar{z}\equiv z/c$ is $z$ normalized by $c$. We will use the bar symbol
to denote length quantities normalized by $c$. With the generatrix, we can
derive the surface element to be

$\displaystyle d\sigma=ac\sqrt{1+(\bar{a}^{2}-1)\bar{z}^{2}}d\bar{z}d\phi$
(43)

where $\bar{a}\equiv a/c$.

The normal vector of the spheroid is related to the gradient of Eq. 40. We
express the unit normal vector by $\bm{n}=\tilde{n}_{i}\bm{b}_{i}$:

$\displaystyle\begin{pmatrix}\tilde{n}_{1}\\\ \tilde{n}_{2}\\\
\tilde{n}_{3}\end{pmatrix}=\frac{1}{\sqrt{1+(\bar{a}^{2}-1)\bar{z}^{2}}}\begin{pmatrix}\sqrt{1-\bar{z}^{2}}\cos\phi\\\
\sqrt{1-\bar{z}^{2}}\sin\phi\\\ \bar{a}\bar{z}\end{pmatrix}\text{ .}$ (44)

We use a tilde to denote components using $\bm{b}_{i}$ as basis vectors.

The radius vector of the surface from $O$ can be decomposed to

$\displaystyle\bm{r}=\bm{g}+\bm{s}$ (45)

where $\bm{g}$ points from the center of mass to the geometric center $G$ and
$\bm{s}$ points from $G$ to the surface. The components of $\bm{s}$ using
$\bm{b}_{i}$ basis vectors are

$\displaystyle\begin{pmatrix}\tilde{s}_{1}\\\ \tilde{s}_{2}\\\
\tilde{s}_{3}\end{pmatrix}=\begin{pmatrix}a\sqrt{1-\bar{z}^{2}}\cos\phi\\\
a\sqrt{1-\bar{z}^{2}}\sin\phi\\\ c\bar{z}\end{pmatrix}$ (46)

For the offset vector $\bm{g}$, we first consider the offset along the axis of
symmetry to follow the strict assumption of axisymmetry. If the length between
$O$ and $G$ is $g$, we have $\bm{g}=g\bm{b}_{3}$. We later relax this
assumption and explore the effects of $\bm{g}$ when it is not along the axis
of symmetry.

With the relevant details of the surface defined, we can now solve for
$\bm{K}$, L, and M. When integrating over the surface, we find a common
occurrence of an integral in the form of

$\displaystyle
E[f(\bar{z})]\equiv\int_{-1}^{1}\frac{f(\bar{z})}{\sqrt{1+(\bar{a}^{2}-1)\bar{z}^{2}}}d\bar{z}.$
(47)

where $f(\bar{z})$ is some function that depends on $\bar{z}$. Although
$\bar{z}$ serves as a dummy variable for integration under this definition, it
originates from the coordinate along the symmetry axis of the spheroid as we
integrate over the surface. The definition of $E[f(\bar{z})]$ implicitly
depends on $\bar{a}$ for easier notation. Appendix C shows the analytical
solutions to the $E$-integrals used in this paper.

For $\bm{K}$, we make use of Eq. 45 for $\bm{r}$ and obtain:

$\displaystyle\bm{K}=\bm{g}\times\oint_{S}\bm{n}d\sigma+\oint_{S}\bm{s}\times\bm{n}d\sigma=0.$
(48)

Note that the integral in the first term, $\oint_{S}\bm{n}d\sigma$, is always
$\bm{0}$ for a closed surface444Under the current context, this is true since
we assumed that the entire surface of a grain reflects gas equally through Eq.
35. If certain patches of the grain reflect gas differently, then there is an
interesting possibility that the integral may not have to be $\bm{0}$. . The
latter term integrates to 0 for a spheroid.

For L, using Eq. 45 gives four terms:

$\displaystyle\textbf{{L}}=\bm{g}\times\oint_{S}\textbf{{N}}d\sigma(\bm{g}\times)+\bm{g}\times\oint_{S}\textbf{{N}}(\bm{s}\times)d\sigma+\oint_{S}\bm{s}\times\textbf{{N}}d\sigma(\bm{g}\times)$
$\displaystyle+\oint_{S}\bm{s}\times\textbf{{N}}(\bm{s}\times)d\sigma.$ (49)

It is convenient to know that the integral of N expressed in $\bm{b}_{i}$
basis vectors is

$\displaystyle\oint_{S}\textbf{{N}}d\sigma=\pi
c^{2}\begin{pmatrix}\bar{a}E[1-\bar{z}^{2}]&0&0\\\
0&\bar{a}E[1-\bar{z}^{2}]&0\\\ 0&0&2\bar{a}^{3}E[\bar{z}^{2}]\end{pmatrix}.$
(50)

For the first term of Eq. 3.2, using the $\bm{b}_{i}$ basis vectors, we find
that

$\displaystyle\bm{g}\times\oint_{S}\textbf{{N}}d\sigma(\bm{g}\times)=-\pi
c^{4}\bar{a}\bar{g}^{2}E[1-\bar{z}^{2}]\begin{pmatrix}1&0&0\\\ 0&1&0\\\
0&0&0\end{pmatrix}$ (51)

where $\bar{g}\equiv g/c$. Each element in the second term integrates to $0$.
Since N is a symmetric 2-rank tensor and both $\bm{s}\times$ and
$\bm{g}\times$ are anti-symmetric 2-rank tensors, the third term is equivalent
to the transpose of the second term. Lastly, the fourth term is a symmetric
tensor, which means there are only six independent elements when expressed as
a matrix in general. For the spheroid, four of the independent elements
integrate to 0 leaving only two non-zero elements:

$\displaystyle\oint_{S}\bm{s}\times\textbf{{N}}(\bm{s}\times)d\sigma=-\pi
c^{4}\bar{a}(1-\bar{a}^{2})^{2}E[\bar{z}^{2}-\bar{z}^{4}]\begin{pmatrix}1&0&0\\\
0&1&0\\\ 0&0&0\end{pmatrix}$ (52)

Thus, we find that

$\displaystyle\tilde{L}_{11}=\tilde{L}_{22}=-\pi
c^{4}\bar{a}\bigg{[}\bar{g}^{2}E[1-\bar{z}^{2}]+(1-\bar{a}^{2})^{2}E[\bar{z}^{2}-\bar{z}^{4}]\bigg{]}$
(53)

while the rest of the elements are zero.

Lastly, for M, using the $\bm{b}_{i}$ basis vectors, we have

M
$\displaystyle=\bm{g}\times\oint_{S}\textbf{{N}}d\sigma+\oint_{S}\bm{s}\times\textbf{{N}}d\sigma$
$\displaystyle=\pi c^{3}\bar{a}\bar{g}E[1-\bar{z}^{2}]\begin{pmatrix}0&-1&0\\\
1&0&0\\\ 0&0&0\end{pmatrix}$ (54)

where we found that the second term is 0.

We can now derive the torque from Eq. 37. Given the assumed
$\bm{A}=A\bm{e}_{3}$, we can express $\bm{A}=\tilde{A}_{i}\bm{b}_{i}$ and find
that

$\displaystyle\bm{A}=-A\sin\theta\bm{b}_{1}+A\cos\theta\bm{b}_{3}$ (55)

using Eq. 39. Since we only consider rotation around $\bm{e}_{2}$ and
$\bm{e}_{2}=\bm{b}_{2}$, the angular velocity is simply
$\bm{\omega}=\omega\bm{e}_{3}=\omega\bm{b}_{2}$. Applying to Eq. 37, the
torque is thus

$\displaystyle\tilde{\Gamma}_{2}=$ $\displaystyle-\rho_{g}v_{\text{th}}\pi
c^{4}\omega\bar{a}\bigg{[}\bar{g}^{2}E[1-\bar{z}^{2}]+(1-\bar{a}^{2})^{2}E[\bar{z}^{2}-\bar{z}^{4}]\bigg{]}$
$\displaystyle-\rho_{g}v_{\text{th}}A\pi
c^{3}\bar{a}\bar{g}E[1-\bar{z}^{2}]\sin\theta$ (56)

and we find that $\tilde{\Gamma}_{1}=\tilde{\Gamma}_{3}=0$. Note that
$\Gamma_{2}=\tilde{\Gamma}_{2}$ also because $\bm{e}_{2}=\bm{b}_{2}$. Thus, we
see that the torque is only limited along $\bm{e}_{2}$ which is what we expect
when we only consider $\bm{\omega}$ along $\bm{e}_{2}$.

The moment of inertia around the axis passing through the geometric center $G$
and perpendicular to the axis of symmetry is

$\displaystyle
I_{\perp}=\frac{4\pi}{15}c^{5}\rho_{s}\bar{a}^{2}(1+\bar{a}^{2})$ (57)

for a homogenous spheroid with material density $\rho_{s}$. However, to shift
the center of mass $O$ away from the geometric center $G$, the spheroid cannot
be homogenous. Nevertheless, based on the double-sphere model, the necessary
shift is relatively small compared to the size of the grain, so the homogenous
sphere remains a good approximation even in this case. Using $I_{\perp}$ and
the parallel axis theorem, we have

$\displaystyle
I_{s}=\frac{4\pi}{15}c^{5}\rho_{s}\bar{a}^{2}(1+\bar{a}^{2}+5\bar{g}^{2})$
(58)

where we can see that $I_{s}$ corresponds to $I_{\perp}$ when $\bar{g}\ll 1$.
As an approximation, we simply adopt $I_{s}$ as the moment of inertia for the
spheroid. The subscript “s” stands for “spheroid.”

The equation of motion has the same form as the double-sphere model:

$\displaystyle I_{s}\dot{\omega}+D_{s}\omega+P_{s}\sin\theta=0$ (59)

where

$\displaystyle D_{s}$ $\displaystyle=\rho_{g}v_{\text{th}}\pi
c^{4}\bar{a}\bigg{[}\bar{g}^{2}E[1-\bar{z}^{2}]+(1-\bar{a}^{2})^{2}E[\bar{z}^{2}-\bar{z}^{4}]\bigg{]}$
$\displaystyle P_{s}$ $\displaystyle=\rho_{g}v_{\text{th}}A\pi
c^{3}\bar{a}\bar{g}E[1-\bar{z}^{2}]$

meaning the spheroid will also behave as a damped oscillator.

We find that the potential energy is $U_{s}=-P_{s}\cos\theta$, which is again
proportional to $-\cos\theta$ and means the final alignment angle is at
$\theta=0$ for $P_{s}>0$ or at $\theta=\pi$ for $P_{s}<0$. Since
$E[1-\bar{z}^{2}]$ is always positive for all $\bar{a}>0$, the sign of $P_{s}$
comes entirely from the sign of $g$ and is independent of $\bar{a}$. In other
words, both prolate and oblate grains will become aligned with their axis of
symmetry along the gas flow, such that the direction from the center of mass
to the geometric center follows the direction of the gas flow
($\bm{g}\parallel\bm{A}$).

Given Eq. 59 we can explore the damping time $t_{d}$ and oscillation time
$t_{o}$ analytically like in the previous section. We find that

$\displaystyle t_{d,s}$
$\displaystyle=\frac{\rho_{s}c}{\rho_{g}v_{\text{th}}}\hat{t}_{d,s},$
$\displaystyle\hat{t}_{d,s}$
$\displaystyle\equiv\frac{8}{15}\frac{\bar{a}(1+\bar{a}^{2}+5\bar{g}^{2})}{\bar{g}^{2}E[1-\bar{z}^{2}]+(1-\bar{a}^{2})^{2}E[\bar{z}^{2}-\bar{z}^{4}]}.$
(60)

The damping time $t_{d,s}$ is once again related to the stopping time modified
by a dimensionless quantity that depends on the structure of the grain
$\hat{t}_{d,s}$ (through $\bar{a}$ and $\bar{g}$). Fig. 6a shows how
$\hat{t}_{d,s}$ depends on $\bar{a}$ and $\bar{g}$. If we have a sphere
($\bar{a}=1$) rotating about its geometric center ($\bar{g}=0$),
$\hat{t}_{d,s}\rightarrow\infty$ (or equivalently, $D_{s}\rightarrow 0$)
meaning it will not be damped by the gas at all. The result makes sense since
the normal direction at each point of the sphere is always directed radially
from the center of the sphere and thus the gas cannot produce any torque as
the sphere spins. However, if the center of mass of the sphere does not
correspond to the geometrical center ($\bar{g}\neq 0$), then $t_{d,s}$ is no
longer infinite and its rotation can be damped. If we have a homogenous
spheroid ($\bar{a}\neq 1$) rotating about its geometric center ($\bar{g}=0$),
$\hat{t}_{d,s}$ is also finite since the normal directions of most of the
surface elements are no longer strictly in the radial direction from $O$ and
gas can produce a torque. In general, as the non-sphericity ($|\bar{a}-1|$)
increases, $\hat{t}_{d,s}$ decreases. In other words, more non-spherical
grains (for both prolates and oblates) are better damped by the gas and are
aligned faster. With $\bar{g}=0.01$ and $\bar{a}=0.9$, $\hat{t}_{d,s}\sim 80$.

Figure 6:  Dimensionless factors for the strictly axisymmetric spheroid as a
function of $\bar{g}$ (the normalized distance from the center of mass to the
spheroid center, $g/c$) and $\bar{a}$ (the aspect ratio, $a/c$). Panel a: the
damping time dimensionless factor, $\hat{t}_{d,d}$. Panel b: the oscillation
time dimensionless factor, $\hat{t}_{o,d}$. Black solid lines correspond to
the spherical case, $\bar{a}=1$. Dashed lines correspond to prolates
($\bar{a}<1$), while dash-dotted lines correspond to oblates ($\bar{a}>1$).

The oscillation time is

$\displaystyle t_{o,s}$
$\displaystyle=\sqrt{\frac{\rho_{s}c^{2}}{\rho_{g}v_{\text{th}}A}}\hat{t}_{o,s}$
$\displaystyle\hat{t}_{o,s}$
$\displaystyle\equiv\sqrt{\frac{16\pi^{2}}{15}\frac{\bar{a}(1+\bar{a}^{2}+5\bar{g}^{2})}{\bar{g}E[1-\bar{z}^{2}]}}$
(61)

which is also similar to the expression for the double-sphere. Fig. 6b shows
how $\hat{t}_{o,s}$ depends on the structure. For a given $\bar{a}$,
increasing $\bar{g}$ decreases $t_{o,s}$ for both prolates and oblates, which
means the key to allowing a spheroid to oscillate lies in the offset between
the center of mass from the geometric center. For $\bar{g}\rightarrow 0$,
$\hat{t}_{o,s}\rightarrow\infty$ (or equivalently, $P_{s}\rightarrow 0$)
meaning the grain will not oscillate without the offset. Fig. 6b also shows
that increasing $\bar{a}$ (more oblong) increases $\hat{t}_{o,d}$. However,
the large difference is in part due to keeping the length along the axis of
symmetry, $c$, fixed (long axis for a prolate, but short axis for an oblate).
With $\bar{g}=0.01$ and $\bar{a}=0.9$, $\hat{t}_{o,s}\sim 35$.

###  3.3 Quasi-Axisymmetric Spheroid

Previously, we have enforced strict axisymmetry which requires the center of
mass to be along the axis of symmetry. As a result, both prolate and oblate
cases will be aligned with their axis of symmetry along the gas flow. Since
polarization observations of disks favor effectively prolate grains, we would
naively continue our discussion with prolate grains. However, if oblate grains
are aligned along some long axis and the direction of the short axis is random
as an ensemble, they can also appear effectively prolate around the azimuth of
a disk if the dust-gas drift velocity is azimuthal. We are particularly
interested in verifying if oblate grains can be aligned to the gas drift along
the long axis if $\bm{g}$ is along the long axis.

The benefits of axisymmetry lie in the integration of the surface and $\bm{g}$
itself does not participate in the integration. In other words, we can allow
the center of mass to be anywhere in the spheroid, while keeping most of the
derivation the same. If $\bm{g}$ is not along the axis of symmetry, the grain
is no longer strictly axisymmetric in terms of the distribution of mass.
However, since the surface remains axisymmetric, we describe the model as a
“quasi-axisymmetric” model. We show that both the shape parameter $\bar{a}$
and the offset vector $\bm{g}$ determine the final alignment direction. In the
case where $\bm{g}$ lies along the long axis of an oblate, we find that the
long axis of the oblate will be aligned with the gas flow.

Following the same notation above, consider $\bm{g}$ that forms an angle
$\psi$ from $\bm{b}_{3}$. Fig. 7a shows a schematic of a quasi-axisymmetric
oblate. We pick $\bm{b}_{1}$ to be in the same plane as $\bm{g}$ and
$\bm{b}_{3}$, which gives

$\displaystyle\begin{pmatrix}\tilde{g}_{1}\\\ \tilde{g}_{2}\\\
\tilde{g}_{3}\end{pmatrix}=\begin{pmatrix}g\sin\psi\\\ 0\\\
g\cos\psi\end{pmatrix}$ (62)

using $\bm{b}_{i}$ as basis vectors. The gas flow remains in the $\bm{e}_{3}$
direction. We can choose $\bm{e}_{1}$ such that $\bm{g}$ lies in the plane
formed by $\bm{e}_{1}$ and $\bm{e}_{3}$. $\bm{e}_{2}$ remains parallel to
$\bm{b}_{2}$.

Figure 7:  Panel a: a schematic of the coordinate system used to describe the
quasi-axisymmetric oblate following the same coloring scheme and notation as
Fig. 5. Without strictly requiring $\bm{g}$ to be parallel to $\bm{b}_{3}$, it
can form an angle $\psi$ from $\bm{b}_{3}$. The green dotted line helps
visualize the angle $\theta$ formed by $\bm{b}_{3}$ and $\bm{e}_{3}$. The blue
dotted line is in the direction of $\bm{g}$ which helps visualize $\psi$.
Panel b: a schematic of a quasi-axisymmetric oblate with
$\bm{g}\parallel\bm{b}_{1}$ (or equivalently $\psi=90^{\circ}$) when reaching
alignment at $\theta=-90^{\circ}$.

We first evaluate $\bm{K}$, L, and M using Eq. 62. We find that
$\bm{K}=\bm{0}$. For L, one only needs to re-evaluate the first term that
depends on $\bm{g}$ and is non-zero. We find that

$\displaystyle\tilde{L}_{11}$ $\displaystyle=-\pi
c^{4}\bar{a}\bar{g}^{2}E[1-\bar{z}^{2}]\cos^{2}\psi$
$\displaystyle\tilde{L}_{13}=\tilde{L}_{31}$ $\displaystyle=-\pi
c^{4}\bar{a}\bar{g}^{2}E[1-\bar{z}^{2}]\cos\psi\sin\psi$
$\displaystyle\tilde{L}_{22}$ $\displaystyle=-\pi
c^{4}\bar{a}\bar{g}^{2}\bigg{[}E[1-\bar{z}^{2}]\cos^{2}\psi+2\bar{a}^{2}E[\bar{z}^{2}]\sin^{2}\psi\bigg{]}$
$\displaystyle\tilde{L}_{33}$ $\displaystyle=-\pi
c^{4}\bar{a}\bar{g}^{2}E[1-\bar{z}^{2}]\sin^{2}\psi$ (63)

while the rest of the elements of L are 0. When $\psi=0$, we recover Eq. 53.
For M, we find that

$\displaystyle\tilde{M}_{12}$ $\displaystyle=-\pi
c^{3}\bar{a}\bar{g}E[1-\bar{z}^{2}]\cos\psi$ $\displaystyle\tilde{M}_{21}$
$\displaystyle=\pi c^{3}\bar{a}\bar{g}E[1-\bar{z}^{2}]\cos\psi$
$\displaystyle\tilde{M}_{23}$ $\displaystyle=-2\pi
c^{3}\bar{a}^{3}\bar{g}E[\bar{z}^{2}]\sin\psi$ $\displaystyle\tilde{M}_{32}$
$\displaystyle=\pi c^{3}\bar{a}\bar{g}E[1-\bar{z}^{2}]\sin\psi$ (64)

while the rest of the elements of M are $0$. When $\psi=0$, we recover Eq.
3.2. Lastly, the moment of inertia remains equal to Eq. 58.

The equation of motion becomes

$\displaystyle
I_{s}\dot{\omega}+D_{q}\omega+P_{q}\sin\theta+P_{c}\cos\theta=0$ (65)

where the coefficients are

$\displaystyle D_{q}$ $\displaystyle=\rho_{g}v_{\text{th}}\pi
c^{4}\bar{a}\bigg{[}(1-\bar{a}^{2})^{2}E[\bar{z}^{2}-\bar{z}^{4}]$
$\displaystyle+\bar{g}^{2}E[1-\bar{z}^{2}]\cos^{2}\psi$
$\displaystyle+2\bar{a}^{2}\bar{g}^{2}E[\bar{z}^{2}]\sin^{2}\psi\bigg{]}$
$\displaystyle P_{q}$ $\displaystyle=\rho_{g}v_{\text{th}}A\pi
c^{3}\bar{a}\bar{g}E[1-\bar{z}^{2}]\cos\psi$ $\displaystyle P_{c}$
$\displaystyle=2\rho_{g}v_{\text{th}}A\pi
c^{3}\bar{a}^{3}\bar{g}E[\bar{z}^{2}]\sin\psi\text{ .}$

The subscript “q” refers to the quasi-axisymmetric model, while the subscript
“c” for $P_{c}$ refers to the coefficient of its cosine term (“$P_{q,c}$”
would be too redundant). The existence of the extra $\cos\theta$ term in Eq.
65 alters the dynamical behavior giving a potential energy of

$\displaystyle U_{q}=-P_{q}\cos\theta+P_{c}\sin\theta.$ (66)

We can derive the alignment angle by finding $\theta_{\text{align}}$ such that
$U$ is a minimum giving

$\displaystyle\theta_{\text{align}}$
$\displaystyle=\text{atan2}\bigg{(}\frac{-P_{c}}{P_{q}}\bigg{)}$
$\displaystyle=\text{atan2}\bigg{(}\frac{-2\bar{a}^{2}E[\bar{z}^{2}]\sin\psi}{E[1-\bar{z}^{2}]\cos\psi}\bigg{)}$
(67)

where atan2 is the arctangent that gives an angle in the correct quadrant by
taking into account the sign of the arguments. The alignment angle does not
depend on the length of $\bm{g}$, but only on $\psi$ and $\bar{a}$. Fig. 8a
shows the alignment angle as a function of $\psi$ for different values of
$\bar{a}$. If $\psi=0$, i.e., $\bm{g}$ follows the axis of symmetry, then
$P_{c}=0$ and $\theta_{\text{align}}=0$ which recovers the previous results.
If $\psi=\pi/2$, i.e., $\bm{g}$ follows the axis perpendicular to the axis of
symmetry for both prolates and oblates, then $P_{q}=0$ and
$\theta_{\text{align}}=-90^{\circ}$. This means that an oblate will be aligned
to $\bm{A}$ along its long axis if its center of mass is shifted along the
long axis. A schematic of this conceptually important case is shown in Fig. 7b
for an oblate. In fact, if the center of mass is preferentially shifted toward
the long axis of the oblate ($\psi>\sim 45^{\circ}$), then the grain is
preferentially aligned along the long axis which is permitted by disk
polarization when considering the ensemble. Eq. 3.3 also suggests a prolate
will be aligned to $\bm{A}$ along its short axis if $\bm{g}$ is shifted along
that direction. We discuss the implications in Sec. 5.

We can also derive the angle, $\theta_{g}$, that $\bm{g}$ forms with $\bm{A}$
which is shown in Fig. 8b. Evidently, $\theta_{g}$ does not entirely follow
$\bm{A}$, but is modified by the shape of the grain. Nevertheless, it is
within $\sim 30^{\circ}$ for the range of $\bar{a}$ explored here.

Figure 8:  Panel a: the angle of alignment $\theta_{\text{align}}$ as a
function of $\psi$ and for different values of $\bar{a}$ (a prolate means
$\bar{a}<1$, while an oblate means $\bar{a}>1$). Panel b: the angle from
$\bm{g}$ when at alignment.

##  4 Alignment in a Disk

Sec. 2 and 3 have established how grains are aligned and how the alignment
timescales depend on environmental factors. This section will now implement
the axisymmetric spheroid model in a simple disk environment to assess the
alignment timescales. Observations of polarization permit either prolate
grains aligned with the long axis toroidally around the disk or oblate grains
also aligned with the long axis toroidally around the disk. For brevity, we
will only apply the results from the strictly axisymmetric prolate case. We
first describe how the gas and dust velocities depend on the properties of the
disk and grains in Sec. 4.1. In Sec. 4.2, we implement the velocity profiles
and analyze the alignment timescales across the disk.

###  4.1 Disk Velocity

We assume that the dynamics of gas are not affected by the dust. We use
($R,\Phi$) to denote the cylindrical radius and azimuth of the disk. Let the
Keplerian frequency in the midplane be $\Omega_{K}\equiv\sqrt{GM_{*}/R^{3}}$
where $M_{*}$ is the stellar mass. It is also convenient to define the
Keplerian velocity $v_{K}\equiv\Omega_{K}R$. We define a pressure scale height
as $H\equiv c_{s}/\Omega_{K}$ and let $P_{g}$ be the gas pressure. The
azimuthal velocity of the gas in the midplane follows from the radial force
balance:

$\displaystyle v_{g,\Phi}=v_{K}\sqrt{1-\eta}$ (68)

where the dimensionless factor

$\displaystyle\eta\equiv-\frac{1}{\Omega_{K}^{2}\rho_{g}R}\frac{\partial
P_{g}}{\partial R}=-\bigg{(}\frac{H}{R}\bigg{)}^{2}\frac{\partial\ln
P_{g}}{\partial\ln R}$

represents the level of deviation from $v_{K}$ due to the pressure gradient
(Nakagawa et al., 1986; Armitage, 2015). It is convenient to define the
pressure gradient as $\beta\equiv\partial\ln P_{g}/\partial\ln R$. We assume
the gaseous disk is accreting and parameterize an inward flow by

$\displaystyle v_{g,R}=-\alpha\bigg{(}\frac{H}{R}\bigg{)}^{2}v_{K}$ (69)

where $\alpha$ is a dimensionless parameter555Note that the $\alpha$ here is
simply a parameterization for the inward flow since we do not presume the
origin. Thus, it is not immediately related to the typical $\alpha$-viscosity
prescription from Shakura & Sunyaev (1973). Nevertheless, if viscosity serves
as the origin of the inward flow, then our parameter $\alpha$ should only
differ from the $\alpha$-viscosity by a factor of 1.5. .

The azimuthal velocity of the dust is (Armitage, 2015):

$\displaystyle v_{d,\Phi}=v_{g,\Phi}-\frac{1}{2}\text{St}v_{d,R}$ (70)

where $v_{d,R}$ is the radial velocity of the dust and St is the Stokes number
of the dust. The radial velocity of the dust is

$\displaystyle v_{d,R}=\frac{v_{g,R}-\text{St}\eta v_{K}}{\text{St}^{2}+1}.$
(71)

With these expressions, we can derive the aerodynamic flow of gas experienced
by the dust $\bm{A}$ through:

$\displaystyle A_{R}$ $\displaystyle\equiv
v_{g,R}-v_{d,R}=-\frac{\alpha\text{St}+\beta}{\text{St}+\text{St}^{-1}}\bigg{(}\frac{H}{R}\bigg{)}^{2}v_{K}$
(72) $\displaystyle A_{\Phi}$ $\displaystyle\equiv
v_{g,\Phi}-v_{d,\Phi}=\frac{1}{2}\frac{-\alpha+\beta\text{St}}{\text{St}+\text{St}^{-1}}\bigg{(}\frac{H}{R}\bigg{)}^{2}v_{K}$
(73)

in the radial and azimuthal directions, respectively. The flow speed is simply
$A=\sqrt{A_{R}^{2}+A_{\Phi}^{2}}$. We can also obtain the flow direction
$\chi$ which determines the alignment direction of the grain:

$\displaystyle\chi=\text{atan2}\bigg{[}\frac{1}{2}\frac{-\alpha+\beta\text{St}}{-(\alpha\text{St}+\beta)}\bigg{]}$
(74)

where $\chi=0$ means the dust feels an outward radial gas flow, while
$\chi=90^{\circ}$ means the dust feels an azimuthally directed flow along the
direction of disk rotation (i.e., a tailwind if the dust orbits azimuthally).
Since prolate grains are aligned along the long axis, the polarization angle
will be parallel to $\chi$, i.e., the polarization angle has a $180^{\circ}$
degeneracy and only ranges from $0^{\circ}$ to $180^{\circ}$.

Fig.9 shows the gas flow direction as a function of St for different levels of
$\alpha$. When $\alpha=0$ (Fig. 9, first column), the gas does not accrete
($v_{g,R}=0$), while $v_{g,\Phi}$ is always sub-Keplerian due to the pressure
gradient $\beta$ from Eq. 68 (Fig. 9a, d). The velocity of the dust depends on
St. When $\text{St}\ll 1$, the dust is completely coupled to the gas and the
dust velocity equals the gas velocity. The gas flow, $\bm{A}$, is thus close
to $\bm{0}$ (Fig. 9g). As St increases to $\sim 1$, $v_{d,\Phi}$ begins to
increase (becomes less negative) and $v_{d,R}$ drops quickly. This is the
well-known scenario where dust grains move inwards quickly when $\text{St}\sim
1$ (i.e., the radial drift; Armitage 2015). The gas flow is predominantly in
the radial direction making $\chi\sim 0^{\circ}$ for relatively small
particles with $\text{St}\ll 1$ (Fig. 9j). When $\text{St}\gg 1$, the grains
decouple from the gas and maintain a Keplerian orbit ($v_{d,\Phi}=v_{K}$,
$v_{d,R}=0$). The gas flow becomes opposite to the azimuthal direction making
$\chi\sim-90^{\circ}$ and the dust feels an azimuthal headwind (Fig. 9j).
Thus, for prolate grains aligned to the gas flow along the long axis, the
polarization direction should be parallel to the radial direction when
$\text{St}\ll 1$ and parallel to the azimuthal direction when $\text{St}\gg
1$.

Figure 9:  Radial and azimuthal velocity profiles as a function of St for
different levels of gas accretion, parameterized by $\alpha$. The left to
right columns correspond to $\alpha=0$ (no gas accretion), 0.01, and 0.1. The
adopted pressure gradient, $\beta$, is -2.75. Top row: the azimuthal velocity
subtracted by the Keplerian velocity $v_{K}$ normalized by $(H/R)^{2}v_{K}$.
The black and green solid lines are the profiles for the dust and gas,
respectively. The horizontal dotted line corresponds to $v_{\Phi}$ with
Keplerian rotation. Second row: the radial velocity for the dust and gas
(black and green lines, respectively). The horizontal dotted line corresponds
to $v_{R}=0$. Third row: the gas flow which is the gas velocity relative to
the dust velocity. The radial and azimuthal components of the gas flow,
$A_{R}$ and $A_{\Phi}$, are shown in solid and dashed lines, respectively. The
horizontal dotted line corresponds to no relative velocity. Fourth row: the
gas flow direction $\chi$ in degrees. $\chi=0^{\circ}$ means radial, outward
flow, while $\chi=90^{\circ}$ means azimuthal flow in the rotation direction.
The horizontal dotted lines mark $0^{\circ}$ and $-90^{\circ}$ for visual
guidance.

When $\alpha$ increases (Fig. 9, second and third column), $v_{g,R}$ becomes
non-zero and travels inward (Fig. 9e, f). However, $v_{g,\Phi}$ does not
change, since it remains completely determined by the pressure gradient (Fig.
9b, c). The main effect on the gas flow is when $\text{St}\gg 1$ where the
dust becomes Keplerian with $v_{d,R}\rightarrow 0$. With a non-zero $v_{g,R}$,
the dust eventually matches $v_{g,R}$ as the St increases (Fig. 9e, f) and
$v_{d,R}$ approaches zero from below, which leads to a characteristic Stokes
number $\text{St}_{c}$ in the high St regime where $A_{R}=0$ exactly:

$\displaystyle\text{St}_{c}=-\frac{\beta}{\alpha}.$ (75)

The value of $\text{St}_{c}$ is formally $\infty$ when $\alpha=0$ and of order
$10^{2}$ when $\alpha=0.01$. It becomes $\sim 10$ when $\alpha$ increases to
0.1. For grains with $\text{St}>\text{St}_{c}$, the alignment angle $\chi$
becomes more negative than $-90^{\circ}$, although the deviation remains
relatively small, so that the alignment is approximately azimuthal.

The exploration shows that grains with large St can produce the azimuthal
direction of polarization, and increasing $\alpha$ lessens the strict need for
very large St. At face value, this may serve as evidence of aerodynamically
large grains that are fairly decoupled from the gas. However, even with
$\alpha=0.1$, which would require a fairly large turbulence (Rosotti, 2023) or
magnetic-wind driven accretion (e.g. Suriano et al., 2018), the required St is
$\sim 20$ which is larger than what is typically considered.

The issue may depend on the assumed pressure gradient. Dust rings are fairly
common (e.g. Andrews et al., 2018) and they can serve as traps for large
grains (e.g. Pinilla et al., 2012; Dullemond et al., 2018). Gas kinematics
also infer pressure bumps (Teague et al., 2018). With quickly varying pressure
gradients along the radius across rings, $\beta$ may be closer to 0, for
example, $\beta=0$ at pressure maxima and minima. Given how common dust rings
are, we also explore the case when $\beta=0$ and show that small St can also
produce gas flow parallel to the azimuthal direction.

Figure 10:  The velocity profiles when $\beta=0$ and $\alpha=0.01$ plotted in
the same way as Fig. 9.

In the limit of $\beta=0$, we have
$\chi\rightarrow\text{atan2}[(-1)/(-2\text{St})]$ as long as $\alpha$ is non-
zero. When $\text{St}\ll 1$, $\chi\rightarrow-90^{\circ}$, which would produce
polarization parallel to the azimuthal direction. Fig. 10d shows $\chi$ in the
limit of $\beta=0$. The reason is the following. With $\beta=0$, the azimuthal
velocity of the gas is in Keplerian orbit, i.e., $v_{g,\Phi}=v_{K}$ (Fig.
10a). As the grains accrete with the gas from outer radii, they become super-
Keplerian and feel a headwind towards $\chi=-90^{\circ}$ (Fig. 10c). In the
case of large grains, $\text{St}\gg 1$, they reach purely Keplerian orbits and
only feel the flow from the radial direction due to the accreting gas.
Although the flow speed is smaller by 1 order of magnitude compared to Fig.
9h, it only lengthens the oscillation time $t_{o}$ and does not alter the
damping time $t_{d}$. This limiting case highlights the fact that the
direction of the gas flow can be complicated, especially in highly structured
gas disks, and requires further investigation.

###  4.2 Alignment timescale

We make a simple prescription of a disk applying the velocity field formulated
above. We assume a stellar mass of $M_{*}=0.5$ $M_{\odot}$. The surface
density is a power-law truncated at an inner radius $R_{a}$ and an outer
radius $R_{b}$:

$\displaystyle\Sigma(R)=\Sigma_{1\text{au}}\bigg{(}\frac{R}{1\text{au}}\bigg{)}^{-p}$
(76)

where $\Sigma_{1\text{au}}$ is the surface density at 1 au. We adopt
$R_{a}=0.1$ au, $R_{b}=100$ au, and $p=1$. We assume a total disk mass that is
$5\%$ of the stellar mass, $M_{\text{disk}}=0.05$ $M_{*}$, which gives
$\Sigma_{1\text{au}}\sim 350$ g cm-2.

The temperature of the disk follows

$\displaystyle T(R)=200\bigg{(}\frac{R}{\text{au}}\bigg{)}^{-0.5}.$ (77)

The pressure scale height is $H\equiv c_{s}/\Omega_{K}$. The gas density in
the midplane is

$\displaystyle\rho_{g}(R)\equiv\frac{\Sigma}{\sqrt{2\pi}H}$ (78)

With this setup, the pressure gradient is $\beta=-2.75$. We adopt
$\alpha=0.01$. For the prolate grain, we adopt $c=1$ mm, $\bar{a}=0.9$, and
$g=0.01$.

Fig. 11c shows the timescales as a function of radius in the disk. We find
that the damping time $t_{d}$ is around 4 days ($\sim 10^{-2}$ year) at $R=1$
au, and reaches a thousand years at $R=100$ au. Also, $t_{d}$ is mostly
smaller than the orbital time $t_{K}$ at all radii up to 70 au, meaning grain
oscillations are mostly damped within a fraction of an orbit. The increase in
$t_{d}$ with increasing radius is due to the decrease in both the midplane
density and temperature (Eq. 3.2). Given that $t_{d}$ is comparable to or less
than $t_{K}$, we reason that the alignment mechanism can operate in disks even
for a grain with a small offset between its center of mass and its geometric
center of $g=0.01$.

The oscillation time $t_{o}$ is around a few minutes ($\sim 10^{-5}$ year) at
$R=1$ au and increases to an hour ($\sim 10^{-4}$ year) at $R=100$ au.
Different from $t_{d}$, the oscillation time $t_{o}$ does not increase
strongly with radius, because it is only proportional inversely to the square
root of $\rho_{g}$ and in addition, benefits from an increased level of gas
flow at larger radii. From the adopted parameters, the St is $\sim 10^{-3}$ at
1 au and increases with increasing radius to $\sim 0.1$ because the midplane
density decreases while the grain size is fixed (Fig. 11a). The increasing St
allows the grain to better decouple and experience stronger flow speed $A$
(Fig. 11b).

As discussed previously for the disk and grain parameters chosen in this
illustrative example, the direction of the flow is predominantly in the radial
direction in the simple power-law prescription of a disk (Fig. 11b). More
comprehensive calculations of the gas and dust velocity in structured disks
are necessary to address the question of the alignment direction.

Figure 11:  The resulting power-law profiles of a disk. Panel a: the St. Panel
b: the dashed green and orange lines are the absolute values of the radial and
azimuthal components of the gas flow, respectively. The black line is the
magnitude of the gas flow $A$. Panel c: the solid line is the damping time
$t_{d}$, the dashed line is the oscillation time $t_{o}$, and the dotted line
is the Keplerian orbital time $t_{K}$.

##  5 Discussion

###  5.1 Comparison to other mechanical alignment mechanisms

The badminton birdie-like alignment mechanism utilizes the grains’
interactions with the ambient gas to achieve alignment, which makes it a type
of mechanical (or aerodynamic) alignment. A popular type of mechanical
alignment is alignment due to helicity, which is the ability to spin up a
grain in the presence of gas flow (Lazarian & Hoang, 2007b). We showed that
even without helicity, a grain can be aligned along the gas flow with the
birdie-like alignment mechanism. These are two completely different effects.
Our general torque equation, Eq. 37, likely captures the helicity effect
through the diagonal terms of M. While those terms are 0 for spheroids (Eq.
3.2), if they are non-zero, the grains will obtain a torque in the same
direction of $\bm{A}$ in grain body frame which corresponds to the definition
of helicity. We leave a detailed exploration of the case of non-zero helicity
to a future study.

Eq. 37 also expects a driving torque from specular reflection if $\bm{K}$ is
non-zero. Recall that for spheroids, $\bm{K}=\bm{0}$. However, if $\bm{K}$
exists (depending on the grain geometry), then gas pressure can generate a
torque in the grain frame and spin up the grain. The torque does not require a
gas flow, making it different from the flow-induced restoring torque from the
birdie-like alignment process. This is similar to the systematic torques
proposed by Purcell (1979), but they differ in physical origin. Interestingly,
the driving torque goes as $\propto c_{s}^{2}$ which is much larger than the
restoring torque (related to M) which is $\propto v_{\text{th}}A$. $\bm{K}$
would have to be very small to not be effective. If $\bm{K}$ can increase the
rotational energy of the grain beyond the escape energy ($2P_{s}$; see Sec. 2
and 3), birdie-like alignment effect may be overwhelmed by the rapid spin.
Nevertheless, since $\bm{K}$ is determined entirely from the surface features
of a grain which we expect to be random in general, we speculate that even if
$\bm{K}$ is the dominant component for certain grains, the alignment direction
should also be random as an ensemble and produce a net zero polarization. We
leave a more comprehensive exploration of the non-zero $\bm{K}$ case to a
future study.

Another type of mechanical alignment is the Gold alignment mechanism (Gold,
1952a, b). We can describe the polarization from the Gold alignment mechanism
through two main parts: the first being the actual Gold alignment part and the
second being the polarization from the projection of grains once they are
Gold-aligned. The first part utilizes prolate grains under the presence of a
flow of gas particles. Since the gas particles predominantly impact the grains
along the long axis (the other two axes are negligible), the angular momentum
is entirely parallel to the short axis, and the angular momentum is mostly
confined to a plane perpendicular to the gas flow. As a collection, the
angular momentum directions are random around the direction of the flow and,
for any instance, the directions of the long axes of the grains are random
around their own angular momentum. The second part is the polarization that is
a natural result of projection of the collection of prolate grains with this
configuration. Polarization mostly comes from the grains whose angular momenta
are perpendicular to the observer as long as the flow is not directed to the
observer.

Our formulation of Eq. 37 does not immediately produce the results of Gold
alignment. The grain structure from Gold (1952a) corresponds to a prolate
which has $\bm{K}=\bm{0}$ and also the center of mass is the same as the
geometric center, making $\textbf{{M}}=\textbf{{0}}$. As a result, the torque
from $\bm{A}$ is always $\bm{0}$ regardless of how large $\bm{A}$ is.
Therefore, we would expect that spinning grains damped by the gas should
eventually stop at an orientation that is entirely determined by the initial
conditions (see Eq. 17). One may notice that from Eq. 17 the spin axis will be
around $\bm{e}_{2}$, which is defined to be perpendicular to $\bm{A}$, and we
may obtain polarization from the Gold-projection effect (second part).
However, the polarization is an artificial result because we have assumed that
grains can only oscillate around $\bm{e}_{2}$. In 3D, with
$\textbf{{M}}=\textbf{{0}}$, the orientation of angular momenta and the
resulting final orientation should be entirely random with respect to
$\bm{A}$.

One possible reconciliation is that we have assumed that the speed of each
patch of the surface is small compared to the sound speed of gas (Appendix B).
Numerous studies have found that the Gold alignment requires supersonic drifts
(Gold, 1952a; Purcell, 1969; Lazarian, 1994). Thus, it appears consistent with
Gold alignment that Eq. 37 should not produce a torque for any $\bm{A}$
without the offset between the center of mass and geometric center. We suspect
that Gold alignment corresponds to the supersonic version of Eq. 37, but leave
the verification to a future investigation.

###  5.2 Does birdie-like alignment work in the ISM or protostellar envelopes?

We estimate the damping and oscillation timescales for a grain in ISM
conditions. Using Eq. 3.2, the damping time is

$\displaystyle t_{d,s}\sim 3.5\times
10^{7}\text{years}\bigg{(}\frac{\rho_{s}}{3\text{g
cm}^{-3}}\bigg{)}\bigg{(}\frac{c}{0.1\mu\text{m}}\bigg{)}\bigg{(}\frac{20\text{cm}^{-3}}{n_{g}}\bigg{)}\bigg{(}\frac{0.3\text{km
/s}}{v_{\text{th}}}\bigg{)}$ (79)

when adopting $\bar{a}=0.9$ and $\bar{g}=0.01$; again, $c$ is the grain size,
which we normalize by the classic size of 0.1 $\mu$m for the ISM. Using Eq.
3.2, the oscillation time is

$\displaystyle t_{o,s}\sim 110\text{hours}\sqrt{\frac{\rho_{s}}{3\text{g
cm}^{-3}}}\bigg{(}\frac{c}{0.1\mu\text{m}}\bigg{)}\sqrt{\frac{20\text{cm}^{-3}}{n_{g}}}\sqrt{\frac{0.3\text{km/s}}{v_{\text{th}}}}\sqrt{\frac{1\text{cm/s}}{A}}$

also using $\bar{a}=0.9$ and $\bar{g}=0.01$. We can easily see that the low-
density environment means the grains will oscillate rapidly, while it takes a
long time to damp out the oscillations.

For the ISM, however, there can be other processes that can systematically
spin up the grain to rates comparable to thermal rotation, like radiative
torques (Draine & Weingartner, 1996; Hoang & Lazarian, 2009). We can
parameterize the rotational energy from other spin-up processes as $fkT$ where
$f$ is a multiplication factor and $kT$ is the thermal energy. If the grains
are in thermal equilibrium, we would expect $f=3/2$. We provide an estimation
of how much flow is necessary to trap the grain into oscillation using
$2P_{s}>fkT$. The relation gives a threshold flow speed $A_{t,s}$:

$\displaystyle A>A_{t,s}$
$\displaystyle\equiv\frac{fv_{\text{th}}}{n_{g}c^{3}}\hat{A}_{t,s}\text{ ,}$
$\displaystyle\hat{A}_{t,s}$
$\displaystyle\equiv\frac{1}{16\bar{a}\bar{g}E[1-\bar{z}^{2}]}$ (80)

where $\hat{A}_{t,s}$ is a dimensionless factor and the physical combination
is the characteristic threshold drift speed. For $\bar{a}=0.9$ and
$\bar{g}=0.01$, we get $\hat{A}_{t,d}\sim 5$.

Adopting the same parameters used for Eq. 79, we get an impractical
$A_{t,d}=8f\times 10^{18}$ cm/s. In contrast, $A_{t,d}\sim 0.5$ cm/s for a
disk using the conditions from Eq. 2.3 with the same $\bar{a}$ and $\bar{g}$.
The main issue is that the threshold flow speed is very sensitive to the grain
size $c$, with $A_{t,s}\propto c^{-3}$, meaning it is much harder to trap
small grains. The secondary issue is that the density is much lower in the ISM
than in protoplanetary disks. Thus, the birdie-like mechanism is unlikely to
compete with other mechanisms that can provide thermal rotational energy to
grains in the ISM.

We next assess the dense cores of molecular clouds, since dust polarization is
routinely measured in dense cores and protostellar envelopes (e.g. Stephens et
al., 2013; Maury et al., 2018; Galametz et al., 2018; Le Gouellec et al.,
2019; Yen et al., 2020; Pattle et al., 2021; Cox et al., 2022; Le Gouellec et
al., 2023; Lin et al., 2024b; Huang et al., 2024). We find that the damping
time is

$\displaystyle t_{d,s}\sim 7\times
10^{4}\text{years}\bigg{(}\frac{\rho_{s}}{3\text{g
cm}^{-3}}\bigg{)}\bigg{(}\frac{c}{1\mu\text{m}}\bigg{)}\bigg{(}\frac{10^{5}\text{cm}^{-3}}{n_{g}}\bigg{)}\bigg{(}\frac{0.3\text{km/s}}{v_{\text{th}}}\bigg{)}$
(81)

where we adopted $c=1$ $\mu$m to account for some grain growth from classic
ISM values of 0.1 $\mu$m (e.g. Pagani et al., 2010) and we assumed
$\bar{a}=0.9$ and $\bar{g}=0.01$. The timescale appears less than or
comparable to the lifetimes of the Class 0/I sources embedded in envelopes
(Evans et al., 2009; Williams & Cieza, 2011). Applying the same conditions,
the oscillation time is

$\displaystyle t_{o,s}\sim 16\text{hours}\sqrt{\frac{\rho_{s}}{3\text{g
cm}^{-3}}}\bigg{(}\frac{c}{1\mu\text{m}}\bigg{)}\sqrt{\frac{10^{5}\text{cm}^{-3}}{n_{g}}}\sqrt{\frac{0.3\text{km/s}}{v_{\text{th}}}}\sqrt{\frac{1\text{cm/s}}{A}}.$

The value for $A$ is unclear, but it can be much less than the adopted value
if the small grains are well coupled to the gas. Following Eq. 5.2, we find
that $A_{t,s}\sim 1.5f\times 10^{12}$ cm/s. The large threshold suggests that
birdie-like alignment is unlikely to trap grains in oscillation in dense cores
of molecular clouds or protostellar envelopes if other processes can drive
thermal rotation unless the grains in such regions have already grown much
larger than 1 $\mu$m (as indicated by multi-frequency dust continuum
observations in some cases, e.g., Kwon et al. 2009)

###  5.3 Connection to observations

The drag force from gas plays an important role in the disk dynamics of dust
and is a key ingredient to explain dust rings due to pressure bumps (e.g.
Pinilla et al., 2012; Dullemond et al., 2018), the streaming instability (e.g.
Youdin & Goodman, 2005; Squire & Hopkins, 2020), or pebble accretion (Ormel,
2017). The birdie-like alignment mechanism makes continuum polarization
observations direct evidence of gas drag and offers a possibility to probe
parts of the dust kinematics. While the Doppler shifts of molecular line
emission can trace the gas kinematics, empirical constraints on the dust
kinematics are harder to come by. Our work adds to the idea that continuum
polarization can measure the relative velocities between the gas and dust
(Kataoka et al., 2019; Mori et al., 2019; Tang et al., 2023), but we offer a
new physical justification that grains can be aligned along the long axis even
with highly subsonic drift.

Several disks show polarization that is elliptical in the polarization
orientation with HL Tau as the most prominent example seen from the
multiwavelength polarization (Stephens et al., 2017; Lin et al., 2024a) and
the high angular resolution polarization (Stephens et al., 2023). If
polarization is indeed due to birdie-aligned grains, the $\bm{A}$-field should
be in the azimuthal direction with a $180^{\circ}$ degeneracy.

Intriguingly, there are two disks with a spiral pattern, namely AS 209 (Mori
et al., 2019) and GG Tau (Tang et al., 2023) that are relevant to our
discussion. The polarization patterns for both are nearly azimuthal, but with
a slight deviation by $4.5^{\circ}$ and $7.1^{\circ}$, respectively. The
direction of deviation for both follows an inward spiral. For both sources,
the conventional alignment of grains with helicity cannot easily explain the
inward spiral and they favor a scenario where prolate grains drift inwards and
in the azimuthal direction, which is better explained by the Gold mechanism
(Mori et al., 2019; Tang et al., 2023). The immediate benefit from the birdie-
like alignment is that while the Gold mechanism requires supersonic drift, the
birdie-like alignment only requires subsonic drift and reproduces the same
orientation that the Gold mechanism offers. Whether or not the birdie-like
alignment can produce the inward spiral depends on the proper treatment of the
gas and dust velocity in a disk, or around a ring for these two sources in
particular, which is beyond the scope of this paper.

###  5.4 Caveats and Future Developments

There are several assumptions made in this exploratory work. First, we assume
that angular momentum is predominantly around $\bm{e}_{2}$ (the direction of
the gas flow-induced torque) when, generally, we might expect angular momentum
around the other two axes as well. In the general case, the rotating motion
around $e_{2}$ (and $b_{2}$ considered in this paper) will become precession
in 3D. For example, if the grain contains significant angular momentum around
$\bm{b}_{3}$, $\Gamma_{2}$ will induce precession around $\bm{e}_{3}$.
However, we can expect that the damping torque should inevitably diminish the
spin and precession at timescales of order $t_{d}$ (which is rather short in
the high-density disk environment where the birdie alignment is most
applicable) and reach alignment akin to what is explored in this work. Note
that a spheroid will not exhibit drag when spinning around its axis of
asymmetry, but it should exist for realistic, irregular grains. Another side
effect of the assumption is that grains may artificially exhibit polarization
simply from the Gold-projection effect in the unlikely limiting case of
$\textbf{{M}}=\textbf{{0}}$ as described in Sec. 5.1.

Another assumption is that the level of $\bm{A}$ is fixed as the grain
oscillates. However, the geometric cross-section in the direction of $\bm{A}$
varies as the grain oscillates. Since we know the force on each patch of the
surface, it is a simple exercise to derive the total force on the axisymmetric
spheroid using Eq. 35, Eq. 36, and Eq. 45:

$\displaystyle\bm{F}=\oint_{S}d\bm{F}=-\rho_{g}v_{\text{th}}\oint_{S}\textbf{{N}}d\sigma(\bm{\omega}\times\bm{g}-\bm{A}).$
(82)

For a spheroid, we know $\oint_{S}\textbf{{N}}d\sigma$ from Eq. 50. For the
case of a sphere, $\bm{g}=\bm{0}$ and $\bar{a}=1$, we recover the Epstein drag
of a sphere (Eq. 3; note that $\bm{u}=-\bm{A}$ for that equation). We can
identify in Eq. 82 that $\bm{\omega}\times\bm{g}-\bm{A}$ is the velocity of
the center of the spheroid relative to the gas as it oscillates around the
center of mass. The relation means that the force that determines the
translational motion of the grain can also depend on the rotational motion and
a fully self-consistent study of the velocity of the grain in a disk requires
incorporating the equation of rotational motion. Nevertheless, since we are
particularly interested in the scenario when $\bm{\omega}$ is small (near
alignment) and we found that $\bm{g}$ only needs to be $\sim 1\%$ of the
length scale for alignment, we do not expect the flow to vary much with
oscillation.

The last assumption of importance is the grain structure. We utilized two
cases of an axisymmetric spheroid: the strictly axisymmetric case where the
offset between the geometric and mass centers is along the axis of symmetry
and the quasi-axisymmetric case where the offset can be anywhere in the
spheroid. The assumption allowed us to study alignment analytically. To
explain the disk polarization observations, the center of mass should be
shifted along the long axis of the prolate or oblate instead of the short axis
even though birdie-like alignment does permit prolates or oblates to be
aligned with $\bm{A}$ along the short axes. In other words, observations
suggest that the predominant grain structure favors the offset vector $\bm{g}$
to be along the long axes. It is unclear how grain growth can lead to this
structure, but we can view this as a result that the grain growth mechanism
should strive to explain. Our formulation of Eq. 37 allows one to precalculate
the dynamical properties of the grain ($\bm{K}$, L, and M) independent of
environmental properties or dynamical state. Eq. 37 should now allow a direct
connection between grain structure and its alignment permitting studies of
several geometries quickly. This is a very attractive possibility as we can
study the alignment of a population of grains each with a unique geometry,
allowing the prediction of the degree of alignment and thus formulate testable
predictions of polarization for various grain growth mechanisms. We leave a
more detailed exploration of these assumptions to future studies.

##  6 Conclusions

Disk scale polarization has shown evidence of effectively prolate grains
aligned azimuthally. However, the alignment mechanism is unclear. In this
paper, we demonstrate the possibility that grains can be aligned by the gas
flow when the center of mass is offset from its geometric center, as in the
case of a badminton birdie. Our main results are summarized as follows:

  1. 1.

To build our physical intuition, we first utilize a simple model of a grain
that is composed of two spheres connected by a rigid, massless rod. We show
that the grain behaves as a damped oscillator under the presence of a systemic
flow of gas, $\bm{A}$, if the two spheres of the grain are not identical. The
oscillation is due to restoring torques generated by $\bm{A}$ and the
asymmetry between the spheres, i.e., the “flow-induced restoring torques”
which resists the angular displacement from the direction of alignment. The
oscillation can be described by a potential well near the direction of
alignment. As the gas damps the grain oscillation, the grain reaches alignment
such that the geometric center follows the direction of the flow with respect
to its center of mass. We derive the damping time and oscillation time and
show that the oscillation time is much shorter than the damping time for
typical circumstellar disk conditions (i.e., underdamped oscillations). The
damping time scale is set by the dust stopping time and the dust properties,
particularly the degrees of grain elongation and asymmetry.

  2. 2.

We derive a formula for the torque of a smooth body by considering subsonic
relative motion and specular reflection of gas on each patch of the surface
(Eq. 37). We capture the surface properties of the grains relevant for the
grain dynamics by $\bm{K}$, L, and M which characterizes the driving torque
from gas pressure, the damping torque as the grain spins, and the torque that
arises from the existence of $\bm{A}$, respectively. These quantities allow us
to quantify the grain rotational motion without predetermining environmental
properties, like the gas density or temperature, and the dynamic state, like
the orientation or angular velocity.

  3. 3.

We apply the torque equation to an axisymmetric spheroid (prolate and oblate)
and obtain the damped oscillator equation of motion along with analytical
solutions to the damping time and oscillation time. By considering a grain
whose center of mass is shifted along the axis of symmetry, which is
characterized by $\bm{g}$ (the offset vector from the center of mass to the
center of the spheroid), we find that the final alignment direction is such
that $\bm{g}$ is parallel to $\bm{A}$. In other words, for a prolate, the long
axis is along the gas flow, or for an oblate, the short axis is along the gas
flow.

  4. 4.

We further consider a grain whose center of mass can be anywhere within the
spheroidal surface and also derive analytical solutions to the damping time
and oscillation time. We show that the final alignment direction depends on
both $\bm{g}$ and the shape of the spheroid. In particular, if the center of
mass of an oblate is shifted along the long axis, then the oblate will be
aligned along the long axis. Since disk polarization observations require
effectively prolate grains, the birdie-like alignment mechanism requires the
offset vector $\bm{g}$ to be along the long axis for both prolates and
oblates.

  5. 5.

We implement a simple disk model and calculate the gas and dust velocity
fields. Using the power-law disk model with a pressure gradient, we expect the
gas flow to be in the radial direction for small Stokes number, St, and in the
azimuthal direction for large St. However, in the center of pressure bumps or
gaps, where the pressure gradient can be small, the gas flow can be in the
azimuthal direction for small St and in the radial direction for large St.
Using just the power-law disk model to study the radial dependence of the
mechanism, we show that the damping time is mostly less than the Keplerian
time $t_{K}$, but it can be comparable to or greater than $t_{K}$ at large
radii near the disk outer edge. Typically, it takes a fraction of an orbit or
only several orbits to align a grain in a disk. The very short alignment
timescale compared to the lifetime of the disk makes birdie-like alignment a
very promising mechanism to explain the inferred toroidally aligned,
effectively prolate grains from several disk polarization observations and
deserves further development.

## Acknowledgements

Z.-Y.D.L. acknowledges support from NASA 80NSSC18K1095, the Jefferson Scholars
Foundation, the NRAO ALMA Student Observing Support (SOS) SOSPA8-003, the
Achievements Rewards for College Scientists (ARCS) Foundation Washington
Chapter, the Virginia Space Grant Consortium (VSGC), and UVA research
computing (RIVANNA). ZYL is supported in part by NASA 80NSSC20K0533 and NSF
AST-2307199. L.W.L. acknowledges support from NSF AST-1910364 and NSF
AST-2307844. MFL acknowledges support from the European Research Executive
Agency HORIZON-MSCA-2021-SE-01 Research and Innovation programme under the
Marie Skłodowska-Curie grant agreement number 101086388 (LACEGAL). MFL is also
grateful for the warmth of the Star Formation group of the Institut de
Ciències de l’Espai and the Universidad de Barcelona during his 2024 research
stay.

## Data Availability

Data underlying this article is available from the corresponding author upon
request.

## References

  * Alves et al. (2018) Alves F. O., et al., 2018, A&A, 616, A56
  * Andrews et al. (2018) Andrews S. M., et al., 2018, ApJ, 869, L41
  * Armitage (2015) Armitage P. J., 2015, arXiv e-prints, p. arXiv:1509.06382
  * Aso et al. (2021) Aso Y., Kwon W., Hirano N., Ching T.-C., Lai S.-P., Li Z.-Y., Rao R., 2021, ApJ, 920, 71
  * Bacciotti et al. (2018) Bacciotti F., et al., 2018, ApJ, 865, L12
  * Cohen et al. (2015) Cohen C., Darbois Texier B., Quéré D., Clanet C., 2015, New Journal of Physics, 17, 063001
  * Cox et al. (2022) Cox E. G., et al., 2022, ApJ, 932, 34
  * Dent et al. (2019) Dent W. R. F., Pinte C., Cortes P. C., Ménard F., Hales A., Fomalont E., de Gregorio-Monsalvo I., 2019, MNRAS, 482, L29
  * Dolginov (1972) Dolginov A. Z., 1972, Ap&SS, 18, 337
  * Dolginov & Mytrophanov (1976) Dolginov A. Z., Mytrophanov I. G., 1976, Ap&SS, 43, 257
  * Draine & Weingartner (1996) Draine B. T., Weingartner J. C., 1996, ApJ, 470, 551
  * Draine & Weingartner (1997) Draine B. T., Weingartner J. C., 1997, ApJ, 480, 633
  * Drążkowska et al. (2023) Drążkowska J., et al., 2023, in Inutsuka S., Aikawa Y., Muto T., Tomida K., Tamura M., eds, Astronomical Society of the Pacific Conference Series Vol. 534, Protostars and Planets VII. p. 717 (arXiv:2203.09759), doi:10.48550/arXiv.2203.09759
  * Dullemond et al. (2018) Dullemond C. P., et al., 2018, ApJ, 869, L46
  * Epstein (1924) Epstein P. S., 1924, Physical Review, 23, 710
  * Evans et al. (2009) Evans Neal J. I., et al., 2009, ApJS, 181, 321
  * Galametz et al. (2018) Galametz M., et al., 2018, A&A, 616, A139
  * Girart et al. (2006) Girart J. M., Rao R., Marrone D. P., 2006, Science, 313, 812
  * Girart et al. (2018) Girart J. M., et al., 2018, ApJ, 856, L27
  * Gold (1952a) Gold T., 1952a, MNRAS, 112, 215
  * Gold (1952b) Gold T., 1952b, Nature, 169, 322
  * Harrison et al. (2019) Harrison R. E., et al., 2019, ApJ, 877, L2
  * Harrison et al. (2021) Harrison R. E., et al., 2021, ApJ, 908, 141
  * Harrison et al. (2024) Harrison R. E., Lin Z.-Y. D., Looney L. W., Li Z.-Y., Yang H., Stephens I. W., Fernández-López M., 2024, ApJ, 967, 40
  * Hoang & Lazarian (2009) Hoang T., Lazarian A., 2009, ApJ, 697, 1316
  * Hoang et al. (2022) Hoang T., Tram L. N., Minh Phan V. H., Giang N. C., Phuong N. T., Dieu N. D., 2022, AJ, 164, 248
  * Huang et al. (2024) Huang B., et al., 2024, ApJ, 963, L31
  * Kataoka et al. (2015) Kataoka A., et al., 2015, ApJ, 809, 78
  * Kataoka et al. (2016a) Kataoka A., Muto T., Momose M., Tsukagoshi T., Dullemond C. P., 2016a, ApJ, 820, 54
  * Kataoka et al. (2016b) Kataoka A., et al., 2016b, ApJ, 831, L12
  * Kataoka et al. (2017) Kataoka A., Tsukagoshi T., Pohl A., Muto T., Nagai H., Stephens I. W., Tomisaka K., Momose M., 2017, ApJ, 844, L5
  * Kataoka et al. (2019) Kataoka A., Okuzumi S., Tazaki R., 2019, ApJ, 874, L6
  * Kirchschlager & Bertrang (2020) Kirchschlager F., Bertrang G. H. M., 2020, A&A, 638, A116
  * Kirchschlager et al. (2019) Kirchschlager F., Bertrang G. H. M., Flock M., 2019, MNRAS, 488, 1211
  * Kwon et al. (2009) Kwon W., Looney L. W., Mundy L. G., Chiang H.-F., Kemball A. J., 2009, ApJ, 696, 841
  * Kwon et al. (2019) Kwon W., Stephens I. W., Tobin J. J., Looney L. W., Li Z.-Y., van der Tak F. F. S., Crutcher R. M., 2019, ApJ, 879, 25
  * Lazarian (1994) Lazarian A., 1994, MNRAS, 268, 713
  * Lazarian (2007) Lazarian A., 2007, J. Quant. Spectrosc. Radiative Transfer, 106, 225
  * Lazarian & Efroimsky (1999) Lazarian A., Efroimsky M., 1999, MNRAS, 303, 673
  * Lazarian & Hoang (2007a) Lazarian A., Hoang T., 2007a, MNRAS, 378, 910
  * Lazarian & Hoang (2007b) Lazarian A., Hoang T., 2007b, ApJ, 669, L77
  * Lazarian et al. (2015) Lazarian A., Andersson B. G., Hoang T., 2015, in , Polarimetry of Stars and Planetary Systems. p. 81, doi:10.48550/arXiv.1511.03696
  * Le Gouellec et al. (2019) Le Gouellec V. J. M., et al., 2019, ApJ, 885, 106
  * Le Gouellec et al. (2023) Le Gouellec V. J. M., Maury A. J., Hull C. L. H., 2023, A&A, 671, A167
  * Lee et al. (2018) Lee C.-F., Li Z.-Y., Ching T.-C., Lai S.-P., Yang H., 2018, ApJ, 854, 56
  * Lesur et al. (2023) Lesur G., et al., 2023, in Inutsuka S., Aikawa Y., Muto T., Tomida K., Tamura M., eds, Astronomical Society of the Pacific Conference Series Vol. 534, Protostars and Planets VII. p. 465
  * Lin et al. (2020a) Lin Z.-Y. D., Li Z.-Y., Yang H., Looney L., Lee C.-F., Stephens I., Lai S.-P., 2020a, MNRAS, 493, 4868
  * Lin et al. (2020b) Lin Z.-Y. D., Li Z.-Y., Yang H., Looney L., Stephens I., Hull C. L. H., 2020b, MNRAS, 496, 169
  * Lin et al. (2022) Lin Z.-Y. D., Li Z.-Y., Yang H., Stephens I., Looney L., Harrison R., Fernández-López M., 2022, MNRAS, 512, 3922
  * Lin et al. (2024a) Lin Z.-Y. D., et al., 2024a, MNRAS, 528, 843
  * Lin et al. (2024b) Lin S.-J., et al., 2024b, ApJ, 961, 117
  * Liu et al. (2024) Liu Y., Takahashi S., Machida M., Tomisaka K., Girart J. M., Ho P. T. P., Nakanishi K., Sato A., 2024, ApJ, 963, 104
  * Maury et al. (2018) Maury A. J., et al., 2018, MNRAS, 477, 2760
  * Mishchenko & Travis (1994) Mishchenko M. I., Travis L. D., 1994, Optics Communications, 109, 16
  * Mori & Kataoka (2021) Mori T., Kataoka A., 2021, ApJ, 908, 153
  * Mori et al. (2019) Mori T., Kataoka A., Ohashi S., Momose M., Muto T., Nagai H., Tsukagoshi T., 2019, ApJ, 883, 16
  * Nakagawa et al. (1986) Nakagawa Y., Sekiya M., Hayashi C., 1986, Icarus, 67, 375
  * Nguyen Tat et al. (2023) Nguyen Tat T., Diep P. N., Hoang T., Tram L. N., Bich Ngoc N., Phuong N. T., Truong B., 2023, arXiv e-prints, p. arXiv:2401.00220
  * Ohashi et al. (2018) Ohashi S., et al., 2018, ApJ, 864, 81
  * Ohashi et al. (2023) Ohashi S., et al., 2023, ApJ, 954, 110
  * Ormel (2017) Ormel C. W., 2017, in Pessah M., Gressel O., eds, Astrophysics and Space Science Library Vol. 445, Formation, Evolution, and Dynamics of Young Solar Systems. p. 197, doi:10.1007/978-3-319-60609-5_7
  * Pagani et al. (2010) Pagani L., Steinacker J., Bacmann A., Stutz A., Henning T., 2010, Science, 329, 1622
  * Pattle et al. (2021) Pattle K., et al., 2021, ApJ, 907, 88
  * Pinilla et al. (2012) Pinilla P., Birnstiel T., Ricci L., Dullemond C. P., Uribe A. L., Testi L., Natta A., 2012, A&A, 538, A114
  * Purcell (1969) Purcell E. M., 1969, Physica, 41, 100
  * Purcell (1979) Purcell E. M., 1979, ApJ, 231, 404
  * Rosotti (2023) Rosotti G. P., 2023, New Astron. Rev., 96, 101674
  * Sadavoy et al. (2018) Sadavoy S. I., et al., 2018, ApJ, 869, 115
  * Sadavoy et al. (2019) Sadavoy S. I., et al., 2019, ApJS, 245, 2
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Squire & Hopkins (2020) Squire J., Hopkins P. F., 2020, MNRAS, 498, 1239
  * Stephens et al. (2013) Stephens I. W., et al., 2013, ApJ, 769, L15
  * Stephens et al. (2014) Stephens I. W., et al., 2014, Nature, 514, 597
  * Stephens et al. (2017) Stephens I. W., et al., 2017, ApJ, 851, 55
  * Stephens et al. (2020) Stephens I. W., Fernández-López M., Li Z.-Y., Looney L. W., Teague R., 2020, ApJ, 901, 71
  * Stephens et al. (2023) Stephens I. W., et al., 2023, Nature, 623, 705
  * Suriano et al. (2018) Suriano S. S., Li Z.-Y., Krasnopolsky R., Shang H., 2018, MNRAS, 477, 1239
  * Takahashi et al. (2019) Takahashi S., Machida M. N., Tomisaka K., Ho P. T. P., Fomalont E. B., Nakanishi K., Girart J. M., 2019, ApJ, 872, 70
  * Tang et al. (2023) Tang Y.-W., et al., 2023, ApJ, 947, L5
  * Tazaki et al. (2017) Tazaki R., Lazarian A., Nomura H., 2017, ApJ, 839, 56
  * Teague et al. (2018) Teague R., Bae J., Bergin E. A., Birnstiel T., Foreman-Mackey D., 2018, ApJ, 860, L12
  * Waterman (1971) Waterman P. C., 1971, Phys. Rev. D, 3, 825
  * Williams & Cieza (2011) Williams J. P., Cieza L. A., 2011, ARA&A, 49, 67
  * Yang (2021) Yang H., 2021, ApJ, 911, 125
  * Yang & Li (2020) Yang H., Li Z.-Y., 2020, ApJ, 889, 15
  * Yang et al. (2016a) Yang H., Li Z.-Y., Looney L., Stephens I., 2016a, MNRAS, 456, 2794
  * Yang et al. (2016b) Yang H., Li Z.-Y., Looney L. W., Cox E. G., Tobin J., Stephens I. W., Segura-Cox D. M., Harris R. J., 2016b, MNRAS, 460, 4109
  * Yang et al. (2017) Yang H., Li Z.-Y., Looney L. W., Girart J. M., Stephens I. W., 2017, MNRAS, 472, 373
  * Yang et al. (2019) Yang H., Li Z.-Y., Stephens I. W., Kataoka A., Looney L., 2019, MNRAS, 483, 2371
  * Yang et al. (2024) Yang H., Fernández-López M., Li Z.-Y., Stephens I. W., Looney L. W., Lin Z.-Y. D., Harrison R., 2024, ApJ, 963, 134
  * Yen et al. (2020) Yen H.-W., et al., 2020, ApJ, 893, 54
  * Youdin & Goodman (2005) Youdin A. N., Goodman J., 2005, ApJ, 620, 459
  * Zamponi et al. (2024) Zamponi J., Maureira M. J., Liu H. B., Zhao B., Segura-Cox D., Ko C.-L., Caselli P., 2024, A&A, 682, A56
  * Zhang et al. (2023) Zhang S., Zhu Z., Ueda T., Kataoka A., Sierra A., Carrasco-González C., Macías E., 2023, ApJ, 953, 96
  * van de Hulst (1957) van de Hulst H. C., 1957, Light Scattering by Small Particles

##  Appendix A Additional Example of a Phase Portrait

Fig. 4 from Sec. 2 presented a phase portrait showing the damped oscillation
behavior of a double-sphere model of a grain. The typical disk conditions mean
grains undergo highly underdamped oscillations. To illustrate how a better-
damped phase portrait appears, we adopt $n_{g}=10^{12}$ cm-3,
$v_{\text{th}}=0.2$ km/s, $A=5\times 10^{-5}$ cm/s, $l=1$ mm, $\epsilon=1$,
and $\kappa=1.01$ and the results are shown in Fig. 12. With a much higher gas
density and smaller velocity difference, the trajectories now clearly spiral
toward the attractors (final alignment points), even starting from initially
unbound states considered in this plot. For reference, using Eq. 2.3, the
damping ratio $\zeta_{d}\sim 0.14$ in this illustrative case.

Figure 12:  The phase portrait plotted in the same way as Fig. 4. In this
example, the grain is more damped.

##  Appendix B Derivation of the drag force on a surface element

In this appendix, we will use a distribution of velocities of the gas for a
moving observer and derive the transfer of momentum from the gas onto a patch
of surface after integrating over the velocity distribution. The derivation is
based on Epstein (1924). We will redefine the symbols utilized in this
appendix for clarity and the same symbols in the main text may not correspond
to the definition here.

We begin by considering the Maxwell velocity distribution for the gas which
gives the number of molecules having velocity components in the cartesian
coordinates between $\xi_{0}$, $\eta_{0}$, $\zeta_{0}$ and $\xi_{0}+d\xi_{0}$,
$\eta_{0}+d\eta_{0}$, $\zeta_{0}+d\zeta_{0}$ through:

$\displaystyle
N_{\xi_{0},\eta_{0},\zeta_{0}}d\xi_{0}d\eta_{0}d\zeta_{0}=N\bigg{(}\frac{h}{\pi}\bigg{)}^{3/2}e^{-h(\xi_{0}^{2}+\eta_{0}^{2}+\zeta_{0}^{2})}d\xi_{0}d\eta_{0}d\zeta_{0}$
(83)

where $N$ is the number of molecules per unit volume. $h$ is defined by

$\displaystyle h\equiv\frac{m}{2kT}$

where $m$ is the mass of each molecule, $k$ is the Boltzmann constant, and $T$
is the temperature.

Consider a traveling observer with a cartesian coordinate system with axes
$x$, $y$, and $z$. The observer travels relative to the rest frame of the gas
with a speed of $V$ and components $\alpha V$, $\beta V$, and $\gamma V$ along
the $x$, $y$, and $z$ axes, respectively. Here, $\alpha$, $\beta$, and
$\gamma$ are direction cosines (not to be confused with $\alpha$ and $\beta$
from Sec. 4). Let $\xi$, $\eta$, and $\zeta$ be the velocity components of a
molecule parallel to $x$, $y$, and $z$, respectively, seen by the observer.
The transformation between the velocity in the observer frame to the velocity
in the gas rest frame follows

$\displaystyle\xi=\xi_{0}-\alpha V\text{ , }$
$\displaystyle\eta=\eta_{0}-\beta V\text{ , }$
$\displaystyle\zeta=\zeta_{0}-\gamma V\text{ . }$ (84)

Since $V$ is constant, $d\xi=d\xi_{0}$, $d\eta=d\eta_{0}$, and
$d\zeta=d\zeta_{0}$. Thus, the velocity distribution seen by the observer is

$\displaystyle N_{\xi,\eta,\zeta}d\xi d\eta d\zeta=$ $\displaystyle
N\bigg{(}\frac{h}{\pi}\bigg{)}^{3/2}e^{-h[(\xi+\alpha V)^{2}+(\eta+\beta
V)^{2}+(\zeta+\gamma V)^{2}]}d\xi d\eta d\zeta$ (85)

We approximate the equation by considering the distribution to the first power
of $V$ which gives

$\displaystyle
N_{\xi,\eta,\zeta}=N\bigg{(}\frac{h}{\pi}\bigg{)}^{\frac{3}{2}}\bigg{[}1-2hV(\alpha\xi+\beta\eta+\gamma\zeta)\bigg{]}e^{-h(\xi^{2}+\eta^{2}+\zeta^{2})}\text{
.}$ (86)

We now wish to calculate the number of impinging particles onto some surface
element $d\sigma$. Let the $x$-axis be the direction that is normal to the
surface element. Within some unit time $dt$, the impinging particles with
velocity ($\xi$, $\eta$, $\zeta$) will be those that are enclosed within a
cylinder with $d\sigma$ as the base and $-\xi dt$ as the height:

$\displaystyle n_{\xi,\eta,\zeta}d\xi d\eta d\zeta dtd\sigma=$ (87)
$\displaystyle-N\bigg{(}\frac{h}{\pi}\bigg{)}^{\frac{3}{2}}\bigg{[}\xi-2hV\xi(\alpha\xi+\beta\eta+\gamma\zeta)\bigg{]}e^{-h(\xi^{2}+\eta^{2}+\zeta^{2})}d\xi
d\eta d\zeta dtd\sigma\text{ .}$

Each imping particle carries a momentum of

$\displaystyle\bm{p}_{i}=m\begin{pmatrix}\xi\\\ \eta\\\ \zeta\end{pmatrix}$
(88)

in the observer frame. Since we know the number of impinging particles onto
the surface from Eq. 87, we can obtain the momentum from the impacts of all
velocities

$\displaystyle\bm{p}^{(i)}_{\sigma,t}dtd\sigma$
$\displaystyle=\int_{-\infty}^{0}d\xi\int_{-\infty}^{\infty}d\eta\int_{-\infty}^{\infty}d\zeta\bm{p}_{i}n_{\xi,\eta,\zeta}dtd\sigma$
(89) $\displaystyle=-mN\Bigg{[}\begin{pmatrix}\frac{1}{4h}\\\ 0\\\
0\end{pmatrix}+\frac{V}{2\sqrt{\pi h}}\begin{pmatrix}2\alpha\\\ \beta\\\
\gamma\end{pmatrix}\Bigg{]}dtd\sigma$ (90)

where $\bm{p}^{(i)}_{\sigma,t}$ is the momentum per unit area and per unit
time.

To calculate the momentum of the emerging gas, we assume specular reflection
in which, after the impact, the velocity of the molecule normal to the surface
flips by a negative sign. The momentum of the emerging molecule is thus
$\bm{p}_{e}=m(-\xi,\eta,\zeta)$. Following the same exercise as above,

$\displaystyle\bm{p}^{(e)}_{\sigma,t}dtd\sigma$
$\displaystyle=\int_{-\infty}^{0}d\xi\int_{-\infty}^{\infty}d\eta\int_{-\infty}^{\infty}d\zeta\bm{p}_{e}n_{\xi,\eta,\zeta}dtd\sigma$
(91) $\displaystyle=-mN\Bigg{[}\begin{pmatrix}-\frac{1}{4h}\\\ 0\\\
0\end{pmatrix}+\frac{V}{2\sqrt{\pi h}}\begin{pmatrix}-2\alpha\\\ \beta\\\
\gamma\end{pmatrix}\Bigg{]}.$ (92)

We can obtain the force onto the surface by looking at the momentum that is
deposited by gas within the unit time:

$\displaystyle d\bm{F}$
$\displaystyle\equiv\frac{[\bm{p}^{(i)}_{\sigma,t}-\bm{p}^{(e)}_{\sigma,t}]dtd\sigma}{dt}$
(93) $\displaystyle=-mN\bigg{(}\frac{1}{2h}+\frac{2V\alpha}{\sqrt{\pi
h}}\bigg{)}\bm{n}d\sigma$ (94)

where $d\bm{F}$ is the force on $d\sigma$ and $\bm{n}$ is the unit direction
along the $x$-axis, which was defined to be in the normal direction of the
surface. One can quickly realize that the force is only in the normal
direction because the momentum of the molecules along the $y$\- or
$z$-directions has not changed. Also, the direction of the force is anti-
parallel to $\bm{n}$.

We can re-express the equation by connecting it to a few physical quantities.
$1/2h$ is simply the isothermal sound speed squared $c_{s}^{2}\equiv kT/m$.
Also, $2/\sqrt{\pi h}$ is the average speed $v_{\text{th}}\equiv\sqrt{8kT/(\pi
m)}$. The quantity $mN$ is the mass density $\rho_{g}$ as defined in the main
text. Furthermore, since $V\alpha$ is the velocity component along the
$x$-axis, we can express it as $V\alpha=\bm{V}\cdot\bm{n}$. When expressed in
the dot product form and applied to $\bm{n}d\sigma$ for Eq. 93, one can
quickly see that $(\bm{V}\cdot\bm{n})\bm{n}$ is simply a vector from $\bm{V}$
that is projected onto the normal direction. We can express this relation
through a projection tensor $\textbf{{N}}\equiv\bm{n}\bm{n}$ giving
$(\bm{V}\cdot\bm{n})\bm{n}=\textbf{{N}}(\bm{V})$. Implementing these new
relations, we have

$\displaystyle
d\bm{F}=-\rho_{g}\bigg{[}c_{s}^{2}\bm{n}+v_{\text{th}}\textbf{{N}}(\bm{V})\bigg{]}d\sigma$
(95)

which matches Eq. 35. Note that the first term is simply the pressure which
does not depend on $\bm{V}$. The second term controls the effect of the motion
of the surface. When $\bm{V}$ is in the same direction as $\bm{n}$, the side
of the surface that experiences reflection feels an extra force in the
direction opposite to $\bm{n}$. When $\bm{V}$ is in the opposite direction of
$\bm{n}$, the surface feels a reduced force.

##  Appendix C Analytical Solutions to the Integrated Constants

In the main text, Eq. 47 defines an integration form that frequently appears
in the equations of motion in Sec. 3. The actual quantities that appeared are
$E[\bar{z}^{2}]$, $E[1-\bar{z}^{2}]$, and $E[\bar{z}^{2}-\bar{z}^{4}]$. Note
that $E[f(\bar{z})+g(\bar{z})]=E[f(\bar{z})]+E[g(\bar{z})]$ where $f$ and $g$
are some functions of $\bar{z}$. The property means we only need to solve for
$E[1]$, $E[\bar{z}^{2}]$, and $E[\bar{z}^{4}]$. Here we provide the analytical
solutions.

For convenience, we can define the elongation factor for prolates (when
$\bar{a}<1$) as

$\displaystyle e\equiv\sqrt{1-\bar{a}^{2}}$ (96)

and the elongation factor for oblates (when $\bar{a}>1$) as

$\displaystyle\varepsilon\equiv\sqrt{\bar{a}^{2}-1}$ (97)

(not to be confused with $\epsilon$ in Sec. 2). Both are positive quantities.
We find that

$E[1]=\begin{cases}\frac{2\arcsin{(e)}}{e}&\text{ , when }\bar{a}<1\\\
2&\text{ , when }\bar{a}=1\\\
\frac{2\operatorname{arcsinh}{(\varepsilon)}}{\varepsilon}&\text{ , when
}\bar{a}>1\end{cases}$ (98)

where $\operatorname{arcsinh}$ is the inverse hyperbolic sine. Furthermore,

$E[\bar{z}^{2}]=\begin{cases}-\frac{\sqrt{1-e^{2}}}{e^{2}}+\frac{\arcsin{(e)}}{e^{3}}&\text{
, when }\bar{a}<1\\\ \frac{2}{3}&\text{ , when }\bar{a}=1\\\
\frac{\sqrt{1+\varepsilon^{2}}}{\varepsilon^{2}}-\frac{\operatorname{arcsinh}{(\varepsilon)}}{\varepsilon^{3}}&\text{
, when }\bar{a}>1.\\\ \end{cases}$ (99)

Lastly, we can obtain

$E[\bar{z}^{4}]=\begin{cases}-\frac{(2e^{2}+3)\sqrt{1-e^{2}}}{4e^{4}}+\frac{3\arcsin{(e)}}{4e^{5}}&\text{
, when }\bar{a}<1\\\ \frac{2}{5}&\text{ , when }\bar{a}=1\\\
\frac{(2\varepsilon^{2}-3)\sqrt{1+\varepsilon^{2}}}{4\varepsilon^{4}}+\frac{3\operatorname{arcsinh}{(\varepsilon)}}{4\varepsilon^{5}}&\text{
, when }\bar{a}>1.\\\ \end{cases}$ (100)

With Eq. 98, 99, and 100, we can get $E[1-\bar{z}^{2}]=E[1]-E[\bar{z}^{2}]$
and $E[\bar{z}^{2}-\bar{z}^{4}]=E[\bar{z}^{2}]-E[\bar{z}^{4}]$.