# The Most Likely Transition Path
for a Class of Distribution-Dependent Stochastic Systems

Wei Wei1,222Email<EMAIL_ADDRESS>and Jianyu Hu2,111Email:
<EMAIL_ADDRESS>the corresponding author
1Institute of Natural Sciences, Shanghai Jiao Tong University,
Shanghai 200240, China
2Division of Mathematical Sciences, Nanyang Technological University,
21 Nanyang Link, 637371, Singapore.

###### Abstract

Distribution-dependent stochastic dynamical systems arise widely in
engineering and science. We consider a class of such systems which model the
limit behaviors of interacting particles moving in a vector field with random
fluctuations. We aim to examine the most likely transition path between
equilibrium stable states of the vector field. In the small noise regime, the
action functional does not involve the solution of the skeleton equation which
describes the unperturbed deterministic flow of the vector field shifted by
the interaction at zero distance. As a result, we are led to study the most
likely transition path for a stochastic differential equation without
distribution dependency. This enables the computation of the most likely
transition path for these distribution-dependent stochastic dynamical systems
by the adaptive minimum action method and we illustrate our approach in two
examples.

Keywords and Phrases: Large deviations, interacting particles, adaptive
minimum action method, distribution-dependent stochastic dynamics, the most
likely transition path, McKean-Vlasov stochastic systems

##  1 Introduction

Large systems of interacting particles are now quite common in many fields.
Particles can represent ions and electrons in plasmas in the context of
physics [11]. In the bio-science, they can model the motion of cells or
bacteria and in economics or social sciences, they represent “agents” [9, 10].
But they are analytically complex and computationally expensive as the number
of particles is very large. Instead, a corresponding mean-field limit equation
is often an alternative model to investigate. There are rich literatures on
deriving such mean-field limit equations; see, for example, [5, 6, 7].
Interactions in such equations often appear as the distribution in the drift
term, which makes them known as distribution-dependent equations. The McKean-
Vlasov stochastic differential equation was first studied in [8] and is a
widely applicable one among these distribution-dependent equations.

The motion of the particles modeled by the McKean-Vlasov stochastic
differential equation can be simply described as follows: particles
interacting with each other wander under random fluctuations in a landscape,
which is formed by a vector field. We are interested in how such a particle
transfers from one equilibrium stable state of the landscape to another
equilibrium stable state. Large deviations theory is a useful tool to deal
with this kind of problem. Roughly speaking, the large deviations theory
characterizes how small the probability of a rare event is via an action
functional. It also extends the concept of potential to non-gradient systems
as quasi-potential, which is widely used in investigating the long time
behaviors of a system under small random perturbations. More details of the
large deviations theory may be found in [2, 3].

There are many large deviations results investigating the McKean–Vlasov
stochastic differential equations. Herrmann et al. [4] derived a large
deviation result for a one-dimensional McKean-Vlasov stochastic equation.
Based on this, they investigated the influence of interactions on exit times
and exit locations. Recently, Reis et al. [1] obtained large deviations
results for a more general type of $n$-dimensional McKean-Vlasov stochastic
equations with Gaussian noise. Liu et al.[12] extended the large deviations
results to the McKean–Vlasov stochastic equations with jumps. Although these
equations contain the distribution of the solution in the drift term, the
action functionals they obtained are simplified to only including a Dirac
measure. Therefore, numerical methods built on the large deviations theory can
be applied to investigate the dynamical behaviors of such McKean-Vlasov
stochastic equations.

The minimum action method(MAM) is a widely used numerical method in finding
the most likely transition path. It was proposed to compute the most likely
transition path for stochastic differential equations with multiple metastable
states [16]. The basic idea of MAM is to find the minimizer of the associated
action functional and take this minimizer as the desired most likely
transition path. MAM was improved by a reparametrization [13] to a fixed
finite time scale and numerically solving an Euler-Lagrange equation to find
such a minimizer. A moving mesh method was used in [14, 15] to enhance the
performance of MAM. Authors in [17] employed finite elements approximation to
find the minimizer of the associated action functional.

In this paper, we calculate the action functionals for a class of
McKean–Vlasov stochastic differential equations. We find that the action
functional, which characterizes the probability of rare events, usually
involves a ‘modified deterministic flow’, which is the flow of the vector
field shifted by the interaction at zero distance. When we examine the most
likely transition path between equilibrium stable states of the vector field,
we find that the action functional does not involve this modified
deterministic flow. A numerical method based on the adaptive minimum action
method [14, 15] is used to calculate the most likely transition path.

The paper is arranged as follows. In section 2, we briefly introduce the
background of the McKean–Vlasov stochastic equation and the associated large
deviations theory. Based on the large deviations theory, we analyze the
transition behaviors of particles modeled by the McKean–Vlasov stochastic
equation. In section 3, we will introduce the numerical method for computing
the most likely transition path and present the most likely transition paths
for a bi-stable system under two different interactions.

##  2 Backgrounds

We consider the following $d$-dimensional McKean-Vlasov equation with
$X^{\epsilon}_{0}=x_{1}$ ,

$dX_{t}^{\epsilon}=V\left(X_{t}^{\epsilon}\right)dt-F*u_{t}^{\epsilon}(X_{t}^{\epsilon})dt+\sqrt{\epsilon}dB_{t},$
(2.1)

where $u^{\epsilon}_{t}$ stands for the law of $X^{\epsilon}_{t}$ and solution
$X^{\epsilon}$ is a stochastic process on the given probability space
$(\Omega,\mathcal{F},\mathbb{P})$. The symbol “$*$” is the usual convolution
operator that, for $\omega_{0}\in\Omega$

$\displaystyle F*u_{s}^{\epsilon}(X_{s}^{\epsilon}(\omega_{0}))$
$\displaystyle=\int_{\mathbb{R}^{d}}F(X_{s}^{\epsilon}(\omega_{0})-x)du_{s}^{\epsilon}(dx)$
$\displaystyle=\int_{\omega\in\Omega}F\left(X_{s}^{\epsilon}(\omega_{0})-X_{s}^{\epsilon}(\omega)\right)d\mathbb{P}(\omega)=\mathbb{E}\left(F(X_{s}^{\epsilon}(\omega_{0})-X_{s}^{\epsilon})\right).$

Equation (2.1) can be regarded as the mean-field limit of the following system
of many particles as $N\to\infty$, which is known as the _propagation of
chaos_.

Let $X^{\varepsilon}_{i,N}(0)=x_{1}$ and

$dX^{\varepsilon}_{i,N}(t)=V\left(X^{\varepsilon}_{i,N}(t)\right)dt-\frac{1}{N}\sum_{j=1}^{N}F\left(X^{\varepsilon}_{i,N}(t)-X^{\varepsilon}_{j,N}(t)\right)dt+\sqrt{\varepsilon}dB_{i,N}(t),$
(2.2)

where the Brownian motions $B_{i,N}(t)$ are independent and $F$ represents the
interaction between particles.

The solution $X^{\varepsilon}_{i,N}$ is regarded as the motion of the $N$-th
particle in a random environment that interacts with all other particles
therein. The solution $X^{\epsilon}$ can be seen as the motion of a particle
in a macroscopic system as $N$ in equation (2.2) is the same magnitude as
Avogadro’s constant. The motion of $X^{\epsilon}$ is generated by the vector
field $V$, the interaction $F$, and the random fluctuation $B_{t}$.

We are interested in the transition behaviors of $X^{\epsilon}$ when the
vector field $V$ admits multiple equilibrium states, that is $V=0$ at these
states. Will these states become the metastable states of $X^{\epsilon}$ or
will the interaction $F$ create a different metastable state? What does the
most likely transition path look like? The large deviations theory has shown
its power in tackling such kinds of problems. We will briefly introduce the
large deviations theory in the following subsection.

###  2.1 The large deviations theory

Let $C([0,T])$ denote the space of continuous functions on the interval
$[0,T]$ taking values in $\mathbb{R}^{d}$, equipped with the uniform norm.
Roughly speaking, if the solution $X^{\epsilon}$ to (2.1) satisfies the large
deviations principle with an action functional $I_{T}$ on $C([0,T])$, we can
estimate the probability that $X^{\epsilon}$ stays in a neighborhood of a path
$\varphi$ as follows,

$\mathbb{P}\\{\|X^{\epsilon}-\varphi\|<\delta\\}\sim\exp\left(-\frac{1}{\epsilon}I_{T}(\varphi)\right),\quad\textrm{as
}\epsilon\to 0.$

The above asymptotic results can be used to estimate the probability of events
associated with $X^{\epsilon}$ as $\epsilon$ tends to $0$. To find the most
likely transition path of $X^{\epsilon}$ from point $x_{1}$ to point $x_{2}$
within time $T$, we only need to solve the following constrained minimization
problem. That is

$I_{T}(\varphi^{*})=\min_{\varphi}I_{T}(\varphi),$ (2.3)

where the minimization is constrained by

$\varphi(0)=x_{1},\quad\varphi(T)=x_{2},\quad\varphi\in C([0,T]).$

The function $\varphi^{*}_{T}$ gives the most likely transition path within
time $T$. The infimum value of equation (2.3) over $T\in[0,\infty]$ is called
quasi-potential $\tilde{V}(x_{1},x_{2})$ between $x_{1}$ and $x_{2}$, that is

$\tilde{V}(x_{1},x_{2})=\inf_{T>0}\inf_{\varphi\in\bar{C}_{x_{1}}^{x_{2}}([0,T])}I_{T}(\varphi),$
(2.4)

where $\bar{C}_{x_{1}}^{x_{2}}(0,T)$ is the space of all absolutely continuous
functions that start at $x_{2}$ and end at $x_{2}$. Suppose that there exists
a path $\tilde{\varphi}$ and time $\tilde{T}\in[0,\infty]$ satisfying,

$I_{\tilde{T}}(\tilde{\varphi})=\tilde{V}(x_{1},x_{2}).$

The path $\tilde{\varphi}$ is the most likely transition path from point
$x_{1}$ to $x_{2}$ that we are looking for.

As is shown above, the action functional $I_{T}$ governs the dynamical
behavior of a stochastic system in the context of the large deviations theory.
With results derived in [1], we are able to calculate the action functional
$I_{T}$ for our McKean–Vlasov equation (2.1).

Suppose the coefficient $b(x,t,\mu_{t}):=V(x)-F*\mu_{t}(x)$ satisfies [1,
Assumption 3.2], where $\mu_{t}$ is the law of a random variable. Let $\eta$
be the solution to the following ordinary differential equation,

$\dot{\eta}(t)=V(\eta(t))-F*\delta_{\eta(t)}(\eta(t))=V(\eta(t))-F(0),\quad\eta(0)=x_{1}.$
(2.5)

Equation (2.5) is often regarded as the skeleton equation. According to [1,
Theorem 4.4], we know that $X^{\epsilon}$ satisfies the large deviations
principle in $C([0,T])$ equipped with the topology of the uniform norm with
the action functional

$I^{x_{1}}_{T}(\varphi)=\frac{1}{2}\int_{0}^{T}\left|\dot{\varphi}(s)-V(\varphi(s))+F(\varphi(s)-\eta(s))\right|^{2}ds,$
(2.6)

where $\varphi$ is an absolutely continuous function. $X_{t}^{\epsilon}$
shares the same action functional with the following corresponding stochastic
differential equation if it satisfies the large deviations principle,

$dZ^{\epsilon}_{t}=\left(V(Z^{\epsilon}_{t})-F(Z^{\epsilon}_{t}-\eta(t))\right)dt+\sqrt{\epsilon}dB_{t}.$
(2.7)

Note that the drift term of (2.7) may involve time $t$ explicitly.

###  2.2 Dynamical behaviors of the McKean–Vlasov stochastic differential
equation

We are interested in the transition behavior of the McKean–Vlasov equation
between the equilibrium stable states of the vector field $V$. Based on the
large deviation theory, we will analyze the dynamical behavior of a class of
McKean–Vlasov stochastic differential equations in this subsection.

Let $x_{1}$ and $x_{2}$ be two equilibrium stable states of the vector field
$V$, that is $V(x_{1})=V(x_{2})=\mathbf{0}$. If the interaction $F$ is
$\mathbf{0}$ at point $\mathbf{0}$, equation (2.5) starting at $x_{1}$ reduces
to

$\dot{\eta}(t)=V(\eta(t)),\quad\eta(0)=x_{1}.$ (2.8)

When we look for the most likely transition path from $x_{1}$ to $x_{2}$, the
solution $\eta$ stays at $x_{1}$. Then, the action functional in this case is

$I^{x_{1}}_{T}(\varphi)=\frac{1}{2}\int_{0}^{T}\left|\dot{\varphi}(s)-V(\varphi(s))+F(\varphi(s)-x_{1})\right|^{2}ds=\int_{0}^{T}|\dot{\varphi}-\bar{b}(\varphi(s))|^{2}ds.$

The corresponding stochastic differential equation (2.7) then becomes

$dZ^{\epsilon}_{t}=\left(V(Z^{\epsilon}_{t})-F(Z^{\epsilon}_{t}-x_{1})\right)dt+\sqrt{\epsilon}dB_{t}=\bar{b}(Z_{t}^{\epsilon})dt+\sqrt{\epsilon}dB_{t},\quad
Z_{0}^{\epsilon}=x_{1}.$ (2.9)

For the solution $\xi$ to the following equation,

$\dot{\xi}(t)=V(\xi(t))-F(\xi(t)-x_{1}),\quad 0\leq t<\infty,$ (2.10)

we know that $I_{t}^{x_{1}}(\xi)=0$. Because we can choose $\xi(0)$ as the
starting point and $\xi(t)$ as the end point. This is also the minimum value
of the action functional $I_{t}^{x_{1}}$ on absolutely continuous functions
from $\xi(0)$ to $\xi(t)$. So the quasi-potential between every two different
points on the path $\xi$ is $0$. Then the trajectories of $\xi$ starting at
different points form a set of equipotential lines. In addition, we assume
that the equilibrium stable states $x_{1}$ and $x_{2}$ are also asymptotic
stable states of the system (2.8). Then for the following system

$dY_{t}^{\epsilon}=V(Y_{t}^{\epsilon})dt+\sqrt{\epsilon}dB_{t},$ (2.11)

$x_{1}$ and $x_{2}$ become two metastable states. It means that $Y^{\epsilon}$
starting near one of these two states, will stay a fairly long time near this
state. But for the solution $Z^{\epsilon}$ to the corresponding stochastic
differential equation (2.9), $x_{1}$ and $x_{2}$ may not be metastable states
due to the presence of the interaction term $F$. Instead, we need to find the
equilibrium states of equation (2.10) to determine the metastable states of
$Z^{\epsilon}$. As is shown in the next section, the most likely transition
path will pass by a saddle point and be attracted to a “interaction-generated”
metastable state of $Z^{\epsilon}$.

Note that the integrand $\bar{b}$ does not involve time $t$ explicitly. We can
use a reparametrization method to find the most likely transition path over
time $T\in[0,\infty]$, such as the geometric minimum action method [13]. For
the cases that the drift term $b$ in equation (2.7) involves time $t$
explicitly, the adaptive minimum action method [14] can be used to compute the
most likely transition path.

##  3 Compute the most likely transition path for the McKean-Vlasov stochastic
differential equation

In this section, we will introduce the numerical method based on the adaptive
minimum action method [14] for computing the most likely transition path for
the distribution-dependent system (2.1).

The idea of the adaptive minimum action method(aMAM) is to find the minimizer
of (2.4), by numerically solving an associated Euler-Lagrange equation with a
large enough time $T$ for constrained minimization problem (2.3). The moving
mesh strategy is used to choose a proper mesh grid via a monitor function in
every iteration. The associated Euler-Lagrange equation is then solved on such
mesh grids. Choose a large enough $T$ and the numerical solution is the most
likely transition path we are looking for.

To solve the constrained minimization problem (2.3), we will numerically solve
the associated Euler-Lagrange equation whose solution is regarded as the
steady state of the gradient flow $\delta I_{T}/\delta\varphi$ Specifically,
the Euler-Lagrange equation for action functional (2.6) is

$\left\\{\begin{array}[]{l}\varphi_{tt}-\nabla_{x}b(\varphi,\eta)\varphi_{t}+\big{(}\nabla_{x}b(\varphi,\eta)\big{)}^{\top}(\varphi_{t}-b(\varphi,\eta))-\nabla_{y}b(\varphi,\eta)b(\eta,\eta)=0,\\\
\varphi(0)=x_{1},\quad\varphi(T)=x_{2}\end{array}\right.$ (3.1)

where $\eta$ is the solution of equation (2.5) and $b(x,y)=V(x)-F(x-y)$.

The adaptive mesh grids for solving equation (3.1) are determined by a monitor
function. In this paper, we adopt the monitor function in [14] and it is given
by $w(s)=|b(\varphi(s),\eta(s))|^{r}/C$ with
$C=\int_{0}^{T}|b(\varphi(s),\eta(s))|^{r}ds$. The new variable $\alpha$ is
then given by $\alpha(s)=\int_{0}^{s}w(\tau)d\tau$. Denote
$\tilde{\varphi}(\alpha(s))=\varphi(s)$. Then, with respect to parameter
$\alpha$, the Euler-Lagrange equation (3.1) becomes

$\left\\{\begin{array}[]{l}w^{2}\tilde{\varphi}^{\prime\prime}+\biggl{(}\big{(}\nabla_{x}b(\tilde{\varphi},\eta)\big{)}^{\top}-\nabla_{x}b(\tilde{\varphi},\eta)+w_{\alpha}\biggr{)}w\tilde{\varphi}^{\prime}+\nabla_{x}b(\tilde{\varphi},\eta)b(\tilde{\varphi},\eta)-\big{(}\nabla_{y}b(\tilde{\varphi},\eta)\big{)}^{\top}b(\eta,\eta)=0,\\\
\tilde{\varphi}(0)=x_{1},\quad\tilde{\varphi}(T)=x_{2},\end{array}\right.$

where $\tilde{\varphi}^{\prime}$ and $w_{\alpha}$ is the derivative with
respect to parameter $\alpha$. The action functional (2.6) then becomes

$I_{T}(\varphi)=I(\tilde{\varphi})=\frac{1}{2}\int_{0}^{1}|\dot{\tilde{\varphi}}(\alpha)w(s)-V(\tilde{\varphi}(\alpha))+F(\tilde{\varphi}(\alpha)-\tilde{\eta}(\alpha))|^{2}\frac{d\alpha}{w(s)}.$

The algorithm is given in Algorithm 1.

Algorithm 1 An algorithm for distribution-dependent equations

1:At the k-th iteration, given $(t_{i}^{k},\varphi_{i}^{k},\eta_{i}^{k})$ for
$i=0,1,\cdots,N$, calculate
$w_{i}^{k}=|b(\varphi_{i}^{k},\eta_{i}^{k})|^{r}/C_{k}$ with
$C_{k}=\sum_{i}(t_{i+1}^{k}-t_{i}^{k})|b\big{(}(\varphi_{i+1}^{k}+\varphi_{i}^{k})/2,(\eta_{i+1}^{k}+\eta_{i}^{k})/2\big{)}|^{r}$.

2:Calculate
$\Delta\alpha_{i}^{k}=\frac{(t_{i+1}^{k}-t_{i}^{k})\left|b\left((\varphi_{i+1}^{k}+\varphi_{i}^{k})/2,(\eta_{i+1}^{k}+\eta_{i}^{k})/2\right)\right|^{r}}{C_{k}}.$

3:Set $\Delta\alpha=1/N$ and $\alpha_{i+1}-\alpha_{i}=\Delta\alpha$ for
$i=0,\cdots,N$ and $\alpha_{0}=0$.

4:Interpolate $(\alpha_{i}^{k},\varphi_{i}^{k})$ to get
$(\alpha_{i},\bar{\varphi}_{i}^{k})$, and interpolate
$(\alpha_{i}^{k},t_{i}^{k})$ to get $(\alpha_{i},t_{i}^{k+1},\eta_{i}^{k+1})$,
where $\eta_{i}^{k}$ is generated by the Euler scheme at time $t_{i}^{k}$.

5:Let $\\{\bar{\varphi}_{i}\\}_{i=1}^{N+1}$ be the solution of the following
linear equations  $\left\\{\begin{aligned}
\frac{\bar{\varphi}_{i}-\bar{\varphi}_{i}^{k}}{\Delta\tau}=&(w_{i}^{k})^{2}\frac{\bar{\varphi}_{i+1}-2\bar{\varphi}_{i}+\bar{\varphi}_{i-1}}{\Delta\alpha^{2}}+\bigg{(}\nabla_{y}b(\bar{\varphi}_{i}^{k},\eta_{i}^{k+1})^{\top}-\nabla_{y}b(\bar{\varphi}_{i}^{k},\eta_{i}^{k+1})+(w_{\alpha})_{i}^{k}\textrm{Id}\bigg{)}w_{i}^{k}\bar{\varphi}_{i}^{\prime
k}\\\
&-\nabla_{y}b(\bar{\varphi}_{i}^{k},\eta_{i}^{k+1})b(\bar{\varphi}_{i}^{k},\eta_{i}^{k+1})-\nabla_{x}b(\bar{\varphi}_{i}^{k},\eta_{i}^{k+1})^{\top}b(\eta_{i}^{k+1},\eta_{i}^{k+1})\\\
\bar{\varphi}_{0}=x_{1},\quad&\bar{\varphi}_{N}=x_{2}.\end{aligned}\right.$
where $\bar{\varphi}_{i}^{\prime
k}=(\bar{\varphi}_{i+1}^{k}-\bar{\varphi}_{i-1}^{k})/(2\Delta\alpha)$ and
$(w_{\alpha})_{i}^{k}=(w_{i+1}^{k}-w_{i-1}^{k})/(2\Delta\alpha)$

6:Repeat step 1-5 until a stopping criterion is fulfilled.

We will compute the most likely transition path for the Maier-Stein model with
different interaction terms. That is

$dX_{t}^{\epsilon}=V\left(X_{t}^{\epsilon}\right)dt-F*u_{t}^{\epsilon}(X_{t}^{\epsilon})dt+\sqrt{\epsilon}dB_{t},\quad
X_{0}^{\epsilon}=(-1,0),$

where

$V(u,v)=\left(\begin{array}[]{c}u-u^{3}-\beta uv^{2}\\\
-\left(1+u^{2}\right)v\end{array}\right),\quad\beta=10.$

Figure 1: The most likely transition path for the Maier-Stein model with
interaction $F_{L}$. The orange line is the most likely transition path
connecting points $(-1,0)$ and $(1,0)$. Grey lines are the equipotential lines
of the associated quasi-potential. We set the number of points $N=200$,
iteration $K=5000$ and time $T=20$. The initial path is chosen as
$y=-0.5x^{2}+0.5$.

In Figure 1, we calculate the most likely transition path connecting two
equilibrium stable states $(-1,0)$ and $(1,0)$, under the interaction
$F=F_{L}$, which is

$F_{L}(u,v)=\left(\begin{array}[]{c}u\\\ -v\end{array}\right).$

The corresponding stochastic differential equation (2.9) then becomes

$dZ^{\epsilon}_{t}=b_{L}(Z_{t}^{\epsilon})dt+\sqrt{\epsilon}dB_{t},\quad
Z_{0}^{\epsilon}=(-1,0),$

where

$b_{L}(u,v)=\left(\begin{array}[]{c}-u^{3}-\beta uv^{2}+1\\\
-u^{2}v\end{array}\right).$

The action functional is

$I_{L}(\varphi)=\int_{0}^{T}|\dot{\varphi}-b_{L}(\varphi(s))|^{2}ds.$

There is only one asymptotic stable point $(-1,0)$ for the following equation,

$\dot{\xi_{L}}(t)=b_{L}(\xi_{L}(t)).$

It is also the metastable point of the corresponding stochastic differential
equation. The grey lines in Figure 1 are the orbits of $\xi_{L}$ starting at
different points and these orbits form the equipotential lines of the quasi-
potential. The interaction $F$ eliminates the metastable point $(1,0)$. The
most likely transition path in this case is actually the most likely exit path
from $(-1,0)$ to $(1,0)$.

Figure 2: The most likely transition path for the Maier-Stein model with
interaction $F_{BS}$. The orange line is the most likely transition path
connecting points $(-1,0)$ and $(1,0)$. Grey lines are the equipotential lines
of the associated quasi-potential. The dark grey point is the saddle point of
equation (3.2). We set the number of points $N=400$, iteration $K=15000$ and
time $T=40$. The initial path is chosen as $y=-0.5x^{2}+0.5$.

In Figure 2, we calculate the most likely transition path under the
interaction described by a modified Biot-Savart kernel which is not singular
at the origin. Precisely, we take $F=F_{BS}$ to be

$F_{BS}(u,v)=\frac{1}{2\pi}\left(\begin{array}[]{c}\frac{-v}{u^{2}+v^{2}+\delta}\\\
\frac{u}{u^{2}+v^{2}+\delta}\end{array}\right),\quad\delta=0.01.$

The corresponding stochastic differential equation (2.9) then becomes

$dZ^{\epsilon}_{t}=b_{BS}(Z_{t}^{\epsilon})dt+\sqrt{\epsilon}dB_{t},\quad
Z_{0}^{\epsilon}=(-1,0),$

where

$b_{BS}(u,v)=\left(\begin{array}[]{c}u-u^{3}-\beta
uv^{2}+\frac{v}{2\pi((u+1)^{2}+v^{2}+\delta)}\\\
-(1+u^{2})v-\frac{u+1}{2\pi((u+1)^{2}+v^{2}+\delta)}\end{array}\right).$

The action functional is

$I_{BS}(\varphi)=\int_{0}^{T}|\dot{\varphi}-b_{BS}(\varphi(s))|^{2}ds.$

For the presence of the interaction term $F_{BS}$, there are two asymptotic
stable points $x_{1}=(-1,0)$, $x_{2}\approx(1.14801,0.0318893)$ and a saddle
point $x^{*}\approx(0.009145,-0.152703)$ for the following equation,

$\dot{\xi}_{BS}(t)=b_{BS}(\xi_{BS}(t)).$ (3.2)

The metastable points for the corresponding stochastic differential equation
are then $x_{1}$ and $x_{2}$. The grey lines in Figure 2 are the orbits of
$\xi_{BS}$ starting at different points and these form the equipotential lines
of the quasi-potential. The saddle point $x^{*}$ is the dark grey point in
Figure 2. The interaction $F$ transports the saddle point $(0,0)$ and stable
point $(1,0)$ to $x^{*}$ and $x_{2}$ respectively. Then, the most likely
transition path will pass by the saddle point $x^{*}$ and be attracted to the
neighborhood of $x_{2}$ as is shown in Figure 2.

##  4 Conclusion

We have investigated the most likely transition path for a class of
distribution-dependent stochastic dynamical systems. In the small noise
regime, this problem turns out to be equivalent to computing this transition
path for a stochastic dynamical system without distribution-dependency.
Therefore, we can apply the adaptive minimum action method to calculate the
most likely transition paths, as demonstrated in several examples.

## Acknowledgements

This work was partly supported by the NSFC grants 11771449.

## Data Availability

The data that support the findings of this study are openly available in
GitHub: https://github.com/JayWeiess/MVSDE-aMAM

## References

  * [1] G. Reis, W. Salkeld and J. Tugaut.  Freidlin–Wentzell LDP in path space for McKean–Vlasov equations and the functional iterated logarithm law.  The Annals of Applied Probability, 29:1487 – 1540, 2019
  * [2] A. Dembod and O. Zeltouni.  Large Deviations Techniques and Applications.  Springer, Berlin, Heidelberg, 2010.
  * [3] M. I. Freidlin and A. D. Wentzell.  Random Perturbations of Dynamical Systems.  Springer, Berlin, Heidelberg, 2012
  * [4] S. Herrmann, P. Imkeller and D. Peithmann.  Large deviations and a Kramers’ type law for self-stabilizing diffusions.  The Annals of Applied Probability, 18:1379 – 1423, 2008.
  * [5] F.Bolley, A. Guillin and C. Villani.  Quantitative concentration inequalities for empirical measures on non-compact spaces.  Probability Theory and Related Fields, 137:541 – 593, 2007.
  * [6] F. Cucker and S. Smale.  Emergent behavior in flocks.  IEEE Transactions on Automatic Control, 52:852– 862, 2007.
  * [7] P.E. Jabin.  A review for the mean field limit for Vlasov equations.  Kinetic and Related Models, 7:661– 711, 2014.
  * [8] H. P. McKean.  A class of Markov processes associated with nonlinear parabolic equations.  Proceedings of the National Academy of Sciences, 56: 1907 – 1911, 1966.
  * [9] R. Carmonaand F. Delarue.  Probabilistic Theory of Mean Field Games with Applications I.  Springer, Cham, 2017.
  * [10] R. Carmonaand F. Delarue.  Probabilistic Theory of Mean Field Games with Applications II.  Springer, Cham, 2017.
  * [11] S. B. Pope.  Turbulent Flows.  Cambridge University Press, 2000.
  * [12] W. Liu, Y. Song, J. Zhai and T. Zhang.  Large and moderate deviation principles for McKean-Vlasov SDEs with jumps.  arXiv preprint arXiv:2011.08403 2020
  * [13] M. Heymann and E. Vanden-Eijnden.  The geometric minimum action method: A least action principle on the space of curves.  Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 61: 1052–1117, 2008.
  * [14] Y. Sun and X. Zhou.  An improved adaptive minimum action method for the calculation of transition path in non-gradient systems.  Communications in Computational Physics, 24: 44–68, 2018.
  * [15] X. Zhou, W. Ren and W. E.  Adaptive minimum action method for the study of rare events.  The Journal of Chemical Physics, 128: 104–111, 2008.
  * [16] W. E, W. Ren and E. Vanden-Eijnden.  Minimum action method for the study of rare events.  Communications on Pure and Applied Mathematics, 57: 637–656, 2004.
  * [17] X. Wan, H. Yu and J. Zhai.  Convergence analysis of a finite element approximation of minimum action methods.  SIAM Journal on Numerical Analysis, 56: 1597–1620, 2018.


# Deciphering Star Cluster Evolution by Shape Morphology

Qingshun Hu  Xinjiang Astronomical Observatory, Chinese Academy of Sciences,
Urumqi, Xinjiang 830011, People’s Republic of China  University of Chinese
Academy of Sciences, Beijing 100049, People’s Republic of China  Yu Zhang
<EMAIL_ADDRESS>Xinjiang Astronomical Observatory, Chinese Academy of Sciences,
Urumqi, Xinjiang 830011, People’s Republic of China  University of Chinese
Academy of Sciences, Beijing 100049, People’s Republic of China  Ali Esamdin
<EMAIL_ADDRESS>Xinjiang Astronomical Observatory, Chinese Academy of
Sciences, Urumqi, Xinjiang 830011, People’s Republic of China  University of
Chinese Academy of Sciences, Beijing 100049, People’s Republic of China
Jinzhong Liu  Xinjiang Astronomical Observatory, Chinese Academy of Sciences,
Urumqi, Xinjiang 830011, People’s Republic of China  University of Chinese
Academy of Sciences, Beijing 100049, People’s Republic of China  Xiangyun Zeng
Xinjiang Astronomical Observatory, Chinese Academy of Sciences, Urumqi,
Xinjiang 830011, People’s Republic of China  University of Chinese Academy of
Sciences, Beijing 100049, People’s Republic of China

###### Abstract

We analyze the morphological evolution of open clusters and provide shape
parameters for 265 open clusters. The results show that the overall shape of
sample clusters becomes more elliptical as they grow older, while their core
remains circular or slightly trend to circularize. There is a negative
correlation of the ellipticities with the number of members of the sample
clusters. A significant negative correlation between the overall ellipticities
and masses is also detected for the sample clusters with log(age/year) $\geq$
8, suggesting that the overall shapes of the clusters are possibly influenced
by the number of members and masses, in addition to the external forces and
the surrounding environment. For most young sample clusters, the radial
stratification degree of the short axis direction is greater than that of the
long, implying that the radial stratification degree in the two directions
within the young sample cluster may be unevenly affected by an internal
evolutionary process. Older sample clusters exhibit lower stratification in
the tangential direction, which possibly means those clusters may continue to
survive for a long time at a low level of stratification. Our analysis shows
that the overall shape of the sample clusters may be more susceptible to the
influence of Galactic tides toward the Galactic center than the shear forces
embedded in Galactic differential rotation. By analyzing the distribution of
the ages and number of members of star clusters, we suggest that NGC 6791 may
originate from superclusters.

Galaxy: stellar content — methods: statistical — open clusters and
associations: general

aasjournal

##  1 Introduction

Stars form in groups out of molecular clouds (e.g. Lada & Lada, 2003; Alves &
Bouy, 2012; Megeath et al., 2016). The youngest groups are distributed
collectively as imprints of the usual filamentary cloud structures. Once
surplus gas and dust is being dispersed (e.g. Tutukov, 1978; Hills, 1980;
Elmegreen, 1983; Mathieu, 1983; Adams, 2000; Kroupa et al., 2001; Kroupa,
2002; Goodwin & Bastian, 2006; Moeckel & Bate, 2010; Pelupessy & Portegies
Zwart, 2012; Brinkmann et al., 2017), two-body gravitational relaxation (e.g.
Bergond et al., 2001; Ryon et al., 2017; Mamon et al., 2019) among stellar
members leads to sphericalization, i.e., an overall circular shape, with
lower-mass members occupy larger volumes (“mass segregation”). The least-
massive members are most vulnerable to dynamically ejection from the system
(“stellar evaporation”), further loosening the gravitational binding (e.g.
Gieles & Baumgardt, 2008; Tang et al., 2019). External perturbation such as
the tidal forces, differential rotation, or disk crossing, exacerbates the
disintegration process, and eventually, the members would supply as Galactic
field stars. The remnant systems, globular or open clusters we see today, with
relatively well-determined distances, ages, and metallicities, serve to
diagnose how stars of different masses evolve, and as laboratories to test
stellar dynamics theories.

Open clusters, in particular, with dozens to thousands of member stars each
(Bergond et al., 2001), with ages from millions to billions of years, and
located in the Galactic disk, provide a sample to probe the Galactic structure
and evolution (Balaguer-Núñez et al., 2020). Most star groups are likely
dispersed upon emergence from molecular clouds (Lada & Lada, 2003), and the
majority of the survival star clusters are expected to be dissolved in some
$10^{8}$ years (e.g. Wielen, 1971; Pandey & Mahra, 1986; Bergond et al., 2001;
Bonatto et al., 2006; Yang et al., 2013). Conceivably, longevity is possible
only among systems relatively free of disturbances.

The dynamical evolution of an open cluster, arising from internal or external
influences, manifests itself in its shape changes. As early as more than a
century ago, the flattening of a moving cluster was reasoned (Jeans, 1916),
and such a phenomenon was indeed detected in the Hyades (Oort, 1979). Pandey
et al. (1990) identified the existence of ”corona” features around open
clusters. Such halo structures were characterized by Nilakshi et al. (2002)
based on the photometric data of the Digital Sky Survey. Bergond et al. (2001)
found open cluster to be generally elongated in a direction parallel to the
plane, attributed to tidal forces. Chen et al. (2004), with a sample of 31
open clusters from the 2MASS (Two Micron All Sky Survey) dataset, found
clusters to be elongated even among the youngest systems, reflecting
morphological relics in the parental molecular clouds; namely, young open
clusters inherit the shape of the dense molecular cloud from which they formed
(Getman et al., 2018). With the all-sky data of the compiled catalog of 2.5
million stars ASCC-2.5, Kharchenko et al. (2009) derived the shape parameters,
that is, the apparent ellipticity and orientation of the ellipse, of 650 open
clusters and concluded that the ellipticities of open clusters generally range
from 0.2 to 0.3, which are lower than those of modeled clusters. Zhai et al.
(2017) computed the morphology parameters of 154 open clusters using the WEBDA
database and suggested all clusters in their sample are elongated in shape.
Bhattacharya et al. (2017), by using probabilistic star counting of Pan-STARRS
point sources, detected the core-tail structure in the aged open cluster
Berkeley 17, with the tails preferentially consisting of low-mass members.
With Gaia Data Release 2 (DR2) data, the substructure and halo of the double
cluster $h$ and $\chi$ Persei was diagnosed by (Zhong et al., 2019),
Similarly, tidal tails are found in the nearby clusters Coma Berenices (Tang
et al., 2019) and Blanco 1 (Zhang et al., 2020). The expected avalanche of
high-precision kinematic and photometric data, particularly from Gaia, affords
reliable membership determination in, and hence the studies of the dynamical
evolution of, star clusters.

In this paper, we analyze a sample of 265 open clusters with which membership
is available in the literature to diagnose the possible dominant processes in
the history of the dynamical evolution of the open clusters in general. We
present in Sec. 2 the data used in our study and the methodology to derive the
shape parameters of the clusters. The results and the discussions are
presented in Sec. 3. A summary of this work is given in Sec. 4.

##  2 Data and Method

The data used in this work is from the Gaia DR2 (see, e.g. Gaia Collaboration
et al., 2016, 2018), a dataset of some 1.3 billion stars (e.g. Gaia
Collaboration et al., 2018; Arenou et al., 2018; Cantat-Gaudin et al., 2018)
with precise astrometry at the sub-milliarcsecond level (Lindegren et al.,
2018) and homogeneous photometry at the mmag level (Evans et al., 2018). The
uncertainties of proper motion are about 0.05 mas yr-1 for bright sources (G
$<14$ mag) and 1.2 mas yr-1 for faint ones (G$\sim 20$ mag). The uncertainties
of parallax are about 0.03 mas and 0.7 mas for bright and faint sources,
respectively (Lindegren et al., 2018). Reliable membership identification in
the open clusters, essential in our work, is made possible by the high-quality
Gaia data.

We started with the member catalog of open clusters (Cantat-Gaudin et al.,
2018), which provides the combined spatio-kinematic-photometric membership of
1229 open clusters based on the UPMASK method (Krone-Martins & Moitinho, 2014)
applied on the Gaia DR2. And, we adopted the age parameters of 269 open
clusters from Bossini et al. (2019). After matching these two samples, we
obtained 266 clusters in common. Except for Blanco 1, which is not included in
our analysis because of its high Galactic latitude and too close a distance,
the rest 265 clusters form the sample for our shape analysis, with logarithmic
ages (years) ranging from 7.006 to 9.927, i.e., from relatively young to aged
systems, fulfilling our requirement to diagnose the evolutionary of the open
clusters.

We analyzed the shape of each cluster by using Nonparametric Bivariate Density
Estimation, and Least Square Ellipse Fitting, both available in Python
packages (Astropy, Matplotlib, Numpy, Scipy, and Sympy). Estimating the
stellar number density profile according to the actual density distribution
without model assumptions, the Nonparametric Bivariate Density Estimation is
simple and nonparametric based on kernel density estimation. Its principle, in
detail, is to estimate the probability density function of a random variable
in a nonparametric way, which works for univariate and multivariate data. It
has the advantage of delineating the actual density of the member stars in a
cluster as far as possible. In this work, the kernel of the density estimation
we adopted is Gaussian kernel. The shape parameters are then derived by the
Least Square Ellipse Fitting, including the ellipticity $e$, and the position
angle $q$ between the long principal axis of the ellipse and the Galactic
plane. The value of ellipticity $e$, ranging between 0 and 1, represents the
shape elongation, whereas $q$ refers to the orientation of a cluster in the
disk, ranging from $-90\arcdeg$ to $90\arcdeg$.

Figure 1: The general scheme of fitting cluster density profiles for four open
clusters (NGC 7235, NGC 752, ASCC 16, and Haffner 22). For each cluster in
each panel, the black dots represent the member stars of a cluster with the
gray lines as the stellar number density profiles at 50% of and 5% of the
peak. the fitted ellipses of 50% and 5% of the peak stellar number density of
the clusters are marked in blue and red, respectively. The percentage labels
on the gray lines represent the relative percentage of the actual density
value relative to the peak. The density contours at 5% in the panel are
incomplete because the cluster member stars show a non-dispersive distribution
at the periphery of them.

The shape parameters are derived in two steps. First, the members of each
cluster with center $(l_{0},b_{0})$ given in Cantat-Gaudin et al. (2018) are
plotted in a two-dimensional spherical Galactic coordinate system $(l,b)$. The
density profile is then overplotted with the Nonparametric Bivariate Density
Estimation. Figure 1 illustrates the results for four clusters, NGC 7235, NGC
752, ASCC 16, and Haffner 22. In estimating the densities, this work uses the
maximum and minimum values of cluster member coordinate data as grid
boundaries. The critical factor in the density profile is the grid size that
dictates the degree of smoothness of the density profile. Empirically we found
the number of grid cells of 40000 as a satisfactory compromise between
resolution and possible density fluctuations. Second, with the density
profile, we assign two ellipses, via the Least Square Ellipse Fitting, to
represent the shapes of the core and the overall of a cluster. For this, we
set arbitrarily at 50% and 5% of the peak stellar number density of each open
cluster. These ellipses are parameterized with $e_{core}$, $q_{core}$ for the
inner shape, and $e_{all}$, $q_{all}$ for the overall shape of a cluster.
Meanwhile, the additional parameters such as $a_{core}$ (half-length axes of
inner ellipises), $a_{all}$ (half-length axes of overall ellipises),
$b_{core}$ (half-short axes of inner ellipises), and $b_{all}$ (half-short
axes of overall ellipises) are provided. We estimate the errors of
ellipticities by using the standard deviation of the ellipse fitted. The
median of relative errors of the inner ellipticities ($e_{core}$,
corresponding to the inner shape of the clusters) for the sample clusters is
about 4%, and that of the overall ($e_{all}$, corresponding to the overall
shape of the clusters) about 11%. We note that one cluster, Platais 10, did
not display a reliable solution with our analysis pipeline for the overall
ellipticity due to the scattered distribution of its members; the overall
shape, therefore, was estimated by visual inspection for it. The shape
parameters and additional parameters of our sample are listed in Table 1.

##  3 Results and discussion

As shown in Figure 2, the ellipticities and orientations of the cores and
overalls of the whole sample clusters are presented in the Galactic Coordinate
System. The clusters in our sample are distributed around the Galactic plane
within the range of the Galactic latitude ($b$) from -40o to 40o. The
histograms of the orientations of the cluster’s cores and overalls are
presented in Figure 3, which shows that the orientations of most clusters in
the sample are within $|q|$ $<$ 45o, which is basically consistent with the
finding of Zhai et al. (2017). This means that the sample of open clusters is
more elongated in a direction parallel to the Galactic plane than
perpendicular to it, consistent with the finding of Bergond et al. (2001).

Figure 2: The distributions of the sample clusters in Galactic coordinate
system. The lengths of the blue and red bars are proportional to the
ellipticities of the cores and overalls of the clusters, respectively. The
inclinations of the blue and red bars indicate orientations of the cores and
overalls of the clusters, respectively. The midpoints of the bars are set to
the cluster’s positions. Figure 3: Distribution of the orientations
($q_{core}$ and $q_{all}$) of this work. The blue histogram (left panel)
corresponds to $q_{core}$ and the red histogram (right panel) to $q_{all}$ of
the sample clusters.

###  3.1 Correlations between Ellipticities and the Ages, Number of member
stars, and Masses of the Clusters

Figure 4: The ellipticities ($e_{core}$ and $e_{all}$) as a function of the
ages of the sample clusters. The inner ellipticities ($e_{core}$) are shown in
the top panel, and the overall ($e_{all}$) in the bottom. The gray bars show
the errors of the ellipticities. In each panel, the corresponding linear fit
is shown as a red line.

With the ellipticities ($e_{core}$ and $e_{all}$) of the sample clusters
obtained, we investigate the correlations between the ellipticities and the
ages, the number of the member stars, and the masses of the clusters. Figure 4
shows the ellipticities ($e_{core}$ and $e_{all}$) as a function of the ages
of the sample clusters. The figure indicates that with age, the overall of the
open clusters become more elliptical while their inner cores remains circular
or slightly trend to circularize, which is consistent with what was inferred
by Chen et al. (2004); Zhai et al. (2017). It should be expected that the core
of a star cluster gets gradually circularised with age due to stellar
dynamics, whereas the overall shape of a cluster continues to be elongated by
external forces such as the Galactic tides, the Galactic differential
rotation, or encountering with giant molecular clouds and so on. We should
note that the significance of the linear fit to our sample is not strong,
moreover, the shape parameters fitted possibly deviate from their true value
due to the projection effect. More samples with three-dimensional information
are needed to determine the true morphological evolution of the clusters over
time.

Figure 5: The ellipticities ($e_{core}$ and $e_{all}$) as a function of the
number of the member stars for the sample clusters. The symbols and bars in
the two panels are the same as Figure 4. In each panel, the corresponding
linear fit is shown by a red line.

Figure 5 shows the ellipticities of the sample clusters ($e_{core}$ and
$e_{all}$) as a function of the number of member stars. From this figure, it
can be concluded that both the inner ellipticities and the overall
ellipticities of the clusters are negatively correlated with the number of the
member stars of the clusters, which may imply that the higher the number of
member stars of the clusters, the stronger the cluster’s gravitational binding
to resist the external disturbances, and the more it can keep smaller
ellipticity. Therefore, clusters with more members should be able to go
through a relatively stable shape evolutionary and thus survive for a long
time in an arcane and complex dynamical environment. Conversely, the clusters
with fewer members are more likely to disintegrate into field stars.

Figure 6: The masses versus the ellipticities ($e_{core}$ and $e_{all}$) of
100 clusters in our sample. the subsample is divided into two age groups
(log(age/year) $<$ 8 shown as red open circles and log(age/year) $\geq$ 8 as
blue). In each panel, the corresponding linear fits are shown by the red and
blue lines, respectively. Figure 7: The ellipticities ($e_{core}$ and
$e_{all}$) vs. the heights from the Galactic disk for two age groups from the
whole sample (the left and middle panels), and, the histogram of the heights
of them (the right panel). The vertical axis of the three panels is the
absolute value of the heights. The blue and red circles in the left and the
middle panels correspond to clusters with log(age/year) $\geq$ 8 and
log(age/year) $<$ 8, respectively. The red hatched histogram in the right
panel corresponds to clusters with log(age/year) $<$ 8 and the blue shadows
histogram for clusters with log(age/year) $\geq$ 8.

It is well known that cluster masses are an important factor affecting the
shape of open clusters. Piskunov et al. (e.g. 2007, 2008) has given the masses
of a number of clusters, which are determined based on the tidal radius. To
investigate the correlation between ellipticities and cluster masses, we
obtained the masses of 100 common clusters by cross-matching our sample with
the one given by Piskunov et al. (e.g. 2007, 2008), and then plot the
distribution of the ellipticities versus masses for these 100 clusters, shown
in Figure 6. Previous studies have shown that the majority of the survival
star clusters are expected to be dissolved in some $10^{8}$ years (e.g.
Wielen, 1971; Pandey & Mahra, 1986; Bergond et al., 2001; Bonatto et al.,
2006; Yang et al., 2013). Therefore, we set to log(age/year) = 8 as a grouping
age, then divided the 100 clusters into two age groups, and shows the
distribution between the ellipticity and mass of the two groups in Figure 6.
As shown by the figure, no significant correlations between the inner
ellipticities and masses and between the overall ellipticities and masses are
detected for the clusters with log(age/year) $<$ 8\. However, a negative
correlation between the inner and overall ellipticities with masses is
detected for the clusters with log(age/year) $\geq$ 8, especially more
significant between the overall ellipticities and the masses. This may imply
that older clusters with larger masses may have a large enough self-
gravitational potential to keep their overall ellipticity small and thus
maintain a stable shape evolution, especially when no giant molecular clouds
or other events are encountered.

To investigate why there is no significant correlation between cluster
ellipticity and its mass for log(age/year)$<$ 8 in Figure 7, we divided the
whole sample of 265 clusters into two age groups of log(age/year) $\geq$ 8 and
log(age/year) $<$ 8, and then plotted the ellipticities versus the absolute
values of the heights from the Galactic disk. The left and middle panels of
Figure 7 show the ellipticities ($e_{core}$ and $e_{all}$) versus the heights
for the two age groups, respectively, and the right panel of the figure show
the histogram of the heights of them. As seen in the right panel of Figure 7,
the heights of the older group are clearly more widely distributed than that
of the younger group, with a considerable number of the clusters in the older
group distributed at relatively high heights. This suggests that the younger
group clusters are closer to the disk and maybe more disturbed by the gravity
of the Galactic disk and the surrounding environment, therefore, the mass of
the younger group clusters may not be the dominant factor affecting the
cluster ellipticity, so the correlation between the cluster ellipticity and
mass is not obvious, as shown in Figure 6. However, for the older clusters,
they are more widely distributed in the height from the disk, and a
considerable number of clusters are located in the higher height from the disk
compared to the younger group, thus these clusters may be relatively less
disturbed by the gravity of Galactic disk and the surrounding environment
compared with the young group, so they show a weak negative correlation
between the cluster ellipticity and mass, as shown in Figure 6. There is no
significant difference in the inner ellipticities distribution between the two
age groups (see the left panel of Figure 7). However, the overall
ellipticities distribution of the samples with log(age/year) $\geq$ 8 is
significantly wider than that of the younger group (see the middle panel of
Figure 7), which indicates that the overall ellipticities maybe tend to widen
with time evolution. The above distributions are consistent with the results
of the distribution of the inner ellipticities and the overall ellipticities
with cluster’s ages shown in the top and bottom panels in Figure 4,
respectively.

Table 1: The Parameters of 265 Open Clusters in This Paper (1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | (11) | (12) | (13) | (14) | (15) | (16)
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Name | $l$ | $b$ | log(age/year) | H | Number | $e_{core}$ | $\sigma_{core}$ | $q_{core}$ | $e_{all}$ | $\sigma_{all}$ | $q_{all}$ | $a_{core}$ | $b_{core}$ | $a_{all}$ | $b_{all}$
– | (degree) | (degree) | – | (pc) | – | – | – | (degree) | – | – | (degree) | – | – | – | –
Alessi1 | 123.255 | -13.335 | 8.935 | -162.5 | 63 | 0.423 | 0.034 | 17.388 | 0.219 | 0.196 | -10.205 | 0.310 | 0.179 | 0.691 | 0.540
Alessi2 | 152.333 | 6.365 | 8.759 | 68.5 | 276 | 0.280 | 0.026 | -12.179 | 0.423 | 0.141 | 6.858 | 0.529 | 0.381 | 1.962 | 1.133
Alessi5 | 288.058 | -1.966 | 7.723 | -13.6 | 308 | 0.174 | 0.007 | -29.674 | 0.293 | 0.087 | -1.436 | 0.435 | 0.359 | 1.327 | 0.938
Alessi6 | 313.643 | -5.544 | 8.777 | -85.5 | 294 | 0.097 | 0.010 | -6.524 | 0.256 | 0.196 | 1.431 | 0.280 | 0.252 | 1.078 | 0.802
Alessi8 | 326.501 | 4.259 | 8.137 | 49.2 | 82 | 0.240 | 0.022 | -0.251 | 0.054 | 0.167 | 8.270 | 0.249 | 0.189 | 0.682 | 0.646
Alessi10 | 31.597 | -21.058 | 7.746 | -158.8 | 75 | 0.262 | 0.013 | -6.236 | 0.259 | 0.149 | -10.651 | 0.250 | 0.185 | 0.848 | 0.629
… | … | … | … | … | … | … | … | … | … | … | … | … | … | … | …

###  3.2 The Stratification degree of Hierarchical Structure

Figure 8: The distribution of $\Delta a/a_{core}$ versus $\Delta b/b_{core}$
of the sample clusters. The color shade of the solid circles indicates the
logarithmic ages (years) of the sample clusters. The size of the symbol is
proportional to the number of sample cluster members.

The hierarchical structure of open clusters has been verified (e.g. Raboud &
Mermilliod, 1998a, b; Chen et al., 2004). But its detailed analysis, such as
radial and tangential stratification degree of cluster interior, has been
lacking. In our work, we use shape parameters to study the stratification
degree. We define ($\Delta a/a_{core}$) and ($\Delta b/b_{core}$) as the
radial stratification degree inside the cluster. In detail, ($\Delta
a/a_{core}$) denotes the radial stratification degree in the long-axis
direction, and ($\Delta b/b_{core}$) denotes the radial stratification degree
in the short-axis direction. And $\Delta a$ and $\Delta b$ equal to
$a_{all}-a_{core}$ and $b_{all}-b_{core}$, respectively. Moreover, we define
$q_{all}-q_{core}$ as the tangential stratification within the cluster.

Figure 8 displays the distribution of ($\Delta a/a_{core}$) versus ($\Delta
b/b_{core}$), which illustrates the radial stratification degree within the
cluster. It can be seen from this figure that for those clusters with a high
number of member stars, the radial stratification degree in the short axis
direction is almost the same as that in the long. The radial stratification
range of most clusters with a small number of member stars is larger than that
of clusters with a large number of member stars, both in the short-axis
direction and in the long-axis direction within the clusters. In addition, we
find that the age of the sample clusters is possibly related to the degree of
radial stratification. To explore this issue, we define $p$ as the ratio of
the radial stratification degree in the short axis direction to that in the
long axis direction. In detail, $p$ $>$ 1 shows the radial stratification
degree in the short axis direction is greater than that of the long axis,
while the opposite is true for $p$ $\leq$ 1\. These parameters have complied
in Table 2. For the sample clusters, we find that in about 59% of the
clusters, the radial stratification degree of the short axis is greater than
that of the long. Furthermore, if we count $p$ by different age groups,
finding that in about 67% of young clusters in (log(age/year) $\leq$ 8), the
radial stratification degree in the short axis direction is greater than that
in the long-axis direction, while for the remaining age groups, about half of
the clusters in each age group have greater radial stratification degree in
the short-axis direction than in the long-axis. This manifests that the number
of clusters in the sample with a greater radial stratification in the short-
axis direction than in the long-axis direction is almost the same as the
opposite case. However, for most young sample clusters, the radial
stratification degree of the short axis direction is greater than that of the
long axis.

Strictly speaking, however, we should consider the shape parameters we derive
as lower limits on the true shape parameters of the cluster. As long as a
cluster has the form of a triaxial ellipsoid, we observe the projection of the
cluster rather than its true size. The degree of possible discrepancy depends
on the axial relationship of the ellipsoid and the direction of the axis
relative to the cluster’s line of sight (Piskunov et al., 2008). It is worth
noting that the true shape of clusters in space is currently unknown, and
there may be projection effects in our results. Nonetheless, some plausible
conclusions can be drawn, in particular, for young clusters located at low
Galactic latitudes, the projection effect may have relatively little impact on
these clusters. Therefore, the result of the radial stratification degree in
the short axis direction being greater than that in the long-axis direction
for most young clusters may provide clues to reveal how gas dissipation
affects the morphology of young clusters’ interiors. Since a young cluster
with an age of only a few million years does not have enough time to
completely change its internal structure through dynamical evolution, it is
possible that the internal structure of the cluster may be affected by the gas
dissipation processes, especially in the direction of the short axis of the
cluster interior.

Figure 9: The heights (H) from the Galactic disk versus $q_{all}$-$q_{core}$
of the sample clusters (left panel), as well as the histogram of the
orientation difference ($q_{all}$-$q_{core}$) of them (right panel). The size
of each circle in the left panel indicates the age of a cluster, with its
color displaying the overall ellipticities ($e_{all}$). The vertical dash line
in the left panel represents the Galactic plane, and the horizontal dash line
in this panel indicates the zero stratification degree at the tangential
direction of the cluster interior. The red curve in the right panel is the
gauss fit to the histogram.

We also investigate the tangential stratification in the interior of the
cluster. From the right panel of Figure 9, we can see that the tangential
stratification inside the clusters shows a concentrated distribution. The left
panel of Figure 9 shows that the open clusters with large overall ellipticity
have relatively low tangential stratification, while clusters with small
overall ellipticity near the Galactic plane are mostly in the opposite state.
Older open clusters may exhibit lower stratification in the tangential
direction. However, the number of old clusters is relatively small in the
sample. This trend needs to be further verified by expanding the sample of old
clusters.

The left panel of Figure 9 also shows that the young clusters are mainly
distributed near the disk, and their tangential stratification covers almost
the entire vertical axis. Open clusters still embedded in their natal clouds
are called embedded clusters (Pelupessy & Portegies Zwart, 2012). They evolve
into optically observed open clusters at the cost of losing mass (called
“infant losing weight”; Goodwin & Bastian, 2006). It is important to note that
these young clusters are vulnerable to self-evolution and external forces. As
a result, they tend to exhibit a high degree of stratification. Furthermore,
when approaching the Galactic disk, the external environment becomes more
arcane and complex, and naturally, the degree of stratification should likely
be high. However, the young clusters in the sample are almost uniform in the
distribution of the tangential stratification degree. This needs to be further
investigated in depth. Furthermore, it is intriguing why clusters with large
overall ellipticity have little tangential stratification. Perhaps, this may
be related to the dynamic evolution and the struggle of external forces, the
exact nature of which is also worthy of further investigation.

It is expected that in the future, as the data from the major survey
telescopes become more accurate and homogeneous, we can obtain the real
distribution of star clusters in three-dimensions with the help of high-
precision single-star distances, thus allowing a deeper understanding of the
stratification of open clusters.

Table 2: Grouping Data | logA $>$ 0 | logA $\leq$ 8 | 8 $<$ logA $\leq$ 9 | logA $>$ 9
---|---|---|---|---
$p$ $>$ 1 | 157 | 72 | 74 | 11
$p$ $\leq$ 1 | 108 | 35 | 65 | 8

###  3.3 Influence of Galactic Tides and Differential Rotation on Cluster
Morphology

We use the overall ellipticities of the clusters at two different directions
in the Galactic plane to trace the influence of the Galactic tides and
differential rotation on the cluster morphology. The following coordinate
systems for the location of open clusters are applied: the spherical Galactic
coordinates $(l,b)$; the XYZ coordinates in a Galactic reference frame
centered on the Sun where these position parameters $(X,Y,Z)$ are directly
obtained from Cantat-Gaudin et al. (2018). Since the cluster morphology is
perpendicular to our line of sight, when the sight direction is at $l=90^{o}$
or $270^{o}$, we see their (X, Z) projection (Z-axis from the Galactic center
to North Galactic pole), revealing the full elongation toward the Galactic
center. It is assumed that the elongation of the cluster morphology is mainly
caused by the Galactic tides toward the Galactic center. However, at $l=0^{o}$
or $180^{o}$, we see their (Y, Z) projection, revealing the elongated in the
direction of Galactic differential rotation. The cluster morphology elongated
may be mainly caused by the shear forces that existed at the direction of the
Galactic differential rotation. Therefore, the projection effect is used to
study the influence that the Galactic differential rotation and the Galactic
tides toward the Galactic center exert on the morphology of the clusters.

Figure 10: The distribution of two subsamples elected from the whole sample in
Galactic disk (X-Y plane, top panel) and X-Z plane (bottom panel). The red
open circles in the top and bottom panels indicate the clusters at the
tangential direction ($(80^{o}\leq l\leq 100^{o})$ and $(260^{o}\leq l\leq
280^{o})$) within $|Z|$ $\leq$ 100 pc, and the blue at the radial
($(170^{o}\leq l\leq 190^{o})$ and $(350^{o}\leq l\leq 10^{o})$) within $|Z|$
$\leq$ 100 pc. The sizes of the open circles in the two panels are
proportional to the overall ellipticities of the subsamples. The black dash
line in the bottom panel represents the Galactic disk. The scales of the
circles are shown on the right side of the panel.

We divide the region around the sun into two parts: radial ($(170^{o}\leq
l\leq 190^{o})$ and $(350^{o}\leq l\leq 10^{o})$) and tangential ($(80^{o}\leq
l\leq 100^{o})$ and $(260^{o}\leq l\leq 280^{o})$) parts. And, the clusters in
both parts within $|Z|$ $<$ 100 pc are selected as two subsamples, i.e. the
radial and tangential subsamples. And then the two subsamples are both
displayed in the Galactic plane ($X-Y$ plane) and $X-Z$ plane in Figure 10. It
can be seen from the figure that the overall ellipticities of most of the
clusters in the tangential subsample are larger than that of the clusters in
the radial. We calculate the average of the overall ellipticities of the 265
clusters and that of the two subsamples separately. The average value of the
overall ellipticities of the radial subsample is smaller than that of the
whole sample, while that of the tangential subsample larger than that of the
whole. Then, the ratio of the average differences of overall ellipticities
between the radial subsample and the whole sample to that between the
tangential subsample and the whole is calculated, and it is about -8. This
means that the Galactic tides toward the Galactic center may have more
influence than the shear forces embedded in the Galactic differential rotation
on the overall ellipticities of our sample clusters.

###  3.4 Evolutionary Distribution of the Sample Clusters

Figure 11: The number of the member stars vs. the ages of the sample clusters.
The sizes of the open circles in the panel are proportional to the overall
ellipticities of the sample clusters. The red open triangle in the panel
indicates the cluster (NGC 6791). Two black dash lines are the boundary lines
of different partitions. The two partitions in the left half of the panel are
marked as Region B and Region A from top to bottom, and in the right half are
marked as Region C and Region D respectively. A legend to display the scales
of those open circles is placed on the right side of the panel. Given that the
definition of supercluster masses (massive ($>$ $10^{5}$ $M_{\odot}$) (Emig et
al., 2020)) and assuming that the average mass of every single star is about 1
$M_{\odot}$, the right endpoint of the x-axis in the panel is set to 5.

Figure 11 shows the distribution of the ages and the number of member stars of
the whole sample clusters. The sizes of the open circles in the panel are
proportional to the overall ellipticities of the sample clusters. It is now
generally accepted that the cluster will disintegrate in a few hundred million
years, so we have drawn a horizontal dashed line at the possible upper limit
of the disintegration time (1 Gyr) in this figure. A vertical dashed line is
drawn at the position where the number of cluster members is 1000. The whole
panel in Figure 11 is then divided into four regions, namely Regions A, B, C,
and D. As seen in Figure 11, most of our sample clusters are located in Region
A, which confirms that open clusters are a class of stellar systems with a
certain age range and a range of the number of member stars.

Open clusters lose their member stars as they evolve. Therefore, assuming
different loss rates of member stars during the cluster evolution, the
clusters seem to exhibit a trend of gradually increasing overall ellipticities
from the lower right of Region A to the upper left of Region A, as shown in
Figure 11. This may suggest that the gravitational binding of the cluster
decreases during cluster evolution with the losing of the cluster members, and
the cluster morphology is susceptible to external disturbances, thus showing a
tendency for the overall ellipticities to become larger. We should note that
the loss rate of cluster member stars does not necessarily remain constant
throughout the evolution process. As the number of cluster members decreases
and the surrounding environment changes during evolution, the loss rate may
also change. The loss rate of member stars in cluster evolution needs to be
studied in depth in the future.

The relatively low number of clusters in region B is consistent with the
finding of a lack of older clusters with ages $t$ $>$ 1 $Gyr$ by Piskunov et
al. (2018). It is interesting to note the lack of clusters in the lower-left
corner of Region A, which could be the result of infant mortality (e.g. Lada &
Lada, 2003; Bastian et al., 2005; Fall et al., 2005; Goodwin & Bastian, 2006).
Previous studies have shown that most stars born in clusters do not remain
there past the age of 10 million years, and only about 10% of stars are
observed in gravitationally bound clusters at this age (Lada & Lada, 2003).
This may be one possible reason for the lack of clusters in the lower-left
corner of Region A.

There are two old and rich clusters in Region C. One of them, NGC 6791, is the
most populous and oldest cluster in our sample with a relatively small overall
ellipticity. Dalessandro et al. (2015) estimated the initial cluster mass of
NGC 6791 to be $(1.5-4)$ $\times$ $10^{5}$ $M_{\odot}$ based on analytical
studies and N-body simulations. It is very unlikely that such an old cluster
will maintain a very low member loss rate during its long-term evolution,
suggesting the initial cluster of it should locate in the right lower part of
Region D. Therefore, the superclusters, massive ($>$ $10^{5}$ $M_{\odot}$) and
compact clusters (Emig et al., 2020), may give a possible evolutionary origin
for NGC 6791.

##  4 Summary

In this study, we use the member catalog of open clusters available in the
literature (Cantat-Gaudin et al., 2018) to give the shape parameters of 265
open clusters. Meanwhile, we conclude in the following for the analysis of
cluster morphology evolution by the shape parameters:

1\. The overall shapes of the sample clusters become more elliptical as they
grow older, while their core shapes remains circular or slightly trend to
circularize. The orientations of the most sample clusters are within $|q|$ $<$
45o, which means that the sample clusters are more elongated in a direction
parallel to the Galactic plane than perpendicular to it.

2\. There is a negative correlation between the ellipticities and the number
of member stars of the sample clusters, implying that the higher the number of
member stars of the clusters, the stronger the cluster’s gravitational binding
to resist the external disturbances, and the more it can keep smaller
ellipticity.

3\. No significant correlations between the inner ellipticities and masses and
between the overall ellipticities and masses are detected for the clusters
with log(age/year) $<$ 8\. However, negative correlations of the inner and
overall ellipticities with masses are detected for the clusters with
log(age/year) $\geq$ 8, especially more significant between the overall
ellipticities and the masses, suggesting that the overall shapes of the
clusters may be influenced by the cluster’s mass, in addition to other
factors.

4\. For sample clusters with a high number of member stars, the radial
stratification degree in the short axis direction is almost the same as that
in the long. Both within the clusters in the short-axis direction and long-
axis direction, most sample clusters with fewer members have a larger radial
stratification range than those with more members. The number of clusters in
the sample with a greater radial stratification in the short-axis direction
than in the long-axis direction is almost the same as the opposite case. For
most young sample clusters, the radial stratification degree of the short axis
direction is greater than that of the long, implying that the degree of
stratification in both directions within the young sample cluster might be
unevenly affected by gas dissipation processes. The tangential stratification
inside the clusters of our sample shows a concentrated distribution. The
sample clusters with relatively large overall ellipticity have relatively low
tangential stratification, while the clusters with smaller overall ellipticity
near the Galactic plane are mostly in the opposite state. Older sample
clusters may exhibit lower stratification in the tangential direction,
implying that those clusters may continue to survive for a long time at a low
level of stratification.

5\. Our analysis of two subsamples of clusters distribution on the radial and
tangential directions shows that the overall ellipticities of sample clusters
may be more susceptible to the influence of Galactic tides toward the Galactic
center than the shear forces embedded in Galactic differential rotation.

6\. Sample clusters appear to show a general trend of increasing overall
ellipticity from the lower right of Region A to the upper left of Region A,
which means the loss of cluster members drives the overall cluster morphology
to become elliptical. By analyzing the distribution of the ages and number of
members of star clusters, we suggest that NGC 6791 may originate from
superclusters.

We acknowledge that our sample clusters are still incomplete. A large and
complete sample is needed to obtain more reliable conclusions about the shape
evolution of open clusters. Besides, It is expected that in the future, as the
data from the major survey telescopes become more accurate and homogeneous, we
can obtain the real distribution of star clusters in three-dimensions, thus
allowing a deeper understanding of the morphology evolution of open clusters.

This work was supported by the National Natural Science Foundation of China
under grants U2031209 and U2031204, the Youth Innovation Promotion Association
CAS (grant No.2018080), 2017 Heaven Lake Hundred-Talent Program of Xinjiang
Uygur Autonomous Region of China. The Authors thank to the reviewer for the
very helpful comments and suggestions. We are grateful to Prof. Dr. Wen-Ping
Chen for an in-depth discussion and suggestions on the collaboration. We would
also like to thank Ms. Chunli Feng for touching up the language of the
article. This study has made use of the Gaia DR2, operated at the European
Space Agency (ESA) space mission (Gaia). The Gaia archive website is
https://archives.esac.esa.int/gaia/.

## References

  * Adams (2000) Adams, F. C. 2000, ApJ, 542, 964. doi:10.1086/317052
  * Alves & Bouy (2012) Alves, J. & Bouy, H. 2012, A&A, 547, A97. doi:10.1051/0004-6361/201220119
  * Arenou et al. (2018) Arenou, F., Luri, X., Babusiaux, C., et al. 2018, A&A, 616, A17. doi:10.1051/0004-6361/201833234
  * Balaguer-Núñez et al. (2020) Balaguer-Núñez, L., López del Fresno, M., Solano, E., et al. 2020, MNRAS, 492, 5811. doi:10.1093/mnras/stz3610
  * Bastian et al. (2005) Bastian, N., Gieles, M., Lamers, H. J. G. L. M., et al. 2005, A&A, 431, 905. doi:10.1051/0004-6361:20041078
  * Bergond et al. (2001) Bergond, G., Leon, S., & Guibert, J. 2001, A&A, 377, 462. doi:10.1051/0004-6361:20011043
  * Bhattacharya et al. (2017) Bhattacharya, S., Mishra, I., Vaidya, K., et al. 2017, ApJ, 847, 138. doi:10.3847/1538-4357/aa89e2
  * Bossini et al. (2019) Bossini, D., Vallenari, A., Bragaglia, A., et al. 2019, A&A, 623, A108. doi:10.1051/0004-6361/201834693
  * Bonatto et al. (2006) Bonatto, C., Kerber, L. O., Bica, E., et al. 2006, A&A, 446, 121. doi:10.1051/0004-6361:20053573
  * Brinkmann et al. (2017) Brinkmann, N., Banerjee, S., Motwani, B., et al. 2017, A&A, 600, A49. doi:10.1051/0004-6361/201629312
  * Cantat-Gaudin et al. (2018) Cantat-Gaudin, T., Jordi, C., Vallenari, A., et al. 2018, A&A, 618, A93. doi:10.1051/0004-6361/201833476
  * Chen et al. (2004) Chen, W. P., Chen, C. W., & Shu, C. G. 2004, AJ, 128, 2306. doi:10.1086/424855
  * Dalessandro et al. (2015) Dalessandro, E., Miocchi, P., Carraro, G., et al. 2015, MNRAS, 449, 1811. doi:10.1093/mnras/stv395
  * Elmegreen (1983) Elmegreen, B. G. 1983, MNRAS, 203, 1011. doi:10.1093/mnras/203.4.1011
  * Emig et al. (2020) Emig, K. L., Bolatto, A. D., Leroy, A. K., et al. 2020, ApJ, 903, 50. doi:10.3847/1538-4357/abb67d
  * Evans et al. (2018) Evans, D. W., Riello, M., De Angeli, F., et al. 2018, A&A, 616, A4. doi:10.1051/0004-6361/201832756
  * Fall et al. (2005) Fall, S. M., Chandar, R., & Whitmore, B. C. 2005, ApJ, 631, L133. doi:10.1086/496878
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1. doi:10.1051/0004-6361/201833051
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1. doi:10.1051/0004-6361/201629272
  * Getman et al. (2018) Getman, K. V., Kuhn, M. A., Feigelson, E. D., et al. 2018, MNRAS, 477, 298. doi:10.1093/mnras/sty473
  * Gieles & Baumgardt (2008) Gieles, M. & Baumgardt, H. 2008, MNRAS, 389, L28. doi:10.1111/j.1745-3933.2008.00515.x
  * Goodwin & Bastian (2006) Goodwin, S. P. & Bastian, N. 2006, MNRAS, 373, 752.doi:10.1111/j.1365-2966.2006.11078.x
  * Hills (1980) Hills, J. G. 1980, ApJ, 235, 986. doi:10.1086/157703
  * Jeans (1916) Jeans, J. H. 1916, MNRAS, 76, 567. doi:10.1093/mnras/76.7.567
  * Kharchenko et al. (2009) Kharchenko, N. V., Berczik, P., Petrov, M. I., et al. 2009, A&A, 495, 807. doi:10.1051/0004-6361/200810407
  * Krone-Martins & Moitinho (2014) Krone-Martins, A. & Moitinho, A. 2014, A&A, 561, A57. doi:10.1051/0004-6361/201321143
  * Kroupa et al. (2001) Kroupa, P., Aarseth, S., & Hurley, J. 2001, MNRAS, 321, 699. doi:10.1046/j.1365-8711.2001.04050.x
  * Kroupa (2002) Kroupa, P. 2002, MNRAS, 330, 707. doi:10.1046/j.1365-8711.2002.05128.x
  * Lada & Lada (2003) Lada, C. J. & Lada, E. A. 2003, ARA&A, 41, 57. doi:10.1146/annurev.astro.41.011802.094844
  * Lindegren et al. (2018) Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2. doi:10.1051/0004-6361/201832727
  * Mamon et al. (2019) Mamon, G. A., Cava, A., Biviano, A., et al. 2019, A&A, 631, A131. doi:10.1051/0004-6361/201935081
  * Mathieu (1983) Mathieu, R. D. 1983, ApJ, 267, L97. doi:10.1086/184011
  * Megeath et al. (2016) Megeath, S. T., Gutermuth, R., Muzerolle, J., et al. 2016, AJ, 151, 5. doi:10.3847/0004-6256/151/1/5
  * Moeckel & Bate (2010) Moeckel, N. & Bate, M. R. 2010, MNRAS, 404, 721. doi:10.1111/j.1365-2966.2010.16347.x
  * Nilakshi et al. (2002) Nilakshi, Sagar, R., Pandey, A. K., et al. 2002, A&A, 383, 153. doi:10.1051/0004-6361:20011719
  * Oort (1979) Oort, J. H. 1979, A&A, 78, 312
  * Pandey & Mahra (1986) Pandey, A. K. & Mahra, H. S. 1986, Ap&SS, 126, 167. doi:10.1007/BF00644184
  * Pandey et al. (1990) Pandey, A. K., Mahra, H. S., & Sagar, R. 1990, AJ, 99, 617. doi:10.1086/115355
  * Pelupessy & Portegies Zwart (2012) Pelupessy, F. I. & Portegies Zwart, S. 2012, MNRAS, 420, 1503. doi:10.1111/j.1365-2966.2011.20137.x
  * Piskunov et al. (2007) Piskunov, A. E., Schilbach, E., Kharchenko, N. V., et al. 2007, A&A, 468, 151. doi:10.1051/0004-6361:20077073
  * Piskunov et al. (2008) Piskunov, A. E., Schilbach, E., Kharchenko, N. V., et al. 2008, A&A, 477, 165. doi:10.1051/0004-6361:20078525
  * Piskunov et al. (2018) Piskunov, A. E., Just, A., Kharchenko, N. V., et al. 2018, A&A, 614, A22. doi:10.1051/0004-6361/201732337
  * Raboud & Mermilliod (1998a) Raboud, D. & Mermilliod, J.-C. 1998, A&A, 329, 101
  * Raboud & Mermilliod (1998b) Raboud, D. & Mermilliod, J.-C. 1998, A&A, 333, 897
  * Ryon et al. (2017) Ryon, J. E., Gallagher, J. S., Smith, L. J., et al. 2017, ApJ, 841, 92. doi:10.3847/1538-4357/aa719e
  * Tang et al. (2019) Tang, S.-Y., Pang, X., Yuan, Z., et al. 2019, ApJ, 877, 12. doi:10.3847/1538-4357/ab13b0
  * Tutukov (1978) Tutukov, A. V. 1978, A&A, 70, 57
  * Wielen (1971) Wielen, R. 1971, A&A, 13, 309
  * Yang et al. (2013) Yang, S.-C., Sarajedini, A., Deliyannis, C. P., et al. 2013, ApJ, 762, 3. doi:10.1088/0004-637X/762/1/3
  * Zhai et al. (2017) Zhai, M., Abt, H., Zhao, G., et al. 2017, AJ, 153, 57. doi:10.3847/1538-3881/153/2/57
  * Zhang et al. (2020) Zhang, Y., Tang, S.-Y., Chen, W. P., et al. 2020, ApJ, 889, 99. doi:10.3847/1538-4357/ab63d4
  * Zhong et al. (2019) Zhong, J., Chen, L., Kouwenhoven, M. B. N., et al. 2019, A&A, 624, A34. doi:10.1051/0004-6361/201834334