A DYNAMICAL POTENTIAL–DENSITY PAIR FOR STAR CLUSTERS WITH NEARLY ISOTHERMAL INTERIORS

Although stellar systems can be modeled directly from observations in non-parametric ways, parametrized models are often desirable as sources of intuition, theoretical insight, or computational convenience. Clusters of stars have been modeled with a wide variety of parametrized functions; many of these are summarized in chapters 2 and 4 of Binney & Tremaine (2008). The simplest commonly used analytic models, for example Plummer (1911) and Hernquist (1990), are two parameter fits that do not allow for a variable ratio between outer (or half mass) radius and core size. Since this ratio varies, and can be quite large for realistic star clusters, three parameter models generally can offer much better fits.

In this paper we use a three-parameter potential–density pair which, in projection, fits observations of globular and open clusters quite well.¹ Although its functional form is reminiscent of Dehnen (1993), the surface density is closer to King (1962), a profile widely used for fitting observed cluster surface brightnesses. When deprojected, this model is similar to the physically motivated isothermal sphere but with the additional advantage of analytic tractability: quantities such as density, potential, and binding energy can easily be written in closed form, but cannot for the isothermal sphere or the related model of King (1966). Three important features are shared with the King profile (and the King model of King 1966): a central, constant-density core, an exterior modification that produces finite mass, and an intermediate power law zone inside the halo (or, for King, truncation) radius. In this intermediate zone our profile scales as $\rho \propto {r}^{-2}$ and the King profile as ∝r⁻³. Although globular and nuclear star clusters are our primary motivation, this potential–density pair may be useful for fitting galactic potentials, as between the core and halo radii, it possesses nearly flat rotation curves. Other possible applications include the cored dark matter profiles found in dwarf spheroidal galaxies (Burkert 2015), or the "β-models" used to fit X-ray observations of galaxy clusters (Cavaliere & Fusco-Femiano 1976).

In Section 2, we present the basic formulae for our model: exactly when possible, and in approximation when necessary. In Section 3, we consider two major applications for our model: as a theoretical testbed for studies of cluster core collapse (Section 3.1), and as a tool for fitting observations (Section 3.2). We summarize and outline future uses for this model in Section 4.

We introduce a stellar density profile, which in functional form is related to the well-known Dehnen or "η" profiles (Dehnen 1993; Tremaine et al. 1994):

$$\begin{eqnarray}&&\rho (r)={\displaystyle \frac{{\rho }_{{\rm c}}}{(1+{r}^{2}/{r}_{{\rm c}}^{2})(1+{r}^{2}/{r}_{{\rm h}}^{2})}}.\end{eqnarray} \tag{ 1 }$$

Here ${\rho }_{{\rm c}}$ is the central density, ${r}_{{\rm c}}$ is the core radius, and ${r}_{{\rm h}}$ is the outer halo radius (which approaches the half-mass radius when ${r}_{{\rm c}}\ll {r}_{{\rm h}}$). This density profile has been chosen to produce a flat central core contained interior to a roughly isothermal $\rho \propto {r}^{-2}$ cluster. The halo densities of isolated clusters are known to fall off as $\rho \propto {r}^{-3.5}$ due to outer anisotropy (Spitzer & Hart 1971; Spitzer & Shapiro 1972); we have chosen the somewhat steeper $\rho \propto {r}^{-4}$ scaling for analytic tractability.

We compare this profile to the classic isothermal sphere solution² in Figure 1. Interior to the halo radius, ${r}_{{\rm h}}$, the density profiles are quite similar, provided we set the core radius ${r}_{{\rm c}}={r}_{\mathrm{iso}}\sqrt{2}/3$, where the isothermal sphere core radius ${r}_{\mathrm{iso}}=3\sigma /\sqrt{4\pi G{\rho }_{{\rm c}}}$. This ensures an asymptotic (${r}_{{\rm c}}\ll r\ll {r}_{{\rm h}}$) match between $\rho (r)$ and the isothermal sphere.

Zoom In Zoom Out Reset image size

In projection, our model's surface brightness profile is

$$\begin{eqnarray}&&I(R)={\displaystyle \frac{{{\rm \Sigma }}_{{\rm c}}}{\Upsilon ({r}_{{\rm h}}-{r}_{{\rm c}})}}\left({\displaystyle \frac{{r}_{{\rm h}}}{\sqrt{1+{R}^{2}/{r}_{{\rm c}}^{2}}}}-{\displaystyle \frac{{r}_{{\rm c}}}{\sqrt{1+{R}^{2}/{r}_{{\rm h}}^{2}}}}\right),\end{eqnarray} \tag{ 2 }$$

where R is a projected 2D radial coordinate and Υ is the mass-to-light ratio (assumed to be constant in r). We have defined a central surface density ${{\rm \Sigma }}_{{\rm c}}=\pi {\rho }_{{\rm c}}{r}_{{\rm c}}{r}_{{\rm h}}/({r}_{{\rm h}}+{r}_{{\rm c}})$. For $R\lesssim {r}_{{\rm h}}$, a suitable choice of parameters can make Equation (2) very similar to a King profile (King 1962); we elaborate on this in Section 3.2. The potential is

$$\begin{eqnarray}{\rm \Phi }(r) & = & -{\displaystyle \frac{4\pi {\rho }_{{\rm c}}{{Gr}}_{{\rm c}}^{2}{r}_{{\rm h}}^{2}}{{r}_{{\rm h}}^{2}-{r}_{{\rm c}}^{2}}}\left[{\displaystyle \frac{{r}_{{\rm h}}}{r}}\mathrm{arctan}(r/{r}_{{\rm h}})-{\displaystyle \frac{{r}_{{\rm c}}}{r}}\mathrm{arctan}(r/{r}_{{\rm c}})\right.\\ & & +\left.{\displaystyle \frac{1}{2}}\mathrm{ln}\left({\displaystyle \frac{{r}^{2}+{r}_{{\rm h}}^{2}}{{r}^{2}+{r}_{{\rm c}}^{2}}}\right)\right].\end{eqnarray} \tag{ 3 }$$

Relatedly, the mass enclosed at a radius r is

$$\begin{eqnarray}&&M(r)={\displaystyle \frac{4\pi {r}_{{\rm c}}^{2}{r}_{{\rm h}}^{2}{\rho }_{{\rm c}}}{{r}_{{\rm h}}^{2}-{r}_{{\rm c}}^{2}}}\left[{r}_{{\rm h}}\mathrm{arctan}(r/{r}_{{\rm h}})-{r}_{{\rm c}}\mathrm{arctan}(r/{r}_{{\rm c}})\right].\end{eqnarray} \tag{ 4 }$$

This expression converges as $r\to \infty$, so the total cluster mass is

$$\begin{eqnarray}&&{M}_{\mathrm{tot}}={\displaystyle \frac{2{\pi }^{2}{r}_{{\rm c}}^{2}{r}_{{\rm h}}^{2}{\rho }_{{\rm c}}}{{r}_{{\rm h}}+{r}_{{\rm c}}}}.\end{eqnarray} \tag{ 5 }$$

When ${r}_{{\rm c}}\ll {r}_{{\rm h}}$, the half-mass radius is nearly ${r}_{{\rm h}}$ and the core mass is given by

$$\begin{eqnarray}&&{M}_{{\rm c}}\equiv M({r}_{{\rm c}})\approx {\displaystyle \frac{2(1-\pi /4)}{\pi }}{\displaystyle \frac{{r}_{{\rm c}}}{{r}_{{\rm h}}}}{M}_{\mathrm{tot}}.\end{eqnarray} \tag{ 6 }$$

The escape velocity from any point in the cluster is ${v}_{\mathrm{esc}}(r)=\sqrt{-2{\rm \Phi }(r)}$. At the cluster center, this becomes

$$\begin{eqnarray}&&{v}_{\mathrm{esc}}(0)={\displaystyle \frac{2}{\sqrt{\pi }}}\sqrt{{\displaystyle \frac{{{GM}}_{\mathrm{tot}}}{{r}_{{\rm h}}-{r}_{{\rm c}}}}}\sqrt{\mathrm{ln}({r}_{{\rm h}}/{r}_{{\rm c}})}.\end{eqnarray} \tag{ 7 }$$

The above formulae have been used in the past to study gravitational lensing (e.g., Appendix A of Elíasdóttir et al. 2007); we now derive new properties of this potential–density pair more relevant for star cluster dynamics. The total potential energy of the cluster is $W=-4\pi G{{\displaystyle \int }}_{0}^{\infty }r\rho (r)M(r){dr}$. This integral evaluates to

$$\begin{eqnarray}W & = & {\displaystyle \frac{{W}_{0}}{{r}_{{\rm h}}-{r}_{{\rm c}}}}\left[{r}_{{\rm c}}\mathrm{ln}4+{r}_{{\rm h}}\mathrm{ln}4-2{r}_{{\rm c}}\mathrm{ln}\left(1+{r}_{{\rm h}}/{r}_{{\rm c}}\right)\right.\\ & & -\left.2{r}_{{\rm h}}\mathrm{ln}\left(1+{r}_{{\rm c}}/{r}_{{\rm h}}\right)\right],\end{eqnarray} \tag{ 8 }$$

where we have defined ${W}_{0}\equiv -{{GM}}_{\mathrm{tot}}^{2}{\pi }^{-1}/({r}_{{\rm h}}-{r}_{{\rm c}})$. By the virial theorem, the total cluster binding energy $E=W/2$.

Computation of the velocity dispersion, $\sigma (r)$, is more challenging. The virially averaged value is ${\sigma }_{{\rm v}}\equiv {(-W/{M}_{\mathrm{tot}})}^{1/2}$, but finding the radial dependence is more difficult. Assuming isotropic velocities, the exact three-dimensional dispersion is

$$\begin{eqnarray}{\sigma }^{2}(r) & = & {\displaystyle \frac{3G}{\rho (r)}}{{\displaystyle \int }}_{r}^{\infty }{\displaystyle \frac{M(r\prime )\rho (r\prime )}{r{\prime }^{2}}}{dr}\prime \\ & = & {\displaystyle \frac{12\pi G{\rho }_{{\rm c}}{r}_{{\rm c}}^{2}{r}_{{\rm h}}^{2}}{{r}_{{\rm h}}^{2}-{r}_{{\rm c}}^{2}}}{\displaystyle \frac{({r}^{2}+{r}_{{\rm h}}^{2})({r}^{2}+{r}_{{\rm c}}^{2})}{{r}_{{\rm c}}^{2}-{r}_{{\rm h}}^{2}}}\times \left[Y(r,{r}_{{\rm h}},{r}_{{\rm c}})\right.\\ & & -\left.Y(r,{r}_{{\rm c}},{r}_{{\rm c}})-Y(r,{r}_{{\rm h}},{r}_{{\rm h}})+Y(r,{r}_{{\rm c}},{r}_{{\rm h}})\right],\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&Y(r,{r}_{1},{r}_{2})={{\displaystyle \int }}_{r}^{\infty }{\displaystyle \frac{{r}_{1}}{r{\prime }^{4}}}{\displaystyle \frac{\mathrm{arctan}(r\prime /{r}_{1})}{1+r{\prime }^{2}/{r}_{2}^{2}}}{dr}\prime .\end{eqnarray} \tag{ 10 }$$

Analytic approximations for $\sigma (r)$ are derived in the Appendix for the limit of ${r}_{{\rm c}}\ll {r}_{{\rm h}}$. In this regime, the central velocity dispersion is

$$\begin{eqnarray}{\sigma }_{{\rm c}}^{2}\equiv \underset{r\to 0}{\mathrm{lim}}\;{\sigma }^{2} & \approx & {\displaystyle \frac{6{{GM}}_{\mathrm{tot}}({\pi }^{2}/8-1)}{\pi {r}_{{\rm h}}}}\\ & \approx & 12\pi G{\rho }_{{\rm c}}{r}_{{\rm c}}^{2}({\pi }^{2}/8-1).\end{eqnarray} \tag{ 11 }$$

Interestingly, in this limit, ${\sigma }_{{\rm c}}^{2}/{\sigma }_{{\rm v}}^{2}\approx 1.011$. In general, clusters with ${r}_{{\rm c}}\sim {r}_{{\rm h}}$ are almost perfectly isothermal for $r\lt {r}_{{\rm h}}/2$, while those with smaller core radii see a plateau of elevated velocity dispersion for ${r}_{{\rm c}}\lt r\lt {r}_{{\rm h}}$. In this plateau the dispersion is nearly isothermal, with ${\sigma }_{{\rm p}}^{2}\approx 3{{GM}}_{\mathrm{tot}}/(\pi {r}_{{\rm h}})$ (for ${r}_{{\rm c}}\ll r\ll {r}_{{\rm h}}$). Although the temperature inversion visible for ${r}_{{\rm c}}\lesssim r\lesssim {r}_{{\rm h}}$ is not physical for long-term cluster equilibria, it is characteristic of a post-collapse cluster core undergoing gravothermal expansion. Figure 2 shows exact numerical solutions for $\sigma (r)$, our approximate analytic ${\sigma }_{\mathrm{far}}(r)$ and ${\sigma }_{\mathrm{near}}(r)$ derived in appendix, and other velocities of interest (${v}_{\mathrm{circ}}$, ${v}_{\mathrm{esc}}$).

Zoom In Zoom Out Reset image size

The isotropic distribution function (DF) can be calculated numerically using Eddington's formula, but a closed form does not appear to exist. This is in contrast to the somewhat similar King model (King 1966). It is possible, however, to modify many of the above formulae into accomodating a central black hole (BH) of mass ${m}_{\bull }$. Leaving $\rho (r)$ and I(R) unchanged, and setting ${M}_{\bull }(r)=M(r)+{m}_{\bull }$ and ${{\rm \Phi }}_{\bull }(r)={\rm \Phi }(r)-{{Gm}}_{\bull }/r$, we find

$$\begin{eqnarray}&&{W}_{\bull }=W-{\displaystyle \frac{4\pi {{Gm}}_{\bull }{r}_{{\rm c}}^{2}{r}_{{\rm h}}^{2}{\rho }_{{\rm c}}\mathrm{ln}({r}_{{\rm h}}/{r}_{{\rm c}})}{{r}_{{\rm h}}^{2}-{r}_{{\rm c}}^{2}}}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}{\sigma }_{\bull }^{2}(r) & = & {\sigma }^{2}(r)+3{{Gm}}_{\bull }\left(1+{\displaystyle \frac{{r}^{2}}{{r}_{{\rm c}}^{2}}}\right)\left(1+{\displaystyle \frac{{r}^{2}}{{r}_{{\rm h}}^{2}}}\right)\\ & & \times \left[{\displaystyle \frac{1}{r}}-{\displaystyle \frac{\pi ({r}_{{\rm h}}^{2}+{r}_{{\rm c}}{r}_{{\rm h}}+{r}_{{\rm c}}^{2})}{2{r}_{{\rm c}}{r}_{{\rm h}}({r}_{{\rm h}}+{r}_{{\rm c}})}}+{\displaystyle \frac{{r}_{{\rm h}}^{2}\mathrm{arctan}(r/{r}_{{\rm c}})}{{r}_{{\rm c}}({r}_{{\rm h}}^{2}-{r}_{{\rm c}}^{2})}}\right.\\ & & -\left.{\displaystyle \frac{{r}_{{\rm c}}^{2}\mathrm{arctan}(r/{r}_{{\rm h}})}{{r}_{{\rm h}}({r}_{{\rm h}}^{2}-{r}_{{\rm c}}^{2})}}\right].\end{eqnarray} \tag{ 13 }$$

It is not obvious, but ${\sigma }_{\bull }^{2}-{\sigma }^{2}\propto 1/r$ at large radii due to multiple high-order cancellations in Equation (13).

Under the approximation that ${r}_{{\rm c}}\ll {r}_{{\rm h}}$, we calculate the two-body relaxation time (Binney & Tremaine 2008),

$$\begin{eqnarray}&&{t}_{{\rm r}}(r)\equiv {\displaystyle \frac{6.5\times {10}^{-2}\;{\sigma }^{3}(r)}{{G}^{2}{m}_{*}\rho (r)\mathrm{ln}{\rm \Lambda }}},\end{eqnarray} \tag{ 14 }$$

with the Coulomb logarithm ${\rm \Lambda }\approx 0.4{M}_{\mathrm{tot}}/{m}_{*}$ (Spitzer & Hart 1971). Specifically, the central relaxation time

$$\begin{eqnarray}&&{t}_{{\rm r}}(0)\approx {\displaystyle \frac{0.39}{\mathrm{ln}{\rm \Lambda }}}\sqrt{{\displaystyle \frac{{r}_{{\rm c}}^{3}}{{{GM}}_{\mathrm{tot}}}}}{\displaystyle \frac{{M}_{\mathrm{tot}}}{{m}_{*}}}{\displaystyle \frac{\sqrt{{r}_{{\rm c}}{r}_{{\rm h}}}}{{r}_{{\rm c}}+{r}_{{\rm h}}}}.\end{eqnarray} \tag{ 15 }$$

Approximating the inner regions of the cluster as isothermal lets us write ${t}_{{\rm r}}(r)\approx {t}_{{\rm r}}(0)(1+{r}^{2}/{r}_{{\rm c}}^{2})(1+{r}^{2}/{r}_{{\rm h}}^{2})$, which is valid until the turnover in the velocity dispersion at $r\approx {r}_{{\rm h}}/2$.

3.1. Single Component Cluster

If we apply the potential–density pair to a cluster composed of single, identical stars of mass m_*, neglecting three-body and tidal capture processes, we can describe very simply the dynamics of core collapse. In particular, under the assumption of self-similarity, we can write

$$\begin{eqnarray}&&{\displaystyle \frac{{{dr}}_{{\rm c}}}{{dt}}}\propto -{\displaystyle \frac{{r}_{{\rm c}}}{{t}_{{\rm r}}(0)}}.\end{eqnarray} \tag{ 16 }$$

Self-similar core collapse is expected in single-mass clusters on analytic grounds (Lynden-Bell & Eggleton 1980), and this is verified in numerical simulations (Baumgardt et al. 2003, and references therein). In the context of our model, self-similarity occurs under a simplifying assumption about heat conduction. Specifically, we assume that the cluster lacks internal sources of heat (i.e., no binaries, no BHs, etc.) and set

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{cond}}\equiv {A}_{\mathrm{cond}}{\displaystyle \frac{M({r}_{{\rm c}}){\sigma }_{{\rm c}}^{2}}{2{t}_{{\rm r}}(0)}}=-{\displaystyle \frac{{{dE}}_{{\rm c}}}{{{dr}}_{{\rm c}}}}{\displaystyle \frac{{{dr}}_{{\rm c}}}{{dt}}},\end{eqnarray} \tag{ 17 }$$

where the binding energy of the core ${E}_{{\rm c}}\approx {M}_{{\rm c}}{\sigma }_{{\rm c}}^{2}/2$. This evolution is trivially self-similar for constant density cores as long as we take the core σ to be isothermal (a good approximation for our model), and set ${A}_{\mathrm{cond}}$ to a constant value. This latter step must unfortunately be done artificially; the mild "bump" in our model's σ profile prohibits a self-consistent calculation of outward heat flux.

We find an analytic solution in the limit of ${r}_{{\rm c}}\ll {r}_{{\rm h}}$ and fixed ${r}_{{\rm h}}$:

$$\begin{eqnarray}&&{r}_{{\rm c}}(t)=\sqrt{{r}_{{\rm c}}^{2}(0)-5.1{{tA}}_{\mathrm{cond}}\mathrm{ln}{\rm \Lambda }{\displaystyle \frac{{m}_{*}}{{M}_{\mathrm{tot}}}}\sqrt{{{GM}}_{\mathrm{tot}}{r}_{{\rm h}}}}\end{eqnarray} \tag{ 18 }$$

Thus, the core-collapse timescale in this limit is

$$\begin{eqnarray}{T}_{\mathrm{cc}} & = & {\displaystyle \frac{1}{2{A}_{\mathrm{cond}}}}{t}_{{\rm r},0}(0)\\ & \approx & {\displaystyle \frac{0.19}{{A}_{\mathrm{cond}}\mathrm{ln}{\rm \Lambda }}}\sqrt{{\displaystyle \frac{{r}_{{\rm h}}^{3}}{{{GM}}_{\mathrm{tot}}}}}{\displaystyle \frac{{M}_{\mathrm{tot}}}{{m}_{*}}}{\displaystyle \frac{{r}_{{\rm c}}^{2}}{{r}_{{\rm h}}^{2}}},\end{eqnarray} \tag{ 19 }$$

where ${t}_{{\rm r},0}(0)$ is the initial central relaxation time. Isotropic Fokker–Planck models take $\approx 330$ central relaxation times (Cohn 1980) to reach total core collapse, implying ${A}_{\mathrm{cond}}\approx 1.5\times {10}^{-3}$. Direct N-body simulations are in reasonable agreement with this (Baumgardt et al. 2003).

With this calibration, we can now self-consistently model the evolution of our cluster. Heat flowing outward from the collapsing core will inflate the outer regions (assuming an isolated cluster unaffected by external tides) so as to satisfy ${dW}/{dt}=0$. This produces the differential equation

$$\begin{eqnarray}&&{\displaystyle \frac{{{dr}}_{{\rm h}}}{{dt}}}={\displaystyle \frac{{{dr}}_{{\rm c}}}{{dt}}}{\displaystyle \frac{W-\mathrm{ln}\left[{\displaystyle \frac{1}{2}}(1+{r}_{{\rm h}}/{r}_{{\rm c}})\right]{W}_{0}}{W+\mathrm{ln}\left[{\displaystyle \frac{1}{2}}(1+{r}_{{\rm c}}/{r}_{{\rm h}})\right]{W}_{0}}}.\end{eqnarray} \tag{ 20 }$$

Equations (16) and (20) provide a simple description of a cluster undergoing core collapse while conserving energy. We plot the core radii and outer radii of core collapsing star clusters in Figure 3. We see that ${r}_{{\rm h}}$ evolves very little, justifying the approximation used in Equation (18).

Zoom In Zoom Out Reset image size

Competing against core collapse is the evaporation of the cluster (Ambartsumian 1938), due to losses from the high-velocity tail of the velocity distribution. We assume a Maxwellian velocity distribution $f(v)\propto \mathrm{exp}(-\frac{3{v}^{2}}{2{\sigma }^{2}})$, and use this to define the fraction of stars in the distribution that will escape the cluster,

$$\begin{eqnarray}&&{\zeta }_{\mathrm{esc}}(r)={\displaystyle \frac{{{\displaystyle \int }}_{{v}_{\mathrm{esc}}(r)}^{\infty }f(v){v}^{2}{dv}}{{{\displaystyle \int }}_{0}^{\infty }f(v){v}^{2}{dv}}}\approx \sqrt{{\displaystyle \frac{6}{\pi }}}{\displaystyle \frac{{v}_{\mathrm{esc}}(r)}{\sigma (r)}}{e}^{-\frac{3}{2}\frac{{v}_{\mathrm{esc}}^{2}(r)}{{\sigma }^{2}(r)}}.\end{eqnarray} \tag{ 21 }$$

The approximate equality assumes ${v}_{\mathrm{esc}}(r)\gtrsim \sigma (r)$. Assuming that this fraction is replenished on a local relaxation time, ${t}_{{\rm r}}(r)$, the total rate of stellar evaporation from the cluster is ${\dot{N}}_{\mathrm{esc}}={{\displaystyle \int }}_{0}^{\infty }{\displaystyle \frac{d{\dot{N}}_{\mathrm{esc}}}{{dr}}}{dr}$, where

$$\begin{eqnarray}&&{\displaystyle \frac{d{\dot{N}}_{\mathrm{esc}}}{{dr}}}{dr}=4\pi {r}^{2}{\displaystyle \frac{\rho (r)}{{m}_{*}}}{\displaystyle \frac{{\zeta }_{\mathrm{esc}}(r)}{{t}_{{\rm r}}(r)}}{dr}.\end{eqnarray} \tag{ 22 }$$

We plot $d{\dot{N}}_{\mathrm{esc}}/d\mathrm{ln}r$ in Figure 4. Most of the mass flux out of the cluster comes from radii ${r}_{{\rm c}}\lesssim r\lesssim {r}_{{\rm h}}$. For fixed ${M}_{\mathrm{tot}}$ and ${r}_{{\rm h}}$, decreasing the core radius ${r}_{{\rm c}}$ will marginally increase the total evaporation rate ${\dot{N}}_{\mathrm{esc}}$. Unsurprisingly, the specific evaporation rate is much higher outside the core than in it. If we approximate ${\dot{N}}_{\mathrm{esc}}\approx d\dot{N}/d\mathrm{ln}r{| }_{r={r}_{{\rm h}}}$, we can compute an evaporation timescale ${T}_{\mathrm{evap}}=N/{\dot{N}}_{\mathrm{esc}}$ as

$$\begin{eqnarray}&&{T}_{\mathrm{evap}}\approx {\displaystyle \frac{14}{\mathrm{ln}{\rm \Lambda }}}\sqrt{{\displaystyle \frac{{r}_{{\rm h}}^{3}}{{{GM}}_{\mathrm{tot}}}}}{\displaystyle \frac{{M}_{\mathrm{tot}}}{{m}_{*}}},\end{eqnarray} \tag{ 23 }$$

in the limit where ${r}_{{\rm c}}\ll {r}_{{\rm h}}$. Comparing to Equation (19), we see that for a cluster without internal heat sources, ${T}_{\mathrm{evap}}/{T}_{\mathrm{cc}}\approx 0.1{({r}_{{\rm h}}/{r}_{{\rm c}})}^{2}$: core collapse will proceed relatively unimpeded by evaporation if ${r}_{{\rm c}}\lesssim 0.1{r}_{{\rm h}}$, but significantly greater initial core radii will allow evaporation to dominate at early times.

Zoom In Zoom Out Reset image size

3.2. Parametrized Fitting to Observations

Despite its pleasant analytical features, this model is primarily useful insofar as it can describe realistic star clusters. In this subsection, we compare in greater detail to the well-known King surface brightness profile (King 1962),

$$\begin{eqnarray}&&{I}_{{\rm K}}(R)={\displaystyle \frac{{{\rm \Sigma }}_{{\rm c},{\rm K}}}{{\Upsilon }}}{\left({\displaystyle \frac{1}{\sqrt{1+{R}^{2}/{r}_{{\rm c},{\rm K}}^{2}}}}-{\displaystyle \frac{1}{\sqrt{1+{r}_{{\rm t},{\rm K}}^{2}/{r}_{{\rm c},{\rm K}}^{2}}}}\right)}^{2}.\end{eqnarray} \tag{ 24 }$$

The three free parameters of the King profile are a characteristic surface density ${{\rm \Sigma }}_{{\rm c},{\rm K}}$, the core radius ${r}_{{\rm c},{\rm K}}$, and the truncation radius ${r}_{{\rm t},{\rm K}}$. The numerically computed King density profile differs from Equation (1) outside a shared constant-density core, as they possess different power law indices (${r}^{-2}$ for Equation (1), ${r}^{-3}$ for King 1962). Furthermore, our $\rho (r)$ declines ∝r⁻⁴ beyond the halo radius, while the King profile falls to $\rho =0$ at $r={r}_{{\rm t},{\rm K}}$. Equation (1) is likely more appropriate for clusters which are not tidally truncated, as was seen in the application of untruncated power law models to many LMC young massive clusters (Elson et al. 1987). Applied to galactic globular clusters, these untruncated models generally outperform the King model when clusters can be resolved out to many half-mass radii, although the advantage is less clear than it is for young massive clusters (McLaughlin & van der Marel 2005), and the globulars can also be well-fit by the truncated (though more extended) model of Wilson (1975).

We can relate Equation (1) and the King profile by equating their central surface density, core radius, and enclosed mass. Respectively,

$$\begin{eqnarray}&&{\displaystyle \frac{\pi {r}_{{\rm h}}{r}_{{\rm c}}{\rho }_{{\rm c}}}{{r}_{{\rm h}}+{r}_{{\rm c}}}}={{\rm \Sigma }}_{{\rm c},{\rm K}}{\left(1-{\displaystyle \frac{1}{\sqrt{1+X}}}\right)}^{2}\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{r}_{{\rm c}}={r}_{{\rm c},{\rm K}}\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}{r}_{{\rm h}} & = & {\displaystyle \frac{1}{2\pi }}\left[{\displaystyle \frac{\sqrt{X}}{1+X}}+\mathrm{arctan}(\sqrt{X})-{\displaystyle \frac{2\mathrm{ln}(\sqrt{X}+\sqrt{1+X})}{\sqrt{1+X}}}\right]\\ & & \times {\left(1-{\displaystyle \frac{1}{\sqrt{1+X}}}\right)}^{-2}.\end{eqnarray} \tag{ 27 }$$

In these equations we have used $X\equiv {r}_{{\rm t},{\rm K}}^{2}/{r}_{{\rm c},{\rm K}}^{2}$. Together, the above equations translate the parameters, $\{{{\rm \Sigma }}_{{\rm c},{\rm K}},{r}_{{\rm c},{\rm K}},{r}_{{\rm t},{\rm K}}\}$ of a given King profile into those for our model, $\{{\rho }_{{\rm c}},{r}_{{\rm c}},{r}_{{\rm h}}\}$. The reverse procedure must be done numerically because Equation (27) is nontrivially transcendental.

In Figure 5 we show best-fit parametrizations for both I(R) and ${I}_{{\rm K}}(R)$. We use the sample of 125 Milky Way globular clusters from Trager et al. (1995), and an additional 68 LMC, SMC, and Fornax clusters from McLaughlin & van der Marel (2005, using their V-band profiles). We fit every cluster's tabulated surface brightness profile, weighting our fits appropriately (using the measurement errors at each isophote for McLaughlin & van der Marel 2005, and more ad hoc data weights provided in Trager et al. 1995). We find that within the sample, mean and median r² values for nonlinear model fitting are extremely similar for both I(R) and ${I}_{{\rm K}}(R)$. In particular, our mean r² statistic for the I(R) (${I}_{{\rm K}}(R)$) profile is 0.952 (0.942), and the corresponding median value is 0.974 (0.969). A histogram of the r² distribution is shown in Figure 5.

Zoom In Zoom Out Reset image size

We have presented a novel potential–density pair that is cored, truncated, and approximately isothermal (within the halo radius). Many important properties can be written in closed form. In this paper, we have calculated the binding energy, surface brightness profile, central velocity dispersion, central escape velocity, and central relaxation time. The primary quantity of interest that does not appear reducible to closed form is the stellar DF.

The properties of this model match, surprisingly well, those of observed star clusters. We have fit our surface brightness profile to a large number of galactic globular clusters and found generally high goodnesses of fit, comparable to or slightly better than best-fit King profiles.

Because of this model's simplicity and analytic tractability, it is well-suited for the idealized study of cluster dynamics. In this paper we have presented a closed-form description of the evaporation and core collapse of single-component, single-particle clusters, and we hope to apply this model to other dynamical processes in the future.

Finally, we note that the cluster model introduced here can facilitate inexpensive computation of the dynamics of mass segregated subclusters. In particular, this potential–density pair can provide a dynamic background potential in which to run N-body simulations of the central subcluster. We are currently using this technique to investigate intermediate mass BH formation due to collisional runaways in nuclear star clusters (A. Küpper et al. 2015, in preparation), but it could also be applied to open or globular clusters.

We thank Haldan Cohn, Bence Kocsis, Andreas Küpper, Nathan Leigh, and Scott Tremaine for helpful discussions. We also thank Jens Hjorth for pointing out the use of this potential–density pair in the strong lensing community. N.S. acknowledges support through NSF grant AST-1410950 and from the Alfred P. Sloan Foundation to Brian Metzger. This work was aided by the hospitality of the Aspen Center for Physics.

Here we derive the approximate formulas for the velocity dispersion near (${\sigma }_{\mathrm{near}}$; $r\ll {r}_{{\rm h}}$) and far (${\sigma }_{\mathrm{far}}$; $r\gg {r}_{{\rm c}}$) from the center of the cluster. The integrals $Y(r,{r}_{{\rm c}},{r}_{{\rm c}})$ and $Y(r,{r}_{{\rm h}},{r}_{{\rm h}})$ from Equation (9) evaluate to elementary functions:

$$\begin{eqnarray}Y(r,{r}_{1},{r}_{1}) & = & {\displaystyle \frac{-1}{24{r}^{3}{r}_{1}^{2}}}\left[4\pi {r}_{1}^{3}+4{{rr}}_{1}^{2}-12\pi {r}^{2}{r}_{1}+3{\pi }^{2}{r}^{3}\right.\\ & & +\left(24{r}^{2}{r}_{1}-8{r}_{1}^{3}\right)\mathrm{arctan}{\displaystyle \frac{{r}_{1}}{r}}\\ & & -\left.16{r}^{3}\mathrm{ln}\left(1+{\displaystyle \frac{{r}_{1}^{2}}{{r}^{2}}}\right)-12{r}^{3}{\mathrm{arctan}}^{2}{\displaystyle \frac{r}{{r}_{1}}}\right].\end{eqnarray} \tag{ 28 }$$

$Y(r,{r}_{{\rm h}},{r}_{{\rm c}})$ and $Y(r,{r}_{{\rm c}},{r}_{{\rm h}})$ are dilogarithms, but we can find approximate solutions by Taylor expanding the arctan terms in the integrands of $Y(r,{r}_{{\rm h}},{r}_{{\rm c}})$ and $Y(r,{r}_{{\rm c}},{r}_{{\rm h}})$ in the limit $r\prime \ll {r}_{{\rm h}}$. These integrands are dominated by contributions from $r\prime \approx {r}_{{\rm c}}$, so this expansion is valid when ${r}_{{\rm c}}\ll {r}_{{\rm h}}$. These approximate solutions are

$$\begin{eqnarray}&&Y(r,{r}_{{\rm h}},{r}_{{\rm c}})\approx {\displaystyle \frac{1}{2{r}^{2}}}+\left({\displaystyle \frac{1}{6{r}_{{\rm h}}^{2}}}+{\displaystyle \frac{1}{2{r}_{{\rm c}}^{2}}}\right)\mathrm{ln}\left({\displaystyle \frac{{r}^{2}}{{r}^{2}+{r}_{{\rm c}}^{2}}}\right)\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&Y(r,{r}_{{\rm c}},{r}_{{\rm h}})\\ &&\approx {\displaystyle \frac{1}{6{r}^{3}}}\left[\pi {r}_{{\rm c}}+r-{\displaystyle \frac{3\pi {r}^{2}{r}_{{\rm c}}}{{r}_{{\rm h}}^{2}}}+2{r}_{{\rm c}}\mathrm{arctan}{\displaystyle \frac{{r}_{{\rm c}}}{r}}\left({\displaystyle \frac{3{r}^{2}}{{r}_{{\rm h}}^{2}}}-1\right)\right]\\ &&\quad -\left({\displaystyle \frac{1}{{r}_{{\rm c}}^{2}}}+{\displaystyle \frac{3}{{r}_{{\rm h}}^{2}}}\right)\mathrm{ln}\left(1+{\displaystyle \frac{{r}_{{\rm c}}^{2}}{{r}^{2}}}\right).\end{eqnarray} \tag{ 30 }$$

Equations (9) and (28–30) define ${\sigma }_{\mathrm{near}}(r)$, the approximate $\sigma (r)$ when $r\ll {r}_{{\rm h}}$ and ${r}_{{\rm c}}\ll {r}_{{\rm h}}$. We likewise derive a ${\sigma }_{\mathrm{far}}(r)$ valid for $r\gg {r}_{{\rm c}}$ and ${r}_{{\rm c}}\ll {r}_{{\rm h}}$ by making the approximation that $r\gg {r}_{{\rm c}}$ in Equation (9). Specifically,

$$\begin{eqnarray*}\begin{array}{ccc}{\sigma }_{\mathrm{far}}^{2}(r) & = & {\displaystyle \frac{6{{GM}}_{\mathrm{tot}}({r}^{2}+{r}_{{\rm h}}^{2})({r}^{2}+{r}_{{\rm c}}^{2})}{\pi {r}_{{\rm h}}}}\\ & & \times \left[{\displaystyle \frac{{\pi }^{2}}{8}}+{\displaystyle \frac{\pi }{6{r}^{3}{r}_{{\rm h}}}}+{\displaystyle \frac{1}{{r}^{2}{r}_{{\rm h}}^{2}}}-{\displaystyle \frac{\pi }{2{{rr}}_{{\rm h}}^{3}}}\right.\end{array}\end{eqnarray*}$$

$$\begin{eqnarray} & & -{\displaystyle \frac{{\mathrm{arctan}}^{2}{\displaystyle \frac{r}{{r}_{{\rm h}}}}}{2{r}_{{\rm h}}^{4}}}-{\displaystyle \frac{\mathrm{arctan}{\displaystyle \frac{{r}_{{\rm h}}}{r}}}{3{r}^{2}{r}_{{\rm h}}}}+{\displaystyle \frac{\mathrm{arctan}{\displaystyle \frac{{r}_{{\rm h}}}{r}}}{{{rr}}_{{\rm h}}^{3}}}\\ & & -\left.{\displaystyle \frac{2\mathrm{ln}\left(1+{r}_{{\rm h}}^{2}/{r}^{2}\right)}{3{r}_{{\rm h}}^{4}}}\right].\end{eqnarray} \tag{ 31 }$$

The limit is nontrivial, but cancellations of first through fourth order terms give ${\sigma }_{\mathrm{far}}^{2}\approx 3{{GM}}_{\mathrm{tot}}/(5r)$ when $r\gg {r}_{{\rm h}}$. When ${r}_{{\rm c}}\ll {r}_{{\rm h}}$, ${\sigma }_{\mathrm{near}}(r)$ and ${\sigma }_{\mathrm{far}}(r)$ can be joined piecewise to approximate $\sigma (r)$ for all r.