# The Classification of all weak solutions to $-\Delta u={u^{-\gamma}}$ in the
half-space

Luigi Montoro  ,  Luigi Muglia  and  Berardino Sciunzi ∗

###### Abstract.

We provide the classification of all the positive solutions to $-\Delta
u=\frac{1}{u^{\gamma}}$ in the half space, under minimal assumption.

###### Key words and phrases:

Singular solutions, half-space, classification result.

2020 Mathematics Subject Classification: 35J75, 35A02, 35B09

∗ corresponding author

##  1\. Introduction

We deal with the classification of positive weak solutions to the singular
problem

$\begin{cases}\displaystyle-\Delta
u=\frac{1}{u^{\gamma}}&\text{in}\,\,\mathbb{R}^{N}_{+}\\\
u>0&\text{in}\,\,\mathbb{R}^{N}_{+}\\\
u=0&\text{on}\,\,\partial\mathbb{R}^{N}_{+},\end{cases}$
($\mathcal{P}_{\gamma}$)

where $N\geq 1$, $\gamma>0$ and
$\mathbb{R}^{N}_{+}:=\\{x\in\mathbb{R}^{N}:x_{N}>0\\}$.

It is well known that the solutions, in general, have not $H^{1}$ regularity
up to the boundary (see [12]) and the equation has to be understood in the
following way: $u\in H^{1}_{loc}(\mathbb{R}^{N}_{+})$ and

$\int_{\mathbb{R}^{N}_{+}}(\nabla
u,\nabla\varphi)=\int_{\mathbb{R}^{N}_{+}}\frac{1}{u^{\gamma}}\varphi\qquad\forall\varphi\in
C^{1}_{c}(\mathbb{R}_{+}^{N}).$ (1)

The zero Dirichlet boundary condition also has to be understood in the weak
meaning

$(u-\varepsilon)^{+}\varphi\chi_{\mathbb{R}^{N}_{+}}\in
H^{1}_{0}(\mathbb{R}^{N}_{+})\qquad\forall\,\varphi\in
C^{1}_{c}(\mathbb{R}^{N})\,.$

######  Remark 1.

In some cases, see e.g. [3] for the case of bounded domains, when dealing with
singular equations, it is natural to understand the zero Dirichlet datum
requiring that $u^{\frac{\gamma+1}{2}}\in H^{1}_{0}(\Omega)$. We stress the
fact that such assumption already implies that $(u-\varepsilon)^{+}\in
H^{1}_{0}(\Omega)$ so that our analysis do apply also to this case.

The aim of this paper is to improve our previous results in [13] where we
provided the classification of the solutions in the case $\gamma>1$ assuming a
uniform $L^{\infty}$ bound in the strips
$\\{x\in\mathbb{R}^{N}_{+}\,:\,0<x_{N}<b\\}$. We remove such an assumption as
well as we remove any regularity assumption on the solutions. Furthermore we
deduce that, in the case $0<\gamma\leq 1$ there are no solutions.

Our main result is the following

######  Theorem 2.

Let $\gamma>1$ and let $u$ be a solution to ($\mathcal{P}_{\gamma}$). Then

$u(x)=u(x_{N})$

with

$u(x_{N})=\frac{(\gamma+1)^{\frac{2}{\gamma+1}}}{(2\gamma-2)^{\frac{1}{\gamma+1}}}(x_{N})^{\frac{2}{\gamma+1}}$

or

$u(x_{N})=\lambda^{-\frac{2}{\gamma+1}}v(\lambda x_{N})\quad\lambda>0,$

where $v(t)\in C^{2}(\mathbb{R}_{+})\cap C(\overline{\mathbb{R}_{+}})$ is the
unique positive solution to

$\displaystyle\begin{cases}\displaystyle-v^{\prime\prime}=v^{-\gamma}\quad&t>0\\\
v(0)=0\quad\displaystyle&\lim_{t\to+\infty}v^{\prime}(t)=1.\end{cases}$ (2)

If else $0<\gamma\leq 1$, then ($\mathcal{P}_{\gamma}$) has no solutions.

######  Remark 3.

Theorem 2 concludes the picture for the classification of solutions to
semilinear elliptic equations in half-spaces involving power nonlinearities.
Actually it is well known that, in the case of positive powers ($\gamma<0$
with our notation), solutions do not exist up to a critical value. The case of
negative powers ($\gamma>0$ with our notation) still exhibit non existence of
the solutions for the case $0<\gamma\leq 1$. On the contrary, if $\gamma>1$
solutions exist and are classified by our Theorem 2.

The study of semilinear elliptic problems involving singular nonlinearities
has attracted many studies starting from the seminal papers [7], see also [1,
2, 3, 4, 5, 8, 10, 12, 14] and the references therein. Starting from general
equations of the form $-\Delta u=u^{-\gamma}+f(u)$, the pure singular problem
$-\Delta u=u^{-\gamma}$ in half spaces arises e.g. when performing a blow-up
analysis near the boundary of the domain since the term ${1}/{u^{\gamma}}$ is
the leading part near zero, see [6, 9].

Our proof is based on the following steps:

Step 1: we start proving the local boundedness of the solutions in
$\overline{\mathbb{R}_{+}}$. We follow the clever technique of J. Serrin [15]
to carry out the Moser iteration scheme. Some care is needed since solutions
have not $H^{1}$-regularity up to the boundary.

Step 2: we prove that the solutions are continuous up to the boundary. Thanks
to the local boundedness of the solutions we succeed in the use a barriers
technique to show that, near the boundary, the solution is trapped among to
continuous functions. This is reduced to a sliding method based on comparison
arguments. This technique is generally successfully used in this issue but the
case of the half space is much harder to study because is not possible to
construct directly global barriers all over the domain, unless a-priori bounds
are assumed.

Step 3: we prove that the solution has 1-D symmetry. Here we develop a
completely different proof that the one in [13], based on the use of the
Kelvin transformation. Again the lack of regularity of the solutions up to the
boundary is the crucial issue.

Conclusion: the proof of the main result is consequently reduce to a ODE
analysis.

##  2\. Regularity of the solutions

We start this section studying the local boundedness of the solutions. Let us
consider $H^{1}_{loc}(\mathbb{R}^{N}_{+})$ weak solutions to the problem that
fulfil the Dirichlet datum in the generalized meaning so that

$(u-\varepsilon)^{+}\in H^{1}_{loc}(\overline{\mathbb{R}^{N}_{+}})\,.$

The proof that follows is based on a standard use of the Moser iteration
scheme. There is here an added difficulty given by the singular nature of the
problem that we overcome following closely the technique developed by J.
Serrin in [15].

######  Theorem 4.

If $u\in H^{1}_{loc}(\mathbb{R}^{N}_{+})$ is a weak solution problem
$\mathcal{P}_{\gamma}$, then $u$ is locally bounded.

###### Proof.

For $\varepsilon>0$ fixed, set

$u_{\varepsilon}\,:=\,(u-\varepsilon)^{+}\,+\,\frac{1}{\varepsilon^{\gamma}}\,.$

Furthermore, for $q\geq 1$ and $l>\frac{1}{\varepsilon^{\gamma}}$, we define

$F(u_{\varepsilon})\,:=\,\begin{cases}u_{\varepsilon}^{q}\qquad\qquad\qquad\qquad\qquad
if\quad\frac{1}{\varepsilon^{\gamma}}\leq u_{\varepsilon}\leq l\\\
ql^{q-1}u_{\varepsilon}-(q-1)l^{q}\qquad\quad if\quad u_{\varepsilon}\geq
l\end{cases}$ (3)

and, consequently,

$G(u)\,:=\,F(u_{\varepsilon})F^{\prime}(u_{\varepsilon})-\frac{q}{\varepsilon^{\beta\gamma}}$
(4)

where $\beta$ is such that $2q=\beta+1$. We now indicate by $\eta$ a standard
smooth cut-off function and consider the test function

$\varphi\,:=\,G(u)\eta^{2}\,.$

It is easy to check that this is a suitable test function for our problem.
Note in fact that $(u-\varepsilon)^{+}\in
H^{1}_{loc}(\overline{\mathbb{R}^{N}_{+}})$ by the definition of the Dirichlet
boundary datum and $G(u)$ is linear for $u_{\varepsilon}>l$. Consequently we
deduce that

$\int_{\mathbb{R}^{N}_{+}}\,|\nabla
u|^{2}G^{\prime}(u)\eta^{2}=\int_{\mathbb{R}^{N}_{+}}-2G(u)\eta(\nabla
u\,,\nabla\eta)+\frac{G(u)\eta^{2}}{u^{\gamma}}$

so that

$\begin{split}\int_{\mathbb{R}^{N}_{+}}\,|\nabla
u|^{2}G^{\prime}(u)\eta^{2}&\leq 2\int_{\mathbb{R}^{N}_{+}}G(u)\eta|\nabla
u||\nabla\eta|+\frac{G(u)\eta^{2}}{u^{\gamma}}\\\ &\leq
2\int_{\mathbb{R}^{N}_{+}}G(u)\eta|\nabla
u||\nabla\eta|+\int_{\mathbb{R}^{N}_{+}}\frac{G(u)\eta^{2}}{\varepsilon^{\gamma}}\\\
&\leq 2\int_{\mathbb{R}^{N}_{+}}G(u)\eta|\nabla
u||\nabla\eta|+\int_{\mathbb{R}^{N}_{+}}G(u)\eta^{2}u_{\varepsilon}\,.\\\
\end{split}$

Taking into account the definition of $F(\cdot)$ and $G(\cdot)$, it is easily
seen that

$G^{\prime}(u)\,:=\,\begin{cases}\displaystyle\frac{\beta}{q}\,(F^{\prime})^{2}\qquad\,\,\,\,\,if\quad(u-\varepsilon)<l-\frac{1}{\varepsilon^{\gamma}}\\\
\\\ \,\,\,\,\displaystyle(F^{\prime})^{2}\qquad\quad
if\quad(u-\varepsilon)>l-\frac{1}{\varepsilon^{\gamma}}\end{cases}$

and $G\leq FF^{\prime}$. Thus, setting

$v(x)\,:=\,F(u_{\varepsilon})$

and observing that $u_{\varepsilon}\,F^{\prime}\leq qF$, we get that

$\begin{split}\min\left\\{1,\frac{\beta}{q}\right\\}\int_{\mathbb{R}^{N}_{+}}|\nabla
v|^{2}\eta^{2}&=\min\left\\{1,\frac{\beta}{q}\right\\}\int_{\mathbb{R}^{N}_{+}}|\nabla
u|^{2}(F^{\prime}(u_{\varepsilon}))^{2}\eta^{2}\\\
&\leq\int_{\mathbb{R}^{N}_{+}}|\nabla u|^{2}G^{\prime}(u)\eta^{2}\\\ &\leq
2\int_{\mathbb{R}^{N}_{+}}FF^{\prime}\eta|\nabla
u||\nabla\eta|+q\int_{\mathbb{R}^{N}_{+}}v^{2}\eta^{2}\,.\\\ \end{split}$

Recalling that $\nabla v=F^{\prime}(u_{\varepsilon})\nabla u$, by a simple use
of the Young inequality, we deduce that

$\int_{\mathbb{R}^{N}_{+}}|\nabla v|^{2}\eta^{2}\leq
C\int_{\mathbb{R}^{N}_{+}}|\nabla\eta|^{2}v^{2}+C\,q\int_{\mathbb{R}^{N}_{+}}v^{2}\eta^{2}\,.$
(5)

The Sobolev inequality now provides that

$\left(\int_{\mathbb{R}^{N}_{+}}(v\eta)^{2^{*}}\right)^{\frac{1}{2^{*}}}\leq\left(\int_{\mathbb{R}^{N}_{+}}|\nabla
v|^{2}\eta^{2}+|\nabla\eta|^{2}v^{2}+2|\nabla\eta||\nabla v|\eta
v\right)^{\frac{1}{2}}$

so that, using again the Young inequality and (5), we get

$\left(\int_{\mathbb{R}^{N}_{+}}(v\eta)^{2^{*}}\right)^{\frac{1}{2^{*}}}\leq
Cq\,\left(\int_{\mathbb{R}^{N}_{+}}v^{2}|\nabla\eta|^{2}\right)^{\frac{1}{2}}+Cq\left(\int_{\mathbb{R}^{N}_{+}}v^{2}\eta^{2}\right)^{\frac{1}{2}}\,.$

This estimate recover the ones in the proof of Theorem 1 of [15] (see equation
(19)). Consequently the remaining part of the proof can be carried over
repeating verbatim the argument of Serrin, thus proving the local boundedness
of $u_{\varepsilon}$ that implies the local boundedness of $u$. ∎

Standard elliptic estimates (see [11]) show that the solution is smooth in the
interior of the half space. Furthermore, we are now in position to prove that
any solution is necessarily continuous up to the boundary. Note that this is
the best regularity we may hope for the case of singular equations.

######  Theorem 5.

If $u\in H^{1}_{loc}(\mathbb{R}^{N}_{+})$ is a weak solution problem
$\mathcal{P}_{\gamma}$, then $u\in C^{0}(\overline{\mathbb{R}^{N}_{+}})$.

###### Proof.

The continuity of the solution in the interior of the domain follows by
standard regularity estimates, since the right hand side is bounded there (see
e.g. [11]). From Lemma 4 in [13] we now that there exists a constant
$C=C(\gamma)$ such that

$u\geq Cx_{N}^{\frac{2}{\gamma+1}}\quad\text{in}\,\,\mathbb{R}^{N}_{+}.$

We therefore turn to provide a bound from above. Here the local boundedness of
the solution comes into play. Let us consider $w_{x_{0}}=w\in
C^{2}(B_{3}(x_{0})\setminus\overline{B_{1}(x_{0})})\cap
C^{0}(\overline{B_{3}(x_{0})\setminus\overline{B_{1}(x_{0})}})$ the unique
solution to

$\begin{cases}-\Delta
w=\frac{1}{w^{\gamma}}&\text{in}\,\,B_{3}(x_{0})\setminus\overline{B_{1}(x_{0})}\\\
w>0&\text{in}\,\,B_{3}(x_{0})\setminus\overline{B_{1}(x_{0})}\\\
w=0&\text{on}\,\,\partial\left(B_{3}(x_{0})\setminus\overline{B_{1}(x_{0})}\right).\end{cases}$
(6)

Setting for $k>1$

$w_{k}\,:=\,kw$

it is easy to see that $w_{k}$ is a supersolution to (6). Furthermore

$w_{k}\geq k\underset{\partial
B_{2}(x_{0})}{\min}\,w\qquad\text{on}\quad\partial B_{2}(x_{0})\,.$

We now we adopt the notation $x=(x^{\prime},x_{N})$ (with
$x^{\prime}\in\mathbb{R}^{N-1}$) and we assume, with no loss of generality,
that

$x_{0}=x_{0}(t)\,:=\,(0^{\prime},t)$

with $-2<t<-1$ and set (see Figure (2))

$\Omega_{t}\,:=\,\left(B_{2}(x_{0})\setminus\overline{B_{1}(x_{0})}\right)\cap\mathbb{R}^{N}_{+}$

Exploiting the fact that the solution is locally bounded, see Theorem 4, we
can assume that

$u\leq w_{k}\quad\text{on}\,\,\,\partial\Omega_{t}\qquad\text{for any
}\,\,\,-2<t<-1\,.$

We claim that: $u\leq w_{k}$ in $\Omega_{t}$ for any $-2<t<-1$.
To prove this, let us first note that

$(u-\varepsilon-w_{k})^{+}\chi_{\Omega_{t}}\in H^{1}_{0}(\Omega_{t})\,.$

This is a consequence of the fact that $(u-\varepsilon)^{+}\in
H^{1}_{loc}(\overline{\mathbb{R}^{N}_{+}})$ by the weak definition of the
Dirichlet datum and, for $t<-1$, it follows that $w_{k}\in H^{1}(\Omega_{t})$
by construction. Consequently there exist a sequence of positive functions
$\varphi_{n}\in C^{1}_{c}(\Omega_{t})$ converging to $(u-\varepsilon-
w_{k})^{+}$ in $H^{1}_{0}(\Omega_{t})$. We set

$\psi_{n}\,:=\,\min\\{(u-\varepsilon-w_{k})^{+},\varphi_{n}\\}$

so that the support of $\psi_{n}$ is strictly contained in $\Omega_{t}$ and in
the region where $u\geq w_{k}$. Plugging $\psi_{n}$ as test function in the
problem we deduce that

$\displaystyle\int_{\Omega_{t}}(\nabla(u-\varepsilon-
w_{k})^{+},\nabla\psi_{n})$ $\displaystyle=$
$\displaystyle\int_{\Omega_{t}}(\nabla u-\nabla w_{k},\nabla\psi_{n})$
$\displaystyle=$
$\displaystyle\int_{\Omega_{t}}\left(\frac{1}{u^{\gamma}}-\frac{1}{w_{k}^{\gamma}}\right)\psi_{n}\leq
0\,.$

We can therefore pass to the limit to deduce that

$\int_{\Omega_{t}}|\nabla(u-\varepsilon-w_{k})^{+}|^{2}=0$

thus proving that $u\leq w_{k}+\varepsilon$ in $\Omega_{t}$. The fact that
$\varepsilon$ is arbitrarily chosen implies that $u\leq w_{k}$ in $\Omega_{t}$
for any $t<-1$. and the claim is proved.
By continuity it follows now that

$u\leq w_{k}\quad\text{in}\,\,\Omega_{-1}\,$

By translation invariance we have that the centre $x_{0}$ can be arbitrary
chosen with the $x_{N}$-coordinate equal to $-1$. Collecting this and the
bound from above, it follows that the graph of $u$ is trapped among the graph
of two functions that are continuous up to the boundary. Namely $u$ is
continuous up to the boundary. ∎

##  3\.  1-D symmetry

We start this section recalling the Kelvin transformation of the solutions $u$
given by:

$\hat{u}(x)\,:=\,\frac{1}{|x|^{N-2}}u(\frac{x}{|x|^{2}})\qquad x\neq 0\,.$

It is known that

$\Delta\hat{u}=\frac{1}{|x|^{N+2}}\Delta u|_{\frac{x}{|x|^{2}}}\qquad x\neq
0\,.$

so that

$-\Delta\hat{u}\,=\,\frac{1}{|x|^{\beta_{\gamma}}}\frac{1}{\hat{u}^{\gamma}}\qquad
x\neq 0\qquad\text{for}\,\,\,\beta_{\gamma}=N+2+\gamma(N-2)\,.$

Such equation is classical in the interior of the half space since $u$ is
smooth there. On the contrary the transformation produces a possible
singularity at zero and, in any case, $\hat{u}$ has not $H^{1}$ regularity up
to the boundary as well as the starting solution $u$. Consequently the
equation is still understood as follows:

$\int_{\mathbb{R}^{N}_{+}}\nabla\hat{u}\cdot\nabla\varphi=\int_{\mathbb{R}^{N}_{+}}\frac{1}{|x|^{\beta_{\gamma}}}\frac{1}{\hat{u}^{\gamma}}\varphi\qquad\forall\varphi\in
C^{1}_{c}(\mathbb{R}_{+}^{N}).$ (7)

We are ready to prove

######  Theorem 6.

If $u\in H^{1}_{loc}(\mathbb{R}^{N}_{+})$ is a weak solution problem
$\mathcal{P}_{\gamma}$, then $u$ has 1-D symmetry, namely

$u=u(x_{N}).$

###### Proof.

By Theorem 4 and Theorem 5 we know that $u$ is locally bounded and continuous
up to the boundary. Furthermore $u$ is smooth in the interior of the half
space by standard elliptic estimates. For any $\lambda<0$, we set

$\Sigma_{\lambda}\,:=\,\\{x\in\mathbb{R}^{N}_{+}\,:\,x_{1}<\lambda\\}\,$

adopting here the notation $x=(x_{1},x^{\prime})$ with
$x^{\prime}\in\mathbb{R}^{N-1}$. Furthermore we set

$x_{\lambda}=(2\lambda-x_{1},x^{\prime})$
$\hat{u}_{\lambda}(x)=\hat{u}(x_{\lambda})=\hat{u}(2\lambda-
x_{1},x^{\prime})\quad\text{in}\quad\Sigma_{\lambda}\,.$

By the reflection invariance of the Laplace operator, we deduce that

$\int_{\mathbb{R}^{N}_{+}}\nabla\hat{u}_{\lambda}\cdot\nabla\varphi=\int_{\mathbb{R}^{N}_{+}}\frac{1}{|x_{\lambda}|^{\beta_{\gamma}}}\frac{1}{\hat{u}_{\lambda}^{\gamma}}\varphi\qquad\forall\varphi\in
C^{1}_{c}(\mathbb{R}_{+}^{N}).$ (8)

We stress that the transformation produces a possible singularity for $u$ at
the origin and, consequently, a possible singularity for $\hat{u}$ at
$0_{\lambda}$ the reflected point of the origin, namely
$0_{\lambda}=(-2\lambda,x^{\prime})$. Now we consider the cut-off functions
$0\leq\phi_{R},\varphi_{\varepsilon}\leq 1$ such that

$\begin{cases}\phi_{R}=1&\text{in}\,\,B_{R}(0)\\\
\phi_{R}=0&\text{in}\,\,\mathbb{R}^{N}\setminus B_{2R}(0)\\\
\displaystyle|\nabla\phi_{R}|\leq\frac{C}{R},&\text{in}\,\,B_{2R}(0)\setminus
B_{R}(0)\,.\end{cases}$ (9)

In dimension $N=2$ we will exploit the fact that, by a capacitary argument,
namely considering a suitable truncation of the fundamental solution, it is
possible to assume that

$\int_{B_{R^{2}}(0)\setminus
B_{R}(0)}|\nabla\phi_{R}|^{2}\underset{R\rightarrow\infty}{\longrightarrow}\,0\,.$

For the cut-off function $\varphi_{\varepsilon}$ we assume that:

$\varphi_{\varepsilon}=\begin{cases}0&\text{in}\,\,B_{\varepsilon^{2}}(0_{\lambda})\\\
1&\text{in}\,\,\mathbb{R}^{N}\setminus B_{\varepsilon}(0_{\lambda})\\\
\end{cases}$ (10)

Again by a capacitary argument (and here in any dimension) we may and do
assume that

$\int_{B_{\varepsilon}(0_{\lambda})\setminus
B_{\varepsilon^{2}}(0_{\lambda})}|\nabla\varphi_{\varepsilon}|^{2}\,\underset{\varepsilon\rightarrow
0}{\longrightarrow}\,0\,.$

Given $\tau>0$ we now consider

$w=w_{\varepsilon,\tau,R,\lambda}\,:=\,(\hat{u}-\hat{u}_{\lambda}-\tau)^{+}\varphi_{\varepsilon}^{2}\phi_{R}^{2}\,.$

Since, by Theorem 5 $u$ is continuous up to the boundary, then $\hat{u}$ is
continuous up to the boundary too, for $x\neq 0$. Consequently, recalling the
construction of the cut-off functions, $w$ is a $H^{1}$-function with its
support that is compactly contained in the half space. Consequently, by
standard density arguments, we can use it as test function in (7) and in (8).
Thus we get that

$\displaystyle\int_{\Sigma_{\lambda}}(\nabla\hat{u},\nabla
w)=\int_{\Sigma_{\lambda}}\frac{1}{|x|^{\beta_{\gamma}}}\frac{1}{\hat{u}^{\gamma}}w$
$\displaystyle\int_{\Sigma_{\lambda}}(\nabla\hat{u}_{\lambda},\nabla
w)=\int_{\Sigma_{\lambda}}\frac{1}{|x_{\lambda}|^{\beta_{\gamma}}}\frac{1}{\hat{u}_{\lambda}^{\gamma}}w$

so that

$\displaystyle\int_{\Sigma_{\lambda}}(\nabla(\hat{u}-\hat{u}_{\lambda}),\nabla
w)=\int_{\Sigma_{\lambda}}\left(\frac{1}{|x|^{\beta_{\gamma}}}\frac{1}{\hat{u}^{\gamma}}-\frac{1}{|x_{\lambda}|^{\beta_{\gamma}}}\frac{1}{\hat{u}_{\lambda}^{\gamma}}\right)w\leq
0\,,$

since $\hat{u}\geq\hat{u}_{\lambda}$ on the support of $w$ and
$|x_{\lambda}|\leq|x|$ in $\Sigma_{\lambda}$.

We now set

$z_{\tau}\,:=\,(\hat{u}-\hat{u}_{\lambda}-\tau)^{+}$

and, by direct computation, we deduce that

$\int_{\Sigma_{\lambda}}|\nabla
z_{\tau}|^{2}\varphi_{\varepsilon}^{2}\phi_{R}^{2}\leq
2\int_{\Sigma_{\lambda}}z_{\tau}|\nabla
z_{\tau}|\varphi_{\varepsilon}^{2}\phi_{R}|\nabla\phi_{R}|+2\int_{\Sigma_{\lambda}}z_{\tau}|\nabla
z_{\tau}|\varphi_{\varepsilon}\phi_{R}^{2}|\nabla\varphi_{\varepsilon}|\,.$

A simple application of the Young inequality

$8ab\leq a^{2}+16b^{2}$

allows us to deduce that

$\begin{split}\int_{\Sigma_{\lambda}}|\nabla
z_{\tau}|^{2}\varphi_{\varepsilon}^{2}\phi_{R}^{2}\leq&8\int_{\Sigma_{\lambda}\cap(B_{2R}(0)\setminus
B_{R}(0))}z_{\tau}^{2}\varphi_{\varepsilon}^{2}|\nabla\phi_{R}|^{2}+\\\
&8\int_{\Sigma_{\lambda}\cap B_{2\varepsilon}(0_{\lambda})\setminus
B_{\varepsilon}(0_{\lambda})}z_{\tau}^{2}\phi_{R}^{2}|\nabla\varphi_{\varepsilon}|^{2}\,\\\
&=I_{1}(R)+I_{2}(\varepsilon)\,.\end{split}$

Note now that

$(z_{\tau})^{2}\leq\hat{u}^{2}\leq\frac{C}{|x|^{2(N-2)}}\qquad\text{and}\qquad|\nabla\phi_{R}|\leq\frac{C}{R}$

so that, observing that $\mathcal{L}\,(B_{2R}(0)\setminus B_{R}(0))=O(R^{N})$,
it follows

$\underset{R\rightarrow+\infty}{\lim}I_{1}(R)=0$

if $N>2$. If else $N=2$ we just note that $z_{\tau}$ is bounded far from the
origin. Furthermore, as recalled here above, by a capacitary argument,
$\phi_{R}$ can be chosen in such a way that
$\displaystyle\int_{B_{R^{2}}(0)\setminus B_{R}(0)}|\nabla\phi_{R}|^{2}=o(1)$
so that again we recover that $I_{1}(R)$ goes to zero as $R$ tends to
$\infty$. On the other hand

$\int_{\Sigma_{\lambda}\cap B_{\varepsilon}(0_{\lambda})\setminus
B_{\varepsilon^{2}}(0_{\lambda})}z_{\tau}^{2}\phi_{R}^{2}|\nabla\varphi_{\varepsilon}|^{2}\leq
C\int_{\Sigma_{\lambda}\cap B_{\varepsilon}(0_{\lambda})\setminus
B_{\varepsilon^{2}}(0_{\lambda})}|\nabla\varphi_{\varepsilon}|^{2}=o(1)$

since $z_{\tau}$ is bounded in $B_{\varepsilon}(0_{\lambda})\setminus
B_{\varepsilon^{2}}(0_{\lambda})$. Note that here we exploit the fact that
$z_{\tau}\leq\hat{u}$ by construction. Consequently we also proved that

$\underset{\varepsilon\rightarrow 0}{\lim}\,\,I_{2}(\varepsilon)=0\,.$

Concluding, passing to the limit, we get that

$\int_{\Sigma_{\lambda}}|\nabla z_{\tau}|^{2}=0$

causing that, by the arbitrariness of $\tau>0$, it follows that

$\hat{u}\leq\hat{u}_{\lambda}\qquad\text{in}\quad\Sigma_{\lambda}\quad\text{for
any}\,\,\lambda<0\,.$

Repeating the argument in the opposite $(-x_{1})$-direction we conclude that,
necessarily,

$\hat{u}(x_{1},x^{\prime})=\hat{u}(-x_{1},x^{\prime})\,.$

By the definition of the Kelvin transformation, this information translates
into the same symmetry property for $u$, namely:

$u(x_{1},x^{\prime})=u(-x_{1},x^{\prime})\,.$

Taking into account the translation invariance of the equation $-\Delta
u=u^{-\gamma}$ with respect to the $x_{1}$-direction, it is easy to see that
$u$ has to be symmetric with respect to any hyperplane $\\{x_{1}=\lambda\\}$
which causes that necessarily $u$ is constant in the $x_{1}$-direction.
Exploiting now the rotation invariance of the problem with respect to any
direction parallel to the boundary, we finally deduce that:

$u=u(x_{N})\,.$

∎

###### Proof of Theorem 2.

Once we know that the solutions have 1-D symmetry the proof of the result, in
the case $\gamma>1$, follows by the ODE analysis carried out in [13].

Let us skip therefore to the case $0<\gamma\leq 1$ to prove that, in this
case, there are no global solutions. We analyse the equation

$-u^{\prime\prime}(t)=\frac{1}{u(t)^{\gamma}}\qquad\text{for}\quad t>0$

with $u(0)=0$. We write the identity
$-u^{\prime\prime}(t)u^{\prime}(t)=\frac{u^{\prime}(t)}{u(t)^{\gamma}}$ and
integrate from $1$ to $t>1$ to deduce that

$\frac{1}{2}(u^{\prime}(t))^{2}=\frac{1}{2}(u^{\prime}(1))^{2}+F(u(1))-F(u(t))$
(11)

where $F(\cdot)$ is a diverging function since $\gamma\leq 1$. By convexity we
know that $u^{\prime}(t)$ is a decreasing function. Therefore necessarily
$u^{\prime}(t)>0$ for any $t>0$ because otherwise the solution cannot be
global and the claim is proved. Consequently $u(t)$ is monotone increasing. If
now

$\underset{t\rightarrow+\infty}{\lim}\,u(t)=+\infty$

we easily get a contradiction by (11) and the fact that
$\underset{u\rightarrow+\infty}{\lim}\,F(u)=+\infty$. If else

$\underset{t\rightarrow+\infty}{\lim}\,u=L>0,$

we observe that there would exists a sequence of points $t_{n}$ such that
$\underset{n\rightarrow+\infty}{\lim}\,u^{\prime\prime}(t_{n})=0$, which is
not possible since the right hand side of the equation do not approaches zero
at such a sequence of points. This proves the claim and the theorem. ∎

Data availability statement

All data generated or analyzed during this study are included in this
published article.

Conflict of interest statement

The authors declare that they have no competing interest.

## References

  * [1] D. Arcoya, L. Boccardo, T. Leonori and A. Porretta,  Some elliptic problems with singular natural growth lower order terms.  J. Differential. Equations 249 (11) (2010), 2771–2795.
  * [2] L. Boccardo and F. Murat,  Almost everywhere convergence of the gradients of solutions to elliptic and parabolic equations. _Nonlinear Anal._ 19 (1992), 405–412.
  * [3] L. Boccardo and L. Orsina,  Semilinear elliptic equations with singular nonlinearities.  Calc. Var. Partial Dif. 37 (3-4) (2010), 363–380.
  * [4] B. Brandolini, F. Chiacchio and C. Trombetti,  Symmetrization for singular semilinear elliptic equations.  Ann. Mat. Pura Appl. 4 (2) (2014), 389–404.
  * [5] A. Canino and M. Degiovanni,  A variational approach to a class of singular semilinear elliptic equations.  J. Convex Anal. 11 (1) (2004), 147–162.
  * [6] A. Canino, F. Esposito and B. Sciunzi, On the Höpf boundary lemma for singular semilinear elliptic equations, J. Differential Equations. 266 9 (2019), 5488-5499.
  * [7] M.G. Crandall , P.H. Rabinowitz and L. Tartar ,  On a Dirichlet problem with a singular nonlinearity, Comm. in PDE, 2 2, (1977), 193–222.
  * [8] M. del Pino,  A global estimate for the gradient in a singular elliptic boundary value problem.  _Proc. Roy. Soc. Edinburgh_ Sect. A 122 (1992), no. 3-4, 341?352.
  * [9] F. Esposito and B. Sciunzi , On the Höpf boundary lemma for quasilinear problems involving singular nonlinearities and applications, J. Functional Analysis, 278 4, (2020), 108346.
  * [10] A. Ferone, A. Mercaldo and S. Segura De Leon,  A singular elliptic equation and a related functional.  _ESAIM: COCV_ 27 (2021), 39.
  * [11] D. Gilbarg and N.S. Trudinger, Elliptic Partial Differential Equations of Second Order, Springer, Berlin, 2nd Edition, 1983.
  * [12] A.C. Lazer and P.J. McKenna,  On a singular nonlinear elliptic boundary-value problem.  Proc. Amer. Math. Soc. 111 (1991), 721–730.
  * [13] L. Montoro, L. Muglia and B. Sciunzi,  Classification of solutions to $-\Delta u=u^{-\gamma}$ in the half-space, Math. Ann. , https://doi.org/10.1007/s00208-023-02717-4
  * [14] F. Oliva and F. Petitta,  Finite and infinite energy solutions of singular elliptic problems: existence and uniqueness,  J. Differential Equations 264, n1 (2018), 311–340.
  * [15] J. Serrin, Local behavior of solutions of quasi-linear equations, Acta Math., 111, (1964), 247–302.


# Structure and Internal Kinematics of Nine Inner Milky Way Globular
Clusters111Based on observations made with the NASA/ESA Hubble Space
Telescope, obtained at the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
under NASA contract NAS 5-26555. These observations are associated with
program GO-15065.

Roger E. Cohen  Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218, USA  Andrea Bellini  Space Telescope Science Institute,
3700 San Martin Drive, Baltimore, MD 21218, USA  Mattia Libralato  Space
Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
Matteo Correnti  Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218, USA  Thomas M. Brown  Space Telescope Science Institute,
3700 San Martin Drive, Baltimore, MD 21218, USA  Jason S. Kalirai  Johns
Hopkins University Applied Physics Laboratory, 11101 Johns Hopkins Road,
Laurel, MD 21723, USA

###### Abstract

This study constitutes part of a larger effort aimed at better characterizing
the Galactic globular clusters (GGCs) located towards the inner Milky Way
bulge and disk. Here, we focus on internal kinematics of nine GGCs, obtained
from space-based imaging over time baselines of $>$9 years. We exploit
multiple avenues to assess the dynamical state of the target GGCs,
constructing radial profiles of projected stellar density, proper motion
dispersion, and anisotropy. We posit that two-thirds (6/9) of our target GGCs
are in an advanced dynamical state, and are close to (or have recently
undergone) core collapse, supported by at least two lines of evidence: First,
we find relatively steep proper motion dispersion profiles, in accord with
literature values for core-collapsed GGCs. Second, we find that our sample is,
in the mean, isotropic even out to their half-light radii, although one of our
target clusters (NGC 6380) is tangentially anisotropic at $>$1$\sigma$ beyond
its half-light radius, in accord with theoretical predictions for clusters
evolving in strong tidal fields. Our proper motion dispersion and anisotropy
profiles are made publicly available.

††facilities: HST (ACS/WFC,WFC3/UVIS)

##  1 Introduction

The internal kinematics of Galactic globular clusters (GGCs) provides a
valuable window into their dynamical state. At the same time, clusters serve
as critical testbeds for comparison with simulations of cluster dynamical
evolution. In particular, by combining cluster structural parameters with
radial profiles of proper motion dispersion and anisotropy, Watkins et al.
(2015a) were able to characterize a relationship between anisotropy and
relaxation time, showing that dynamically older clusters were isotropic out to
larger radii. They also found that more concentrated clusters exhibited
steeper slopes in their dispersion profiles. This study has the advantage that
it addresses a statistically significant sample of 22 GGCs with a range of
dynamical ages, providing a solid empirical footing to which other samples may
be compared.

Observational characterization of the internal kinematics of GGCs also
provides important constraints on theoretical predictions. Specifically,
simulated GGCs never reach full energy equipartion (Baumgardt & Makino, 2003;
Trenti & van der Marel, 2013; Bianchini et al., 2016; Torniamenti et al.,
2019), largely consistent with measurements of internal kinematics over a
useful baseline in stellar mass within GGCs (Bellini et al., 2018; Libralato
et al., 2018, 2019; Watkins et al., 2020). In addition, the observed behavior
of GGC radial anisotropy profiles (Watkins et al., 2015a; Bellini et al.,
2017a; Libralato et al., 2019; Raso et al., 2020) is consistent with model
predictions that tidally underfilling clusters generally relax from the inside
out (e.g. Tiongco et al., 2016a; Zocchi et al., 2016).

Such data-model comparisons are feasible largely because the observations
consist of proper motions rather than line-of-sight velocity measurements.
While the availability of both is ideal (i.e. Watkins et al., 2015b; Bellini
et al., 2017a), proper motions have several advantages: First, they describe
stellar motions in two orthogonal components, allowing for a measurement of
anisotropy. Second, they generally yield larger sample sizes, in many cases
over a significant baseline of stellar mass. Third, because proper motions
result from the comparison of entire images rather than targeting individual
sources a priori, information concerning the relative motions of field stars
and background galaxies is obtained as well.

Despite their advantages, the use of relative proper motions for GGC internal
kinematics requires exquisite photometric and astrometric precision, ideally
in more than one observational epoch, which became feasible only recently. For
this reason, the GGCs located towards the inner Milky Way bulge and disk have,
until now, been excluded from any analysis of internal kinematics; the
necessary deep, multi-epoch imaging simply did not exist. Furthermore, in the
rare cases where even single-epoch deep imaging has been analyzed, GGCs
towards the bulge turn out to have structural parameters differing
significantly from their catalog values based on optical integrated light
(e.g. Lanzoni et al., 2010; Cohen et al., 2014; Saracino et al., 2015).
Meanwhile, existing compilations, taken at face value, suggest that inner
Milky Way GGCs are, from a dynamical point of view, a particularly interesting
subset: They are preferentially concentrated, with a high incidence of core-
collapse (Harris, 1996, 2010 edition, hereafter H96; Baumgardt et al. 2019).
Therefore, inner Milky Way GGCs represent a valuable opportunity to build up a
sample of dynamically evolved clusters in a strong tidal field. From the
standpoint of internal kinematics and structure, such a sample is ideal not
only for empirical comparison with a large ensemble of well-studied GGCs
(Watkins et al., 2015a), but also with theoretical predictions regarding
clusters evolving in tidal fields (Baumgardt & Makino, 2003; Vesperini et al.,
2014; Sollima et al., 2015; Bianchini et al., 2017a).

With this in mind, we exploit our recent second-epoch Hubble Space Telescope
(HST) imaging of nine GGCs located towards the inner Milky Way (GO-15065,
PI:Cohen) to study the internal kinematics of these clusters. This imaging was
obtained with the primary goal of cleaning and homogenizing cluster color-
magnitude diagrams (CMDs) of clusters with extant first-epoch archival imaging
to measure global cluster properties (R. E. Cohen et al. 2020, in prep.), but
we found that the precision of the available proper motions is sufficient for
an internal kinematic analysis, at least of a fixed-mass population (see Sect.
4.1).

The remainder of this paper is organized as follows: In the next section, we
summarize our observations and data reduction strategy to yield precision
astrophotometric catalogs. In Sect. 3 we use single-mass radial number density
profiles to measure the structural parameters of our target clusters, and in
Sect. 4 we present radial profiles of proper motion dispersion and anisotropy.
In Sect. 5 we use our results to assess the dynamical state of our target
clusters, and in the final section we summarize our results and discuss
avenues for future investigation.

##  2 Observations and Data Reduction

We make use of multi-epoch HST imaging of nine GGCs located towards the Milky
Way bulge and inner disk. These clusters were selected based on extant first-
epoch HST imaging of sufficient quality to enable precision astrometry
faintward of their main sequence turnoffs (MSTOs). We have obtained second-
epoch photometry of all nine of these clusters with ACS/WFC onboard HST in the
F606W and F814W filters (GO-15065, PI:Cohen), facilitating multi-epoch proper
motion analyses while homogenizing the photometry to the same system as the
large, well-studied sample of 65 GGCs observed in the HST GGC Treasury Survey
(Sarajedini et al., 2007). The observations we use are summarized in Table 1
along with catalog values of the cluster heliocentric distances and iron
abundances, although these should probably only be regarded as indicative (we
return to this point later in Sect. 6) and are included to demonstrate the
approximate range of values spanned by the target clusters. The multi-epoch
imaging we use covers time baselines from 9.1 to 16.0 yr, sufficient not only
for cluster-field separation based on relative proper motions, but also for
analysis of the internal kinematics of our target clusters, which is the focus
of the present study.

The photometric and astrometric reduction techniques we use have been
presented in detail elsewhere (Bellini et al., 2014, 2017b, 2018; Libralato et
al., 2018), and a detailed description of the observations, proper-motion-
cleaned CMDs and differential reddening corrections for our target clusters is
presented in a companion paper (Cohen et al. in prep). Therefore, we summarize
them very briefly here for convenience.

For each target cluster, point spread function fitting (PSF) photometry is
performed separately for each epoch, in multiple passes, as described in
Bellini et al. (2017b) and Bellini et al. (2018), and the resulting
instrumental catalogs are calibrated to the Vegamag system (e.g. Bellini et
al., 2017b; Nardiello et al., 2018; Libralato et al., 2019). Next, per-image
distortion-corrected positions are used to calculate individual stellar proper
motions relative to the bulk cluster motion using a least-squares fit of a
straight line to the position of each star versus time in each coordinate.
This procedure is performed iteratively, and the resulting proper motions are
corrected for any remaining systematic residuals as a function of color,
magnitude, and position, propagating in quadrature the uncertainties on these
corrections (Bellini et al., 2014, 2018; Libralato et al., 2018). The
photometry is corrected for differential reddening, also iteratively, by
shifting each star along the reddening vector to place it on an empirically
determined fiducial sequence, by an amount determined using local, well-
measured proper motion members (e.g. Milone et al., 2012; Bellini et al.,
2017c).

We make several cuts in our astrophotometric catalogs to retain only well-
measured stellar sources (note that an additional set of more stringent cuts
are necessary for internal kinematics, described in Sect. 4.1). These cuts are
applied simultaneously (rather than sequentially) on a per-filter, per-epoch
basis using various photometric diagnostic parameters as follows:

  1. 1.

The fractional flux within the PSF fitting radius from neighbors (prior to
neighbor subtraction) is denoted $o$. We set a fixed upper limit on $o$ on a
per-filter, per-cluster basis by retaining only the stars with the 70% of
smallest $o$ values. In practice, this resulted in maximum values of $o<$0.7
in all cases, a slightly harsher cut than other studies which typically make
cuts at a fixed value of $o<$1 (e.g. Bellini et al., 2017c). However, this
turns out to be inconsequential to our results, and the impact of our $o$ cut
is similar across our sample of target clusters. Since $o$ describes relative,
rather than absolute flux contribution from neighbors, the stars which passed
our other quality cuts described below and failed only our $o$ cuts are
preferentially faint, with 97% of them in the median (and $<$90% in all cases)
lying faintward of the upper main sequence and red giant branch where we
perform our analyses (see Sects. 3.2 and 4.1).

  2. 2.

RADXS is the ratio between the flux outside the PSF fitting radius to the flux
predicted by the model PSF, and can therefore be positive (i.e. for galaxies)
or negative (i.e. for cosmic rays). We retain only sources with
-0.05$\leq$RADXS$\leq$0.05.

  3. 3.

The QFIT value is essentially a linear correlation coefficient between
observed pixel values and those predicted by the PSF model, so that a perfect
fit will have QFIT=1 and worse PSF fits will have 0$<$QFIT$<$1\. Because QFIT
becomes strongly magnitude-dependent towards fainter magnitudes (e.g. fig. 8
of Bellini et al. 2017b), we exclude stars with the lowest 7.5% of QFIT values
at their magnitude, and exclude all stars with QFIT$<$0.7 regardless of
magnitude. As with the cut on $o$, the impact of the QFIT cut does not
appreciably vary over our target cluster sample since we intentionally kept
the observing strategy nearly identical for all of the second-epoch imaging we
use here.

In the left panel of Fig. 1, we show an example differential-reddening-
corrected CMD from our second-epoch imaging of NGC 6355. This CMD shows all of
the sources that passed our quality cuts, with the exception of saturated
stars, which are shown as grey crosses and excluded from any further analysis
here. Of the sources passing our quality cuts, those without measured proper
motions are shown as grey points, most of which lie at the faint end of the
CMD because our second epoch imaging is, by design, deeper than the first
epoch archival imaging. The upper and lower right-hand panels of Fig. 1 refer
to additional more stringent cuts made to select a sample for kinematic
analysis (see Sect. 4.1).

We perform artificial star tests to quantify photometric incompleteness as a
function of color, magnitude, and distance from the cluster center (i.e.
crowding) in our second epoch imaging. For each cluster, 106 artificial stars
are inserted, one at a time, with a realistic luminosity function (i.e.
exponentially increasing towards fainter magnitudes) and spatial distribution
(more crowded towards the cluster center, based on an exponential radial
density profile), and given a color distribution based on the observed CMD.
Artificial star tests are performed using software based on that used for the
ACS GGC Treasury Survey (Anderson et al., 2008), and artificial stars are
considered recovered if they are found within 0.5 pix of their input positions
and pass all of the quality cuts listed in items 1-3 above.

Figure 1: Illustration of stars passing our quality cuts, for the example case of NGC 6355. (Left:) CMD, corrected for differential reddening, showing all stars passing the basic photometric quality cuts described in Sect. 2, except for saturated stars, which are shown as grey crosses. Stars without valid proper motions are shown in grey, whereas proper-motion-selected cluster and field stars are shown using black and green circles respectively. (Upper Right:) Vector point diagram showing relative proper motions of cluster and field stars. (Lower Right:) CMD and proper motion uncertainties of cluster stars passing the basic cuts listed in Sect. 2 (black), zoomed in on the upper main sequence and red giant branch, to highlight the subset of these that pass the more stringent photometric and astrometric quality cuts for internal kinematics described in Sect. 4.1 (blue points), and the mono-mass sample selected within this subset for kinematic analysis based on their proper motion uncertainties (red points).  Table 1: Summary of Observations  Cluster |  $R_{\sun}$aaFrom Valenti et al. (2007, 2010) where available, otherwise H96. | [Fe/H]bbFrom Carretta et al. (2009). | $t_{base}$ | Epoch | GO | Filter | Exposures
---|---|---|---|---|---|---|---
| kpc | dex | yr |  |  |  |
NGC 6256 | 9.1 | -0.62$\pm$0.09 | 9.871 | Aug 2009 | 11628 | WFC3/UVIS F555W | 3$\times$360s
|  |  |  |  | 11628 | WFC3/UVIS F814W | 3$\times$100s
|  |  |  | Jun 2019 | 15065 | ACS/WFC F606W | 1$\times$60s $+$ 4$\times$502s
|  |  |  |  | 15065 | ACS/WFC F814W | 1$\times$60s $+$ 4$\times$505s
NGC 6325 | 7.8 | -1.37$\pm$0.14 | 9.131 | May 2010 | 11628 | WFC3/UVIS F438W | 3$\times$435s
|  |  |  |  | 11628 | WFC3/UVIS F555W | 3$\times$85s
|  |  |  | Jun 2019 | 15065 | ACS/WFC F606W | 1$\times$30s $+$ 4$\times$494s
|  |  |  |  | 15065 | ACS/WFC F814W | 1$\times$30s $+$ 4$\times$493s
NGC 6342 | 8.4 | -0.49$\pm$0.14 | 9.871 | Aug 2009 | 11628 | WFC3/UVIS F438W | 3$\times$420s
|  |  |  |  | 11628 | WFC3/UVIS F555W | 3$\times$80s
|  |  |  | Jun 2019 | 15065 | ACS/WFC F606W | 1$\times$10s $+$ 4$\times$493s
|  |  |  |  | 15065 | ACS/WFC F814W | 1$\times$10s $+$ 4$\times$493s
NGC 6355 | 9.0 | -1.33$\pm$0.14 | 9.881 | Aug 2009 | 11628 | WFC3/UVIS F438W | 3$\times$440s
|  |  |  |  | 11628 | WFC3/UVIS F555W | 3$\times$80s
|  |  |  | Jun 2019 | 15065 | ACS/WFC F606W | 1$\times$30s $+$ 4$\times$495s
|  |  |  |  | 15065 | ACS/WFC F814W | 1$\times$30s $+$ 4$\times$496s
NGC 6380 | 9.2 | -0.40$\pm$0.09 | 9.117 | Mar 2010 | 11628 | WFC3/UVIS F555W | 3$\times$440s
|  |  |  |  | 11628 | WFC3/UVIS F814W | 3$\times$80s
|  |  |  | Apr 2019 | 15065 | ACS/WFC F606W | 1$\times$60s $+$ 4$\times$502s
|  |  |  |  | 15065 | ACS/WFC F814W | 1$\times$60s $+$ 4$\times$502s
NGC 6401 | 7.7 | -1.01$\pm$0.14 | 16.048 | Jul 2003 | 9799 | ACS/WFC F606W | 2$\times$340s
|  |  |  |  | 9799 | ACS/WFC F814W | 2$\times$240s
|  |  |  | Aug 2019 | 15065 | ACS/WFC F606W | 1$\times$60s $+$ 4$\times$487s
|  |  |  |  | 15065 | ACS/WFC F814W | 1$\times$60s $+$ 4$\times$487s
NGC 6453 | 10.7 | -1.48$\pm$0.14 | 9.123 | May 2010 | 11628 | WFC3/UVIS F438W | 3$\times$450s
|  |  |  |  | 11628 | WFC3/UVIS F555W | 3$\times$80s
|  |  |  | Jun 2019 | 15065 | ACS/WFC F606W | 1$\times$40s $+$ 4$\times$498s
|  |  |  |  | 15065 | ACS/WFC F814W | 1$\times$30s $+$ 4$\times$499s
NGC 6558 | 7.4 | -1.37$\pm$0.14 | 15.911 | Sep 2003 | 9799 | ACS/WFC F606W | 1$\times$340s
|  |  |  |  | 9799 | ACS/WFC F814W | 2$\times$10s $+$ 1$\times$340s
|  |  |  | Aug 2019 | 15065 | ACS/WFC F606W | 1$\times$10s $+$ 4$\times$498s
|  |  |  |  | 15065 | ACS/WFC F814W | 1$\times$10s $+$ 4$\times$498s
NGC 6642 | 8.6 | -1.19$\pm$0.14 | 15.144 | Mar 2004 | 9799 | ACS/WFC F606W | 1$\times$10s $+$ 1$\times$340s
|  |  |  |  | 9799 | ACS/WFC F814W | 1$\times$10s $+$ 1$\times$340s
|  |  |  | Apr 2019 | 15065 | ACS/WFC F606W | 1$\times$10s $+$ 4$\times$493s
|  |  |  |  | 15065 | ACS/WFC F814W | 1$\times$10s $+$ 3$\times$493s

##  3 Structural Parameters

To determine the centers of our clusters and construct number density
profiles, we make use of star counts rather than integrated light. This choice
avoids systematics due to shot noise bias from a small number of very luminous
stars, which has been demonstrated to yield erroneous determinations of
cluster centers (i.e. luminosity rather than gravity centers) and artifacts in
radial profiles constructed from integrated light (e.g. Calzetti et al., 1993;
Anderson & van der Marel, 2010). We use our full-frame second-epoch ACS/WFC
imaging (i.e. without requiring proper motion information) for center
determinations and construction of density profiles due to its larger field of
view, higher signal-to-noise, and homogeneity across our target clusters
compared to the heterogenous, shallower archival first-epoch imaging.

###  3.1 Cluster Centers

####  3.1.1 Determination of Center Coordinates

We determine cluster centers by fitting ellipses to isodensity contours, as
this method outperforms symmetry-based “pie-slice”methods, as discussed by
Goldsbury et al. (2010) and Anderson & van der Marel (2010). To mitigate the
effects of incompleteness and differential reddening, we construct a density
map using only sources brighter than the MSTO in each cluster, making a
magnitude cut parallel to the reddening vector222Throughout this study, we
assume $(A_{F606W},A_{F814W})/E(B-V)$ = $(2.876,1.884)$ (Casagrande &
VandenBerg, 2014). Following Goldsbury et al. (2010), we construct an
oversampled density map using spatial bins placed every 2$\arcsec$, counting
the number of sources within 10$\arcsec$ for each bin, and an example density
map is shown in Fig. 2. Eight isodensity contours, shown in red in Fig. 2, are
generated from the resultant density map, evenly spaced between the minimum
non-zero density and maximum density, and ellipses, shown in blue in Fig. 2,
are fit to each of the innermost four isodensity contours (excluding the
central one). The cluster center and its uncertainty we report, given in Table
2, are then the mean and standard deviation of the individual ellipse centers.

Figure 2: An example density map for NGC 6401 demonstrating our ellipse
fitting procedure to determine the cluster center. The density map (shown in
greyscale according to the horizontal color bar at the top of the plot) is
constructed using stars brightward of the MSTO, making a magnitude cut
parallel to the reddening vector (see text for details). Eight contours (shown
in red) are fit to the density map between its minimum (non-zero) and maximum
value. Ellipses (shown in blue) are fit to the innermost four contours
(excluding the central one), and the location and uncertainty of the cluster
center reported in Table 2 is the mean and standard deviation of these
individual ellipse centers.

####  3.1.2 Comparison With Other Center Values

Given that our clusters are relatively concentrated (H96; also see Sect. 3.2),
precise measurements of their centers are crucial to obtaining reliable
density profiles and constraining the structure of the cluster cores. We first
compare our centers with values given in the H96 catalog, which are taken from
Shawl & White (1986) and Picard & Johnston (1995) (4 clusters each) and the
center of NGC 6380 is taken from Djorgovski & Meylan (1993). This comparison
is shown in the right-hand panel of Fig. 3. While there is a statistically
null mean offset across the sample, the individual clusters are offset from
their catalog locations typically by several arcsec. This is not due to
uncertainties in our absolute astrometry, which was calibrated to Gaia DR2
(Gaia Collaboration et al., 2018) with residuals of $<$0.1$\arcsec$, nor is it
due to the absolute motions of the clusters, although the latter contribution
is not entirely insignificant, ranging from $\sim$0.1 to 0.4$\arcsec$ over the
time since the original measurements assuming the absolute proper motions from
Baumgardt et al. (2019). The offsets between our measured centers and the
catalog values closely mirror the results of Goldsbury et al. (2010), who used
a similar technique and photometric data, and found 26/65 (40%) of their
clusters to be discrepant by more than 5$\arcsec$, while we find the same for
two of our nine target clusters, also noting that they found 8/65 (12%) to
differ by more than 10$\arcsec$, while none of our target clusters show such a
large offset. Given the similarity of our results and those of Goldsbury et
al. (2010), in combination with our center uncertainties of $\leq$1$\arcsec$
(see Table 2), we surmise that the original catalog centers (obtained from
optical integrated light in ground-based images, in many cases with
photographic plates and/or seeing worse than 1$\arcsec$) are typically
incorrect by several arcsec.

To test this hypothesis, we independently redetermine our cluster centers and
uncertainties using an alternative technique described and employed, for
example, in Montegriffo et al. (1995) and Lanzoni et al. (2019). This
alternative technique functions via iterative recentering, in which the mean
offset is calculated for each coordinate using stars inside some maximum
radius $R_{max}$ from a trial center, recentering on the mean value in each of
the two coordinates until convergence, indicated by a change of
$<$0.01$\arcsec$ in the center location since the previous iteration.
Uncertainties are obtained by rerunning the centering procedure for different
values of $R_{max}$ (from 15$\arcsec$ to 30$\arcsec$ in steps of 2.5$\arcsec$)
and calculating the standard deviation of the resulting centers. A comparison
between the centers obtained via this iterative recentering method and our
ellipse fitting method is shown in the left panel of Fig. 3, revealing
excellent agreement: The differences in center positions are all
$<$0.7$\arcsec$ (with a median difference of $<$0.4$\arcsec$), and the mean
difference in center locations (given in the left panel of Fig. 3) is not
statistically significant, so we adopt the ellipse fitting centers.

For NGC 6256, our center from ellipse fitting is located 0.35$\arcsec$ from
the center determined by Cadelano et al. (2020a) from an independent analysis
of the first epoch imaging, within their quoted uncertainty of
$\sim$0.4$\arcsec$.

Figure 3: A comparison between the cluster center coordinates we determine
using ellipse fitting to isodensity contours versus those determined by
iterative recentering (left, see text) and those from the catalog of H96
(right). In each panel, the dotted grey lines correspond to a null positional
offset, and the unweighted mean offset across the sample and its standard
deviation are reported. Given the agreement of our centers determined from two
independent methods shown in the left-hand panel, we argue that the larger
discrepancies in the right-hand panel are due to typical errors of a few
arcsec in the original catalog centers, which were based on integrated light
in shallow optical images (see text for details). Note the difference in scale
between the two panels.

###  3.2 Radial Density Profiles

####  3.2.1 Construction of Density Profiles

We measure cluster radial density profiles by using completeness-corrected
star counts, which are robust to biases inherent in integrated light profiles
(e.g. Calzetti et al., 1993). As with the cluster center determination, we
restrict our stellar sample to stars brightward of the MSTO in each cluster,
and make a magnitude cut parallel to the reddening vector. The density profile
is then calculated by dividing the sample into logarithmically-spaced radial
annuli333This serves as a compromise between spatial resolution and number
statistics. We also limit uncertainties by requiring that each radial bin has
at least 3 observed stars per azimuthal sub-sector, which affects at most the
innermost few radial bins., and dividing each annulus into azimuthal sub-
sectors (e.g. Djorgovski, 1988; Ferraro et al., 1999; Lanzoni et al., 2019).
We calculate the completeness-corrected density in each azimuthal sub-sector,
and take the density and its uncertainty in each radial annulus as the mean
and standard deviation of the densities among azimuthal sub-sectors.

Due to the size of the ACS/WFC field of view, our HST imaging only extends to
$<$200$\arcsec$ from the centers of our target clusters. Therefore, we use
star counts from Gaia DR2 in the outer parts of our density profiles to
characterize the background level (e.g. de Boer et al., 2019; Raso et al.,
2020). Special care must be taken to mitigate incompleteness in the Gaia
photometry, so for the Gaia catalog we also make a Wesenheit magnitude cut
(parallel to the reddening vector) brightward of the cluster MSTO444Our Gaia
magnitude cut is set so that the limiting Gaia $G$ Wesenheit magnitude is
13$\leq$$W_{G}$$\leq$14 where $W_{G}$ =
$G-((B_{P}-R_{P})(A_{G}/(A_{BP}-A_{RP})))$ assuming
($A_{BP},A_{G},A_{RP}$)/$E(B-V)$ = (3.374,2.740,2.035) (Casagrande &
VandenBerg, 2014). While we have intentionally set this magnitude cut well
brightward of Gaia’s photometric detection limit, incompleteness even at
bright magnitides also varies with crowding, so we only construct the Gaia
density profile outside a minimum radius $r_{min}$ of 40-80$\arcsec$
(depending on the cluster) from the cluster center. By combining Wesenheit
magnitude cuts with radius cuts in the Gaia sample (and no additional cuts,
i.e. on astrometric parameters, which we do not use), we are able to generate
a photometric sample from Gaia with incompleteness that does not vary with
radius (note that any spatial inhomogeneity in the Gaia catalog that persists
despite our cuts will propagate to uncertainties in the density profile, since
the uncertainties in each radial bin are based on the variance in the
azimuthal sub-sectors).

Importantly, by setting both the Gaia and HST ACS/WFC magnitude cuts
brightward of the cluster MSTO, the Gaia stellar number counts cover the same
narrow cluster mass range as our ACS/WFC sample, so that we can exploit their
range of radial overlap to scale the shallower wide-field (i.e. Gaia) portion
to the HST portion and combine them to increase our radial coverage (e.g.
Miocchi et al., 2013, and references therein). Specifically, the Gaia and
ACS/WFC density profiles are compared in their range of radial overlap
(between $r_{min}$ and the outer limit of the ACS/WFC field of view), and a
mean scale factor (weighted by the observational uncertainties in the profile)
is computed to scale the Gaia density profile to the ACS/WFC density profile.
The fractional uncertainty of this mean scale factor is propagated to the
final observed density profile, so that any mismatch in the shape of the Gaia
and ACS/WFC density profiles in their range of radial overlap would propagate
to the uncertainties in the final profile. Fortunately, this turns out not to
be a significant contributor to the error budget, as the uncertainties in the
scale factor are generally quite small, with a median contribution of 5.4%.

Lastly, after scaling the Gaia portion of the density profile to the ACS/WFC
portion, the background, calculated as the weighted mean of the densitity in
the outermost (flat) portion of the density profile, is subtracted
(propagating in quadrature the resulting background uncertainty). As an
external check on the structural parameters we calculate using Gaia photometry
for the outer portion of the density profile, we also construct and fit (see
Sect. 3.2.2 below) density profiles using $JK_{s}$ near-IR PSF photometry from
the Vista Variables in the Via Lactea survey (Alonso-García et al., 2018),
also using radius and Wesenheit magnitude cuts, for the clusters in common,
and obtain structural parameters that agree with our Gaia-based results to
within their uncertainties. In Fig. 4 we show density profiles of each cluster
before background subtraction in cyan, with the background and its
$\pm$1$\sigma$ uncertainty indicated using solid and dashed horizontal black
lines respectively. The final background-subtracted density profile is shown
using black filled circles.

####  3.2.2 Density Profile Fitting

We fit King (1966) profiles, generated with the publicly available LIMEPY
software package (Gieles & Zocchi, 2017)555https://github.com/mgieles/limepy,
to our projected density profiles to obtain the best-fitting structural
parameters. The best-fitting model parameters and their uncertainties are
obtained using a maximum likelihood approach, and the use of flat priors
renders the posterior probability distribution function (PDF) proportional to
the likelihood, calculated as:

$\mathcal{L}\simeq\rm{exp}\left(-\frac{1}{2}\chi^{2}\right)$ (1)

Where:

$\chi^{2}=\sum_{i}^{N}\frac{(\Delta s_{i})^{2}}{\sigma_{i}^{2}}$ (2)

Here, $\Delta s_{i}$ represents the difference in observed and model-predicted
density at each of $N$ points in the observed profile, and $\sigma_{i}$
represents the uncertainty on the observed density. We maximize the log
likelihood using the emcee affine-invariant Markov chain monte carlo (MCMC)
sampler (Foreman-Mackey et al., 2013), using 25 walkers over 300 burn-in
iterations, followed by 1500 production iterations (sufficient given
autocorrelation times of $<$25 iterations). When fitting the King profiles,
there are only two free parameters, namely the shape parameter $W_{0}$ (which
has a one-to-one relationship with the concentration parameter $c$) and the
radial scale, which can be characterized by any one of either the King radius
$r_{\rm 0}$, the half-light and half-mass radii $r_{\rm hl}$ and $r_{\rm hm}$,
the observational core (i.e. half-power) radius $r_{\rm c,obs}$, or the tidal
radius $r_{\rm t}$ since all of these radii have a fixed relationship to each
other at a given $W_{0}$ (as discussed by e.g. McLaughlin & van der Marel
2005; Chatterjee et al. 2013). The best-fit values (the median over all post-
burnin MCMC iterations) and their uncertainties (16th and 84th percentiles)
for these structural parameters are listed for our target clusters in Table 2.

The best-fitting King (1966) profile for each of our target clusters is shown
as a red line in Fig. 4, and the 1$\sigma$ (2$\sigma$) uncertainty range is
illustrated using dark grey (light grey) shading. While three of our target
clusters (NGC 6355, 6380, 6401) are well-fit by the King (1966) profiles, the
model fits underpredict the central density to some extent in nearly all of
the other cases, and we return to this point in Sect. 5.1.

Table 2: Cluster Centers and Results of King Profile Fits  Cluster | RA (J2015.5) | Dec (J2015.5) | $\sigma_{pos}$ | $W_{0}$ | $r_{\rm 0}$ | $r_{\rm t}$ | $r_{\rm c,obs}$ | $r_{\rm hl}$ | $r_{\rm hm}$ | $c$
---|---|---|---|---|---|---|---|---|---|---
| deg | deg | $\arcsec$ |  | $\arcsec$ | $\arcsec$ | $\arcsec$ | $\arcsec$ | $\arcsec$ |
NGC6256 | 254.886107 | -37.120968 | 0.2 | $9.0^{+0.6}_{-0.6}$ | $10.3^{+1.0}_{-0.9}$ | $1317^{+499}_{-395}$ | $10.1^{+0.9}_{-0.9}$ | $109^{+51}_{-35}$ | $154^{+83}_{-55}$ | $2.12^{+0.15}_{-0.18}$
NGC6325 | 259.496327 | -23.767677 | 0.6 | $9.2^{+0.4}_{-0.5}$ | $7.4^{+1.2}_{-1.2}$ | $1026^{+220}_{-189}$ | $7.3^{+1.2}_{-1.2}$ | $86^{+22}_{-17}$ | $124^{+35}_{-29}$ | $2.15^{+0.11}_{-0.12}$
NGC6342 | 260.291573 | -19.587659 | 0.2 | $8.6^{+0.3}_{-0.3}$ | $6.0^{+0.6}_{-0.7}$ | $588^{+84}_{-79}$ | $5.8^{+0.6}_{-0.6}$ | $47^{+7}_{-6}$ | $64^{+12}_{-10}$ | $2.00^{+0.08}_{-0.10}$
NGC6355 | 260.993533 | -26.352827 | 1.0 | $7.3^{+0.6}_{-0.6}$ | $13.7^{+1.5}_{-1.5}$ | $566^{+220}_{-150}$ | $13.1^{+1.3}_{-1.3}$ | $47^{+16}_{-10}$ | $62^{+20}_{-11}$ | $1.62^{+0.18}_{-0.17}$
NGC6380 | 263.618611 | -39.069530 | 0.7 | $7.3^{+0.4}_{-0.4}$ | $18.8^{+1.5}_{-1.4}$ | $751^{+174}_{-138}$ | $18.0^{+1.3}_{-1.2}$ | $62^{+12}_{-9}$ | $83^{+15}_{-10}$ | $1.60^{+0.11}_{-0.11}$
NGC6401 | 264.652191 | -23.909605 | 0.5 | $7.4^{+0.7}_{-0.7}$ | $12.2^{+1.2}_{-1.1}$ | $533^{+307}_{-169}$ | $11.7^{+1.0}_{-1.0}$ | $44^{+23}_{-11}$ | $58^{+30}_{-13}$ | $1.64^{+0.23}_{-0.20}$
NGC6453 | 267.715508 | -34.598477 | 0.3 | $9.0^{+0.7}_{-0.8}$ | $7.4^{+1.3}_{-1.2}$ | $932^{+367}_{-310}$ | $7.2^{+1.2}_{-1.2}$ | $77^{+38}_{-27}$ | $109^{+62}_{-43}$ | $2.12^{+0.17}_{-0.23}$
NGC6558 | 272.573974 | -31.764508 | 0.4 | $8.3^{+0.6}_{-0.6}$ | $7.1^{+0.8}_{-0.8}$ | $597^{+250}_{-188}$ | $6.9^{+0.7}_{-0.7}$ | $48^{+22}_{-14}$ | $64^{+35}_{-20}$ | $1.93^{+0.18}_{-0.19}$
NGC6642 | 277.975957 | -23.475602 | 0.2 | $8.8^{+0.4}_{-0.4}$ | $4.5^{+0.6}_{-0.7}$ | $503^{+106}_{-91}$ | $4.4^{+0.6}_{-0.7}$ | $41^{+9}_{-7}$ | $57^{+15}_{-12}$ | $2.06^{+0.11}_{-0.10}$
Figure 4: Radial density profiles for our target clusters. Small cyan filled
circles show densities before background subtraction, and the mean background
and its uncertainty is shown using solid and dashed black horizontal lines
respectively. The best fit King (1966) profile is shown in red, and the dark
and light grey shaded regions indicate the 1$\sigma$ and 2$\sigma$
uncertainties from the fit.

##  4 Internal Kinematics

###  4.1 Quality Cuts and Sample Selection

For the analysis of photometry and stellar density profiles in Sects. 3.1-3.2,
we applied a set of basic quality cuts to the photometry (listed in Sect. 2)
to remove sources which are non-stellar, spurious, and/or poorly
characterized. An analysis of the internal kinematics of our target clusters
requires a much more stringent set of quality cuts, as blending and crowding
can systematically bias the resulting kinematic profiles (e.g. Bellini et al.,
2014; Libralato et al., 2018). Therefore we make the following additional
quality cuts on the sample used for a kinematic analysis, in addition to those
listed in Sect. 2:

  1. (i)

We only retain stars with the best (highest) 50% of QFIT values at their
magnitude.

  2. (ii)

We only retain stars that have measured proper motions based on $\geq$85% of
the exposures in which they were detected.

  3. (iii)

We only retain stars that have a reduced $\chi^{2}\leq$ 1.3 in both
coordinates from the proper motion fit.

  4. (iv)

We use a cut in relative proper motion to remove foreground and background
disk and bulge stars.

  5. (v)

We retain only stars with proper motion uncertainties in each coordinate less
than half of the local (measured using the nearest 75 stars in magnitude and
distance from the cluster center) proper motion dispersion $\sigma_{\mu}$.

  6. (vi)

We use a mild color cut around the cluster fiducial sequence to remove
outliers and as an additional safeguard against any field contaminants that
have proper motions similar to the cluster bulk motion.

The impact of these cuts on the stellar sample is illustrated for an example
case in Fig. 1. In the left panel, stars with measured proper motions are
shown as green or black filled circles for field or cluster stars
respectively, based on their relative proper motions shown in the vector point
diagram in the upper right-hand panel of Fig. 1. In the lower right, we show a
CMD zoomed in on the MSTO region, and also plot one-dimensional proper motion
error versus magnitude. Cluster members passing the basic cuts from Sect. 2
are again shown in black, as in the other panels, but now stars that also pass
the more stringent cuts listed in items (i)-(vi) are shown in blue. Lastly, we
place bright and faint magnitude limits on our high-quality sample, shown in
red in Fig. 1. Specifically, we require F814W magnitudes brighter than 0.5 mag
faintward of the MSTO, which has the advantage of providing, again, a mono-
mass sample for kinematic analysis consistently across our target clusters.
Meanwhile, we impose a bright limit by avoiding the magnitude range where
stars are detected in only short exposures, which inflates their proper motion
uncertainties (e.g. Raso et al., 2020). While Fig. 1 represents just one
example out of our nine target clusters, it is a typical case, and it bears
mention that in all cases, the archival first epoch imaging is the limiting
factor in terms of not only photometric depth but also astrometric precision:
The first epoch imaging is undithered and consists of fewer individual
exposures than our second epoch imaging (this is why we require fairly
stringent cuts on the reduced $\chi^{2}$ values from the proper motion fits).

###  4.2 Proper Motion Dispersion

We measure the proper motion dispersion profiles of our target clusters using
our high-quality mono-mass sample (shown in red in the lower right panel of
Fig. 1 for an example case), and divide this sample into radial bins that
serve as a compromise between preserving spatial resolution while remaining
adequately populated following Watkins et al. (2015a). Within each radial bin,
the proper motion dispersion $\sigma_{\mu}$ is calculated by assuming that the
proper motion follows a Gaussian distribution centered around a cluster mean
value with a total observed scatter resulting from the quadrature sum of the
observational uncertainties on individual proper motions plus the proper
motion dispersion (e.g. Walker et al., 2006; van der Marel & Anderson, 2010;
Watkins et al., 2015a; Raso et al., 2020). The only two free parameters are
the mean proper motion and the proper motion dispersion, and we solve for
these two parameters using a maximum likelihood approach with emcee, quoting
the median of the posterior distribution as the best-fit value and the 16th
and 84th percentiles as the uncertainties. This calculation includes the
correction of van de Ven et al. (2006) to account for the underestimate of the
dispersion using maximum likelihood estimators, although this correction is
essentially insignificant and affects our measured $\sigma_{\mu}$ values by
$<$0.5% in the majority of cases (and $<$3% in all cases), well within their
uncertainties.

We then fit a 4th order monotonic decreasing polynomial, constrained to be
flat at zero radius, to proper motion dispersion as a function of radius
following Watkins et al. (2015a), also using emcee. These fits are shown,
along with their uncertainties, in Fig. 5. For comparison, we also overplot
Gaia DR2-based proper motion dispersions from Baumgardt et al. (2019) in blue.
While we caution against blindly extrapolating our fits beyond their range of
validity, they generally compare well with the proper motion dispersions in
Baumgardt et al. (2019). The minority of cases where the Baumgardt et al.
(2019) values are discrepant with our fits beyond their (necessarily larger)
uncertainties (most notably NGC 6401, and to a lesser extent NGC 6256 and NGC
6380) correspond to the more crowded, extincted sightlines within our sample,
making it more difficult to disentangle cluster and field stars using Gaia
DR2, which would cause the resultant velocity dispersions to be overestimated,
consistent with Fig. 5. The results of our polynomial fits are given in Table
3, where the fit values and their uncertainties from our maximum likelihood
fits are given for the cluster center, as well as at the King radius $r_{\rm
0}$666We use $r_{\rm 0}$ here rather than the half-power radius $r_{\rm
c,obs}$ to allow a direct comparison with Watkins et al. (2015a)., 0.5 times
the half-light radius, and the half-light radius $r_{\rm hl}$ from our density
profile fits in Sect. 3.2. Our velocity dispersion (and anisotropy) profiles
are available electronically via Table 4.

Figure 5: Proper motion dispersion profiles versus radius for our target
clusters. Fourth order monotonic, decreasing polynomial fits, constrained to
be flat at the cluster centers, are shown in red, and their 1-$\sigma$
(2-$\sigma$) uncertainties are shown using dark (light) grey shading.
Additional data from the Gaia DR2-based proper motion dispersion profiles of
Baumgardt et al. (2019) are shown in blue, and the orange vertical dashed
lines indicate the value of $r_{\rm hl}$ for each cluster from our density
profile fits in Sect. 3.2.  Table 3: Proper Motion Dispersion at Selected
Clustercentric Radii

Cluster | $\sigma_{\mu}(r=0)$ | $\sigma_{\mu}(r_{\rm 0})$ | $\sigma_{\mu}(0.5r_{\rm hl})$ | $\sigma_{\mu}(r_{\rm hl})$
---|---|---|---|---
| mas yr-1 | mas yr-1 | mas yr-1 | mas yr-1
NGC6256 | $0.171^{+0.006}_{-0.006}$ | $0.168^{+0.004}_{-0.004}$ | $0.134^{+0.003}_{-0.003}$ |
NGC6325 | $0.143^{+0.007}_{-0.007}$ | $0.142^{+0.005}_{-0.005}$ | $0.116^{+0.005}_{-0.006}$ |
NGC6342 | $0.116^{+0.004}_{-0.004}$ | $0.115^{+0.003}_{-0.003}$ | $0.102^{+0.002}_{-0.002}$ | $0.089^{+0.003}_{-0.003}$
NGC6355 | $0.134^{+0.005}_{-0.005}$ | $0.130^{+0.003}_{-0.003}$ | $0.125^{+0.002}_{-0.002}$ | $0.114^{+0.003}_{-0.003}$
NGC6380 | $0.188^{+0.003}_{-0.004}$ | $0.181^{+0.002}_{-0.002}$ | $0.173^{+0.002}_{-0.002}$ | $0.158^{+0.002}_{-0.002}$
NGC6401 | $0.172^{+0.004}_{-0.004}$ | $0.169^{+0.003}_{-0.003}$ | $0.164^{+0.002}_{-0.002}$ | $0.152^{+0.002}_{-0.002}$
NGC6453 | $0.152^{+0.004}_{-0.004}$ | $0.150^{+0.004}_{-0.004}$ | $0.123^{+0.002}_{-0.002}$ | $0.098^{+0.004}_{-0.004}$
NGC6558 | $0.127^{+0.005}_{-0.005}$ | $0.122^{+0.004}_{-0.004}$ | $0.094^{+0.002}_{-0.003}$ | $0.081^{+0.003}_{-0.003}$
NGC6642 | $0.127^{+0.005}_{-0.004}$ | $0.126^{+0.004}_{-0.004}$ | $0.109^{+0.002}_{-0.002}$ | $0.093^{+0.003}_{-0.003}$

Note. — The values given here are determined directly from maximum likelihood
fits of a 4th order decreasing polynomial (constrained to be flat in the
center) to the observed proper motion dispersion profiles as described in
Sect. 4.2.

Table 4: Radial Proper Motion Dispersion and Anisotropy Profiles

Cluster | Bin | $r_{min}$ | $r$ | $r_{max}$ | N | $\sigma_{\mu}$ | $\sigma_{\rm tan}/\sigma_{\rm rad}$
---|---|---|---|---|---|---|---
|  | $\arcsec$ | $\arcsec$ | $\arcsec$ |  | mas yr-1 |
NGC6256 | 1 | 0.646 | 4.055 | 5.960 | 49 | 0.170$\pm$0.012 | 0.856$\pm$0.126
NGC6256 | 2 | 6.064 | 8.273 | 9.988 | 97 | 0.184$\pm$0.009 | 0.951$\pm$0.104
NGC6256 | 3 | 10.000 | 12.857 | 15.216 | 128 | 0.168$\pm$0.008 | 1.033$\pm$0.095
NGC6256 | 4 | 15.238 | 17.645 | 19.948 | 128 | 0.158$\pm$0.007 | 0.994$\pm$0.085
NGC6256 | 5 | 20.144 | 22.831 | 25.563 | 153 | 0.152$\pm$0.006 | 0.947$\pm$0.083

Note. — Table 4 is published in its entirety in machine-readable format. A
portion is shown here for guidance regarding its form and content.

###  4.3 Anisotropy

We use our relative proper motions to characterize trends of anisotropy versus
radius within our target clusters, where anisotropy is defined as the ratio of
the tangential to radial proper motion dispersion $\sigma_{\rm
tan}/\sigma_{\rm rad}$. With this definition, isotropy is indicated by
$\sigma_{\rm tan}/\sigma_{\rm rad}$=1, and values of $\sigma_{\rm
tan}/\sigma_{\rm rad}$ larger (smaller) than unity indicate tangential
(radial) anisotropy. To measure anisotropy as a function of radius, we perform
a linear fit following Watkins et al. (2015a), again using emcee. The best-fit
straight lines are shown in red in Fig. 6, and in Table 5 we give the fit
anisotropy at several characteristic radii for each cluster as in Table 3, and
here we also include the linear gradient of $\sigma_{\rm tan}/\sigma_{\rm
rad}$ with radius in the second column. We find that our target clusters in
general do not show statistically significant evidence for anisotropy,
although we are only able to sample between 0.8 and 2.7$r_{\rm hl}$ depending
on the cluster in question. The only possible exceptions to this trend are NGC
6342, which appears radially anisotropic beyond $r_{\rm hl}$, and NGC 6380,
which is the only cluster showing statistically significant tangential
anistropy (also beyond $r_{\rm hl}$). In this context, we point out that while
Watkins et al. (2015a) do not find any such clear incidence of tangential
anisotropy in their sample, it may actually be expected on theoretical
grounds. In fact, N-body models predict that clusters orbiting in strong tidal
fields should become tangentially anisotropic, at least in their outer parts
(Baumgardt & Makino, 2003; Vesperini et al., 2014; Sollima et al., 2015;
Tiongco et al., 2016b; Bianchini et al., 2017a), because tidal forces will,
over time, preferentially remove stars on radial orbits.

Figure 6: Anisotropy versus radius for our target clusters, defined as $\sigma_{\rm tan}$/$\sigma_{\rm rad}$ so that clusters with $\sigma_{\rm tan}$/$\sigma_{\rm rad}$$<$1 are radially anisotropic and $\sigma_{\rm tan}$/$\sigma_{\rm rad}$$>$1 are tangentially anisotropic, with isotropy indicated by the dashed horizontal line. Linear fits are shown in red, and their 1-$\sigma$ (2-$\sigma$) uncertainty ranges are shown using dark (light) grey shading. The location of $r_{\rm hl}$ is indicated using a dashed vertical orange line as in Fig. 5.  Table 5: Proper Motion Anisotropy at Selected Clustercentric Radii  Cluster | $\delta(\sigma_{\rm tan}/\sigma_{\rm rad})/\delta r$ | $\sigma_{\rm tan}/\sigma_{\rm rad}(r=0)$ | $\sigma_{\rm tan}/\sigma_{\rm rad}(r_{\rm 0})$ | $\sigma_{\rm tan}/\sigma_{\rm rad}(0.5r_{\rm hl})$ | $\sigma_{\rm tan}/\sigma_{\rm rad}(r_{\rm hl})$
---|---|---|---|---|---
| arcsec-1 |  |  |  |
NGC6256 | $-0.001^{+0.001}_{-0.001}$ | $0.99^{+0.05}_{-0.05}$ | $0.98^{+0.04}_{-0.04}$ | $0.95^{+0.03}_{-0.03}$ |
NGC6325 | $0.000^{+0.002}_{-0.002}$ | $0.98^{+0.07}_{-0.07}$ | $0.98^{+0.06}_{-0.06}$ | $0.98^{+0.05}_{-0.04}$ |
NGC6342 | $-0.001^{+0.001}_{-0.001}$ | $0.98^{+0.05}_{-0.05}$ | $0.97^{+0.04}_{-0.04}$ | $0.95^{+0.03}_{-0.03}$ | $0.92^{+0.04}_{-0.04}$
NGC6355 | $0.000^{+0.001}_{-0.002}$ | $0.95^{+0.06}_{-0.06}$ | $0.96^{+0.04}_{-0.04}$ | $0.96^{+0.03}_{-0.03}$ | $0.97^{+0.03}_{-0.03}$
NGC6380 | $0.002^{+0.001}_{-0.001}$ | $0.94^{+0.03}_{-0.03}$ | $0.97^{+0.02}_{-0.02}$ | $0.99^{+0.01}_{-0.02}$ | $1.04^{+0.02}_{-0.02}$
NGC6401 | $0.001^{+0.001}_{-0.001}$ | $1.00^{+0.03}_{-0.04}$ | $1.00^{+0.03}_{-0.03}$ | $1.01^{+0.02}_{-0.02}$ | $1.02^{+0.02}_{-0.02}$
NGC6453 | $0.000^{+0.001}_{-0.001}$ | $1.01^{+0.05}_{-0.05}$ | $1.01^{+0.04}_{-0.04}$ | $1.01^{+0.02}_{-0.02}$ | $1.02^{+0.06}_{-0.05}$
NGC6558 | $-0.001^{+0.002}_{-0.002}$ | $1.05^{+0.06}_{-0.06}$ | $1.04^{+0.05}_{-0.05}$ | $1.03^{+0.03}_{-0.03}$ | $1.01^{+0.05}_{-0.05}$
NGC6642 | $-0.002^{+0.002}_{-0.002}$ | $1.04^{+0.07}_{-0.07}$ | $1.03^{+0.06}_{-0.06}$ | $1.00^{+0.04}_{-0.04}$ | $0.96^{+0.05}_{-0.05}$

###  4.4 Rotation

Any rotation that is present in a given target cluster will be absorbed by the
six-parameter linear transformations used to place stars in individual
exposures onto a distortion-corrected master positional reference frame. Since
we measure proper motions relative to the bulk motion of the cluster, any non-
cluster objects will then have a component imparted to their proper motion
that is identical but opposite in sign as the true rotation of the cluster.
When a suitable background population (i.e. with a proper motion that can be
characterized in a straightforward way) is present, this aspect of the proper
motion analysis can be exploited to measure cluster rotation (e.g. Massari et
al., 2013; Bellini et al., 2017a). While this is not the case for our target
clusters, Heyl et al. (2017) demonstrated that cluster rotation in the plane
of the sky can be detected via skewness in the tangential proper motion
distribution. The skewness can be characterized by the $G_{\rm 1}$ value and
its significance $Z_{\rm G1}$, where symmetric distributions have absolute
values of $G_{\rm 1}$ close to zero, and the significance level $Z_{\rm G1}$
maps directly to $p$-values such that absolute values of $Z_{\rm G1}>$2 imply
a detection at $>$2$\sigma$ (e.g. Libralato et al., 2018; Bellini et al.,
2018). We list the values of $G_{\rm 1}$, $Z_{\rm G1}$, and (for convenience)
the corresponding $p$-value in Table 6. We do not find evidence of rotation in
any of the clusters in our sample, but only in two cases (NGC 6342 and NGC
6558) can we rule out rotation at a statistically significant
($\sim$2$\sigma$) level. Such a lack of detectable rotation is perhaps
expected on both observational and theoretical grounds: Simulations show that
the radius from the cluster center where the rotation signal is strongest is
at 1-2$r_{\rm hm}$ and decreases somewhat over time (Tiongco et al., 2017),
while our high-quality proper motion sample only extends to
$\lesssim$1.5$r_{\rm hm}$ in all but one case (and only extends to $<$$r_{\rm
hm}$ for 3/9 target clusters). On the observational side, the relative
strength of rotation appears to anticorrelate with half-mass relaxation time
(Sollima et al., 2020), rendering any rotational signal very difficult to
detect given the likely old dynamical ages of our target clusters (Log $t_{\rm
rh}\lesssim$9; see Sect. 5.3).

Table 6: Skewness from Tangential Proper Motion Distribution  Cluster | $G_{\rm 1}$ | $Z_{\rm G1}$ |  $p$-value
---|---|---|---
NGC6256 | -0.060 | -1.474 | 0.141
NGC6325 | -0.020 | -0.275 | 0.783
NGC6342 | -0.131 | -2.451 | 0.014
NGC6355 | -0.003 | -0.066 | 0.947
NGC6380 | 0.025 | 0.902 | 0.367
NGC6401 | -0.018 | -0.750 | 0.453
NGC6453 | -0.045 | -1.104 | 0.269
NGC6558 | 0.095 | 1.967 | 0.049
NGC6642 | 0.016 | 0.344 | 0.731

##  5 Discussion: The Dynamical State of the Target Clusters

With radial profiles of projected density, proper motion dispersion, and
anisotropy in hand for our target clusters, we compare our sample to
properties of Milky Way globular clusters at large. Using three lines of
reasoning, detailed below, we contend that our sample, comprised
preferentially of clusters in the inner Milky Way, is highly dynamically
evolved, and that the majority of our target clusters are either undergoing
core collapse or are immediately pre- or post-core-collapse.

###  5.1 Structural Parameters and King Profile Fits

The presence of a power-law (rather than flat) core in the density profile of
a cluster indicates that it is in an advanced dynamical state.
Observationally, this is perhaps best supported by the power-law core of NGC
362 (Dalessandro et al., 2013), since its core-collapsed nature has been
independently confirmed based on internal kinematics alone (Libralato et al.,
2018). In addition, Noyola & Gebhardt (2006) showed that clusters with
significant power-law slopes in their cores also tend to have shorter
relaxation times (i.e. they are dynamically old). On the theoretical side,
simulations have consistently shown that cluster cores show a power-law slope
when they are close to core collapse (e.g. Chatterjee et al., 2013; Zocchi et
al., 2016). However, Vesperini & Trenti (2010) illustrate that even in the
presence of larger sample sizes than we have at our disposal, observational
errors render an assessment of whether a cluster is immediately pre- or post-
core-collapse ambiguous.

While any clear detection of power-law cores in our target cluster density
profiles is masked by observational uncertainties, we argue that six of our
target clusters are immediately pre- or post-core-collapse, if not undergoing
core collapse currently. There are at least three (albeit not strictly
independent) observational characteristics of core-collapsed cluster density
profiles seen in N-body models reproduced by these six clusters: First, our
best-fitting King (1966) profiles underpredict the central density (with the
exception of NGC 6642), in accord with simulations of clusters close to core
collapse (Zocchi et al., 2016). Second, all of these clusters have values of
$r_{\rm c,obs}/r_{\rm hl}\approx$ 0.1, as found from density profile fits to
simulated core-collapsed clusters (Chatterjee et al., 2013). Third, all six of
these clusters have best-fitting values of $W_{0}\gtrsim$8.5, or equivalently,
$c\gtrsim$ 1.97 to within their uncertainties. On the theoretical side, this
is consistent with fits of King (1966) profiles to simulated clusters shortly
before and after core collapse (Zocchi et al., 2016), and on the observational
side, this is in line with the value of $W_{0}$=8.83 measured from the density
profile of the core-collapsed cluster NGC 6752 (Ferraro et al., 2003). Lastly,
we note that while the H96 catalog lists two of our target clusters (NGC 6380
and NGC 6401) as candidate core-collapsed (and the rest as core-collapsed),
all of our target clusters were excluded from the surface brightness profile
fits of McLaughlin & van der Marel (2005) because the King (1966) fits by
Trager et al. (1995) yielded unreliable structural parameters.

###  5.2 Dispersion and Concentration

An empirical relationship between the steepness of the radial velocity
dispersion profile and the cluster concentration was found by Watkins et al.
(2015a). They use the ratio of the velocity dispersion at the King radius
$r_{\rm 0}$ to the velocity dispersion at the half-light radius $r_{\rm hl}$
to characterize the steepness of the velocity dispersion slope as a function
of distance from the cluster center, finding that clusters with higher
concentrations have steeper dropoffs in their velocity dispersion profiles.
Their data are reproduced in Fig. 7, illustrating that core-collapsed clusters
in their sample (shown in black) occupy a particular locus, whereas the
remainder of the clusters (shown in grey) tend to have less steep dropoffs in
their dispersion profiles. In blue, we overplot our target clusters which are
not core-collapsed (NGC 6355, 6380, 6401), while the remainder of our sample,
which we assert are close to core collapse, are shown in red. This subset
nicely coincides with the core-collapsed sample from Watkins et al. (2015a),
and shows even steeper dispersion slopes in a couple of cases (NGC 6453,
6558), albeit with larger uncertainties.

Although NGC 6256 and NGC 6325 are not shown in Fig. 7 because their half-
light radii are beyond the distance for which we were able to extract
dispersion profiles, extrapolating their dispersion profiles out to $r_{\rm
hl}$ yields values of $\sigma_{\mu}(r_{\rm 0})/\sigma_{\mu}(r_{\rm hl})\simeq$
1.6 and 1.4 respectively, which, when combined with their King (1966)
concentrations in Table 2 places them with the other core-collapsed clusters
in Fig. 7.

Figure 7: Ratio of the velocity dispersion at the King radius $r_{\rm 0}$ to
velocity dispersion at the half-light radius $r_{\rm hl}$ as a function of
concentration, after fig. 14 of Watkins et al. (2015a). Core-collapsed and
non-core-collapsed clusters from Watkins et al. (2015a) are shown in black and
grey respectively, and core-collapsed and non-core-collapsed clusters from our
target sample are shown in red and blue respectively (NGC 6256 and NGC 6325
are not shown because their velocity dispersion profiles to not extend to
$r_{\rm hl}$, but see text for further discussion).

###  5.3 Anisotropy and Relaxation Time

Another empirical relationship characterized by Watkins et al. (2015a) is the
correlation between the anisotropy (at the cluster core or half-light radii)
versus relaxation time. In particular, using relaxation times from the H96
catalog, they found that regions in the clusters with relaxation times less
than $\sim$0.55 Gyr (i.e. the cores of their sample clusters) have had time to
become isotropic, whereas most of the clusters in their sample, with longer
half-light relaxation times, have not become relaxed out to their half-light
radii. Taking to heart their suggestion that this relation may be inverted as
a means to constrain the relaxation times (i.e. dynamical state) of other
clusters, we use our linear fits from Sect. 4.3 to calculate the anisotropy
$\sigma_{\rm tan}/\sigma_{\rm rad}$ for our target clusters at their best-fit
$r_{\rm 0}$ and $r_{\rm hl}$ from Table 2. From Table 5, we find that all of
our target clusters have $\sigma_{\rm tan}/\sigma_{\rm rad}$=0.99 at their
core radii to within their uncertainties, and the weighted mean (and its
standard deviation) of $\sigma_{\rm tan}/\sigma_{\rm rad}$ at $r_{\rm 0}$ is
0.98$\pm$0.01, in good agreement with the results of Watkins et al. (2015a).
At the half-light radii, however, we find a mean anistropy value across our
sample of 1.01$\pm$0.01, or 1.00$\pm$0.01 if we extrapolate the linear fits of
$\sigma_{\rm tan}/\sigma_{\rm rad}$ out to $r_{\rm hl}$ for NGC 6256 and NGC
6325. This indicates that unlike the larger, more heterogenous sample of
Watkins et al. (2015a), our inner Milky Way sample is, on average, highly
dynamically evolved. Specifically, since our target clusters are isotropic
even out to their half-light radii, the results of Watkins et al. (2015a, see
their sect. 5.4) imply half-light relaxation times of Log $t_{\rm r,hl}$
$\lesssim$9777This range of values for $t_{\rm r,hl}$ is in good agreement
with the half-mass relaxation times given for our target clusters by Baumgardt
et al. (2019), despite their possible overestimate of proper motion
dispersions from Gaia DR2 discussed in Sect. 4.2.. Lastly, comparison with the
core-collapsed cluster NGC 362 again provides independent support for the
advanced dynamical status of our target clusters, this time via the anisotropy
profile, as Libralato et al. (2018) found that NGC 362 members from the upper
main sequence brightward are isotropic out to beyond $r_{\rm hl}$.

##  6 Conclusions and Future Prospects

We have calculated projected radial profiles of stellar density, proper motion
dispersion, and anisotropy self-consistently for nine inner Milky Way GGCs,
using MSTO-mass stars. The projected profiles reveal that our sample is, in
the mean, dynamically evolved, as revealed by their high concentrations, steep
proper motion dispersion slopes, and lack of anisotropy out to our detection
limits. In particular, three of our target clusters (NGC 6355, NGC 6380, NGC
6401) are well-fit by single-mass King (1966) models, while the remaining six
are likely currently undergoing (or very close to) core collapse, with central
densities underpredicted by the King model fits and structural parameters
consistent with other core-collapsed GGCs.

As we currently have at our disposal the most homogenous sample of dynamically
old GGCs to date, additional epochs of well-dithered high-spatial-resolution
imaging as we have obtained offer tantalizing possibilities to confirm or
refute predictions of N-body models. For example, the detection of tangential
anisotropy in NGC 6380 bolsters predictions that given sufficient tidal
forces, clusters will exhibit tangential anisotropy in their outer regions at
late times. However, these models make more specific testable predictions,
namely that tidal forces are the critical driving factor determining the
extent of this tangential anisotropy (Bianchini et al., 2017a) and that the
tangential anisotropy is largely erased after core collapse (Baumgardt &
Makino, 2003). Additional epochs of well-dithered imaging would also allow for
proper motion dispersion measurements over a significant baseline of stellar
mass to quantify the level of energy equipartion in our target clusters (e.g.
Libralato et al., 2018; Watkins et al., 2020), which our data quality
unfortunately does not currently permit. With this value in hand, the core-
collapsed nature of any cluster can be independently verified using only
internal kinematics (Libralato et al., 2018; Bianchini et al., 2018b), while
allowing a test of the predicted empirical relationship between concentration
and mass segregation (de Vita et al., 2019), especially important in light of
current discrepancies between theory and observations for two well-studied
core-collapsed GGCs (Cadelano et al., 2020b). More generally, simulations
predict that dynamically evolved ($t\gtrsim$10 $t_{\rm rc}$) clusters should
have low mass-to-light ratios of M/L$\lesssim$3 inside $r_{\rm hl}$ (and post-
core-collapse dynamical evolution will act to lower the global mass-to-light
ratio), and further predict that the exact values of M/L are metallicity
dependent (Bianchini et al. 2017b, but see Dalgleish et al. 2020).

Lastly, current and forthcoming near-IR massively multiplexed spectrographs
and integral field units will yield multi-element chemical abundances, but
also precise radial velocities for large samples of cluster members. The
resulting velocity dispersion profiles, in physical units, may be combined
with our proper motion dispersion profiles, in angular units, to estimate
dynamical distances (Watkins et al., 2015b) which are independent of
assumptions on the extinction law and its parameterization. In fact, we have
intentionally left our radial profiles in terms of angular rather than
physical units precisely because the distances to our target clusters remain,
for the most part, uncertain well beyond 10%. This is evidenced by both
comparisons between optical and near-infrared photometric distances,
(H96,Valenti et al. 2007, 2010) and substantial revisions to distances of
bulge GGCs when extinction is properly accounted for and/or additional
standard candles are available (e.g. Alonso-García et al., 2015; Cadelano et
al., 2020a).

It is a pleasure to thank the anonymous referee for a thorough reading of this
manuscript and useful comments. Support for program HST GO-15065 was provided
by NASA through a grant from the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
under NASA contract NAS 5-26555.

## References

  * Alonso-García et al. (2015) Alonso-García, J., Dékány, I., Catelan, M., et al. 2015, AJ, 149, 99
  * Alonso-García et al. (2018) Alonso-García, J., Saito, R. K., Hempel, M., et al. 2018, A&A, 619, A4. doi:10.1051/0004-6361/201833432
  * Anderson et al. (2008) Anderson, J., Sarajedini, A., Bedin, L. R., et al. 2008, AJ, 135, 2055
  * Anderson & van der Marel (2010) Anderson, J. & van der Marel, R. P. 2010, ApJ, 710, 1032
  * Baumgardt & Makino (2003) Baumgardt, H. & Makino, J. 2003, MNRAS, 340, 227
  * Baumgardt et al. (2019) Baumgardt, H., Hilker, M., Sollima, A., et al. 2019, MNRAS, 482, 5138
  * Bellini et al. (2014) Bellini, A., Anderson, J., van der Marel, R. P., et al. 2014, ApJ, 797, 115
  * Bellini et al. (2017a) Bellini, A., Bianchini, P., Varri, A. L., et al. 2017, ApJ, 844, 167
  * Bellini et al. (2017b) Bellini, A., Anderson, J., Bedin, L. R., et al. 2017, ApJ, 842, 6
  * Bellini et al. (2017c) Bellini, A., Anderson, J., van der Marel, R. P., et al. 2017, ApJ, 842, 7
  * Bellini et al. (2018) Bellini, A., Libralato, M., Bedin, L. R., et al. 2018, ApJ, 853, 86
  * Bianchini et al. (2016) Bianchini, P., van de Ven, G., Norris, M. A., et al. 2016, MNRAS, 458, 3644
  * Bianchini et al. (2017a) Bianchini, P., Sills, A., & Miholics, M. 2017, MNRAS, 471, 1181
  * Bianchini et al. (2017b) Bianchini, P., Sills, A., van de Ven, G., et al. 2017, MNRAS, 469, 4359
  * Bianchini et al. (2018b) Bianchini, P., Webb, J. J., Sills, A., et al. 2018, MNRAS, 475, L96
  * Cadelano et al. (2020a) Cadelano, M., Saracino, S., Dalessandro, E., et al. 2020, ApJ, 895, 54
  * Cadelano et al. (2020b) Cadelano, M., Dalessandro, E., Webb, J. J., et al. 2020, arXiv:2009.02333
  * Calzetti et al. (1993) Calzetti, D., de Marchi, G., Paresce, F., et al. 1993, ApJ, 402, L1
  * Carretta et al. (2009) Carretta, E., Bragaglia, A., Gratton, R., et al. 2009, A&A, 508, 695
  * Casagrande & VandenBerg (2014) Casagrande, L., & VandenBerg, D. A. 2014, MNRAS, 444, 392
  * Chatterjee et al. (2013) Chatterjee, S., Umbreit, S., Fregeau, J. M., et al. 2013, MNRAS, 429, 2881
  * Cohen et al. (2014) Cohen, R. E., Mauro, F., Geisler, D., et al. 2014, AJ, 148, 18
  * de Boer et al. (2019) de Boer, T. J. L., Gieles, M., Balbinot, E., et al. 2019, MNRAS, 485, 4906. doi:10.1093/mnras/stz651
  * de Vita et al. (2019) de Vita, R., Trenti, M., & MacLeod, M. 2019, MNRAS, 485, 5752
  * Dalessandro et al. (2013) Dalessandro, E., Ferraro, F. R., Massari, D., et al. 2013, ApJ, 778, 135
  * Dalgleish et al. (2020) Dalgleish, H., Kamann, S., Usher, C., et al. 2020, MNRAS, 492, 3859
  * Djorgovski (1988) Djorgovski, S. 1988, The Harlow-Shapley Symposium on Globular Cluster Systems in Galaxies, 126, 333
  * Djorgovski & Meylan (1993) Djorgovski, S. & Meylan, G. 1993, Structure and Dynamics of Globular Clusters, 50, 325
  * Ferraro et al. (1999) Ferraro, F. R., Paltrinieri, B., Rood, R. T., et al. 1999, ApJ, 522, 983
  * Ferraro et al. (2003) Ferraro, F. R., Possenti, A., Sabbi, E., et al. 2003, ApJ, 595, 179
  * Ferraro et al. (2009) Ferraro, F. R., Beccari, G., Dalessandro, E., et al. 2009, Nature, 462, 1028
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Conley, A., Meierjurgen Farr, W., et al. 2013, Astrophysics Source Code Library
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1
  * Gieles & Zocchi (2017) Gieles, M. & Zocchi, A. 2017, Astrophysics Source Code Library
  * Goldsbury et al. (2010) Goldsbury, R., Richer, H. B., Anderson, J., et al. 2010, AJ, 140, 1830
  * H96 (1996) Harris, W. E. 1996, AJ, 112, 1487 (H96)
  * Heyl et al. (2017) Heyl, J., Caiazzo, I., Richer, H., et al. 2017, ApJ, 850, 186. doi:10.3847/1538-4357/aa974f
  * King (1966) King, I. R. 1966, AJ, 71, 64
  * Lanzoni et al. (2010) Lanzoni, B., Ferraro, F. R., Dalessandro, E., et al. 2010, ApJ, 717, 653
  * Lanzoni et al. (2019) Lanzoni, B., Ferraro, F. R., Dalessandro, E., et al. 2019, ApJ, 887, 176
  * Libralato et al. (2018) Libralato, M., Bellini, A., van der Marel, R. P., et al. 2018, ApJ, 861, 99
  * Libralato et al. (2019) Libralato, M., Bellini, A., Piotto, G., et al. 2019, ApJ, 873, 109
  * Massari et al. (2012) Massari, D., Mucciarelli, A., Dalessandro, E., et al. 2012, ApJ, 755, L32
  * Massari et al. (2013) Massari, D., Bellini, A., Ferraro, F. R., et al. 2013, ApJ, 779, 81. doi:10.1088/0004-637X/779/1/81
  * McLaughlin & van der Marel (2005) McLaughlin, D. E. & van der Marel, R. P. 2005, ApJS, 161, 304
  * Milone et al. (2012) Milone, A. P., Piotto, G., Bedin, L. R., et al. 2012, A&A, 540, A16
  * Miocchi et al. (2013) Miocchi, P., Lanzoni, B., Ferraro, F. R., et al. 2013, ApJ, 774, 151. doi:10.1088/0004-637X/774/2/151
  * Montegriffo et al. (1995) Montegriffo, P., Ferraro, F. R., Fusi Pecci, F., et al. 1995, MNRAS, 276, 739
  * Nardiello et al. (2018) Nardiello, D., Libralato, M., Piotto, G., et al. 2018, MNRAS, 481, 3382
  * Noyola & Gebhardt (2006) Noyola, E. & Gebhardt, K. 2006, AJ, 132, 447. doi:10.1086/505390
  * Picard & Johnston (1995) Picard, A. & Johnston, H. M. 1995, A&AS, 112, 89
  * Raso et al. (2020) Raso, S., Libralato, M., Bellini, A., et al. 2020, ApJ, 895, 15
  * Saracino et al. (2015) Saracino, S., Dalessandro, E., Ferraro, F. R., et al. 2015, ApJ, 806, 152
  * Sarajedini et al. (2007) Sarajedini, A., Bedin, L. R., Chaboyer, B., et al. 2007, AJ, 133, 1658
  * Shawl & White (1986) Shawl, S. J. & White, R. E. 1986, AJ, 91, 312
  * Sollima et al. (2015) Sollima, A., Baumgardt, H., Zocchi, A., et al. 2015, MNRAS, 451, 2185
  * Sollima et al. (2020) Sollima, A., Baumgardt, H., & Hilker, M. 2020, Star Clusters: From the Milky Way to the Early Universe, 351, 516. doi:10.1017/S1743921319007099
  * Tiongco et al. (2016a) Tiongco, M. A., Vesperini, E., & Varri, A. L. 2016, MNRAS, 455, 3693
  * Tiongco et al. (2016b) Tiongco, M. A., Vesperini, E., & Varri, A. L. 2016, MNRAS, 461, 402
  * Tiongco et al. (2017) Tiongco, M. A., Vesperini, E., & Varri, A. L. 2017, MNRAS, 469, 683. doi:10.1093/mnras/stx853
  * Torniamenti et al. (2019) Torniamenti, S., Bertin, G., & Bianchini, P. 2019, A&A, 632, A67
  * Trager et al. (1995) Trager, S. C., King, I. R., & Djorgovski, S. 1995, AJ, 109, 218
  * Trenti & van der Marel (2013) Trenti, M. & van der Marel, R. 2013, MNRAS, 435, 3272
  * Valenti et al. (2007) Valenti, E., Ferraro, F. R., & Origlia, L. 2007, AJ, 133, 1287
  * Valenti et al. (2010) Valenti, E., Ferraro, F. R., & Origlia, L. 2010, MNRAS, 402, 1729
  * van der Marel & Anderson (2010) van der Marel, R. P. & Anderson, J. 2010, ApJ, 710, 1063
  * van de Ven et al. (2006) van de Ven, G., van den Bosch, R. C. E., Verolme, E. K., et al. 2006, A&A, 445, 513
  * Vesperini & Trenti (2010) Vesperini, E. & Trenti, M. 2010, ApJ, 720, L179
  * Vesperini et al. (2014) Vesperini, E., Varri, A. L., McMillan, S. L. W., et al. 2014, MNRAS, 443, L79
  * Walker et al. (2006) Walker, M. G., Mateo, M., Olszewski, E. W., et al. 2006, AJ, 131, 2114
  * Watkins et al. (2015a) Watkins, L. L., van der Marel, R. P., Bellini, A., et al. 2015, ApJ, 803, 29
  * Watkins et al. (2015b) Watkins, L. L., van der Marel, R. P., Bellini, A., et al. 2015, ApJ, 812, 149
  * Watkins et al. (2020) Watkins, L. L., van der Marel, R. P., Bellini, A., et al. 2020, IAU Symposium, 351, 544
  * Zocchi et al. (2016) Zocchi, A., Gieles, M., Hénault-Brunet, V., et al. 2016, MNRAS, 462, 696