Creating realistic world - star cluster?

Let's determine some basic parameters of the cluster.

You say you want the stars to be 1-2 light-years apart. This is a little tricky because star clusters don't necessarily have uniform densities. Globular clusters, for instance, have radial density profiles decreasing from some core density at $R=0$ to zero density somewhere further out. There are several density profiles used to model this. The Plummer model is easily the simplest, as analytic expressions exist for density and potential. The King model (see King (1966), the third in a series of papers) is more accurate, though the density must be solved for numerically, and the Michie model can be even more accurate (see e.g. this thesis), as it can account for dark matter.

I'm going to follow a modified version of the models of Bonatto & Bica (2014), who use a King profile (here applicable to open clusters). The fraction of the total number of stars $N$ within a radius $R$ is given by $$n(R)=\frac{x^2-4u\left(\sqrt{1-x^2}-1\right)+u^2\ln\left(1+x^2\right)}{u^2\ln u^2-(u-1)(3u-1)}\tag{1}$$ where $$x\equiv R/R_c,\quad u\equiv\sqrt{1+\left(R_t/R_c\right)^2}$$ with $R_t$ and $R_c$ being the tidal and core radii of the cluster; the cluster is truncated at $R=R_t$. The authors recommend a ratio $R_t/R_c\approx15$. Therefore, we just have to specify one of those two parameters (we also know that $u=226$). The tidal radius can be computed depending on the mass of the cluster and its distance from the center of the galaxy (see e.g. Wu et al. (2009)). That would require us to specify additional parameters, so I'll assume that the cluster is average, and give it a tidal radius of about 10 parsecs. Thus, the core radius is roughly 0.67 parsecs. This could be an overestimate; Piskunov et al. (2008), for instance, surveyed 236 open clusters in the Milky Way and found that the tidal radius peaked at around 6 parsecs. Still, 10 parsecs is by no means unreasonable. At any rate, I'd prefer a more massive cluster, so more stars may survive for longer periods of time.

Some quick Mathematica gives me

Clear[Rt, Rc]
x[R_, Rc_] := R/Rc
u[Rt_, Rc_] := Sqrt[1 + (Rt/Rc)^2]u[Rt_, Rc_] := Sqrt[1 + (Rt/Rc)^2]
n[R_, Rt_, Rc_] := ((x[R, Rc])^2 - 4*(u[Rt, Rc])*(Sqrt[1 + (x[R, Rc])^2] - 1) + (u[Rt, Rc])^2*Log[1 + (x[R, Rc])^2])/((u[Rt, Rc])^2*Log[(u[Rt, Rc])^2] - (u[Rt, Rc] - 1)*(3*u[Rt, Rc] - 1))
Plot[n[R, 10, 2/3], {R, 0, 10}, AxesLabel -> {"Radius (Parsecs)", "Number of stars"}]

   

The actual number of stars inside $R$ is simply $n$ multiplied by the number of stars in the cluster, $N$. We'll take $N\sim10^3$, a decently sized cluster. How do we now go from $n\cdot N$ to number density, $\mathcal{N}$ (just to make things more confusing)? Well, take the formula for the mass $M$ of a sphere inside a radius $R$ with radial density $\rho(r)$: $$M=\int_0^R 4\pi r^2\rho(r)\mathrm{d}r$$ Simply divide $M$ by the average mass of a star (which I'll get to later, but is unimportant for now) and we get $$n\cdot N=\int_0^R4\pi r^2\mathcal{N}(r)\mathrm{d}r$$ We just differentiate and divide: $$\mathcal{N}(r)=\frac{\mathrm{d}n}{\mathrm{d}r}N\frac{1}{4\pi r^2}\tag{2}$$ In the code, we add

P[R_, Rt_, Rc_] = D[n[R, Rt, Rc], R]*1000*(4*Pi*R^2)^(-1)P[R_, Rt_, Rc_] = D[n[R, Rt, Rc], R]*1000*(4*Pi*R^2)^(-1)
Plot[P[R, 10, 2/3], {R, 0, 10}, AxesLabel -> {"Radius (Parsecs)", "Number density"},ScalingFunctions->"Log"]

where I've replaced $\mathcal{N}$ with $P$, and we get, with a logarithmic $y$-axis,

Note that the number density goes to infinity as $R$ goes towards the center (which we can check by simply evaluating the limit of $\mathcal{N}(R)$ as $R\to0$). However, we know that $n(R_t)$ is finite (although $\lim_{R\to\infty} n(R)=\infty$, which is why we truncated it). I played around with Animate for $N$:

Animate[Plot[P[R, 10, 2/3, Nt], {R, 0, 10}, AxesLabel -> {"Radius (Parsecs)", "Number density"}, ScalingFunctions -> "Log"], {Nt, 10^2, 10^5}, AnimationRunning -> False]

and found that you can get extremely high number densities in the center, if you so desire, which makes sense.

To find the average separation between two stars, we calculate the mean interparticle distance, scaled correctly (i.e. the Wigner"“Seitz radius): $$r_s(R)=\left(\frac{3}{4\pi\mathcal{N}(R)}\right)^{1/3}\tag{3}$$ We add another section of code:

rs[R_, Rt_, Rc_] := (3/(4*Pi*(n[R, Rt, Rc])))^(1/3)
Plot[rs[R, 10, 2/3], {R, 0, 10}, AxesLabel -> {"Radius (parsecs)", "rs"}, ScalingFunctions -> "Log"]

and, plotting this function, get

We see that $r_s$ is always greater than about 2 light-years for $0\leq R\leq R_t$. Towards the center. Therefore, the systems you're interested in should be near the outer edge of the cluster.

Hold on a minute, though. This is a density distribution based entirely on a sphere of stars, with nothing else influencing them. We've guessed that this models an open cluster fairly well, but we also need to assume that the cluster, if left to its own devices, is unstable, and the stars will drift away. Let's try adding a central object, like you suggested - an intermediate-mass black hole (IMBH). This isn't terribly far-fetched; it's an idea that's been explored quite a lot in the case of globular clusters. In open clusters, we have problems because we need to explain how the black hole got there in the first place. However, we can ignore that for now.

If there's a central, massive object in a star cluster, the cluster should exhibit a Bahcall-Wolf cusp. Within a distance equal to one fifth of the black hole's sphere of influence, the density should obey a power law of the form $$\rho_{BW}(R)\propto R^{-7/4}$$ Number density, assuming a roughly homogeneous population (again, more on that later), should obey this power law: $$\mathcal{N}_{BW}(R)\propto R^{-7/4}\tag{4}$$ All we have to do is calculate the sphere of influence, find the number density there according to the King model we computed, and fit the two together.

The sphere of influence, according to one definition, is $$r_h=\frac{GM_{BH}}{\sigma_0}\tag{5}$$ where $G$ is the gravitational constant $M_{BH}$ is the mass of the black hole, and $\sigma_0$ is the central stellar velocity dispersion, which I'll choose to be about .81 km/s - modeled after observations by Geller et al. (2009) of M35, a reasonably analog to our cluster. I'm going to choose $M_{BH}\sim10^4M_\odot$. We then get $r_h\approx 6.25$ parsecs, and the Bahcall-Wolf cutoff radius is roughly 1.25 parsecs.

Let's calculate the number density at $R=\frac{1}{5}r_h$ from the King model. We get 17.44 stars per cubic parsec - quite large! If we then assume that the cusp stars there, we know that $$\mathcal{N}_{BW}(R)=17.44=A(1.25)^{-7/4}$$ for some constant $A$, which we find to be $A=25.77$. All we have to do now is to plot the two together in a piecewise number density function, defined as follows: $$\mathcal{N}(R)= \begin{cases} \frac{\mathrm{d}n}{\mathrm{d}R}N\frac{1}{4\pi R^2},\quad 0\leq R\leq1.25\\ 25.77R^{-7/4},\quad1.25\leq R\leq10\\ \end{cases}\tag{6}$$ Plotting this on with a logarithmic $y$-axis gives us

Here's the code I used to do the calculations and generate the plots (although you can play around with various parameters, like black hole mass and total cluster mass):

Clear[Rt, Rc]
Rt = 10(*Tidal radius*)
Rc = 2/3(*Core radius*)
x[R_, Rc_] := R/Rc
u[Rt_, Rc_] := Sqrt[1 + (Rt/Rc)^2]
n[R_, Rt_, Rc_] := ((x[R, Rc])^2 - 4*(u[Rt, Rc])*(Sqrt[1 + (x[R, Rc])^2] - 1) + (u[Rt, Rc])^2*Log[1 + (x[R, Rc])^2])/((u[Rt, Rc])^2*Log[(u[Rt, Rc])^2] - (u[Rt, Rc] - 1)*(3*u[Rt, Rc] - 1))(*Number of stars in King model*)
P[R_, Rt_, Rc_] = D[n[R, Rt, Rc], R]*1000*(4*Pi*R^2)^(-1)(*Number density from King model*)
rs[R_, Rt_, Rc_] := (3/(4*Pi*(n[R, Rt, Rc])))^(1/3)(*Wigner-Seitz radius*)
PBW[R_, Rt_, Rc_] := 25.77*(R)^(-7/4)(*Bahcall-Wolf cusp*)
np[R_, Rt_, Rc_] := Integrate[(1/1000)*4*Pi*x^2*PBW[x, Rt, Rc], {x, 0, 1.25}] + Integrate[(1/1000)*4*Pi*y^2*P[y, Rt, Rc], {y, 1.25, R}](*Total number of stars in cluster*)
Plot[n[R, 10, 2/3], {R, 0, 10}, AxesLabel -> {"Radius (Parsecs)", "Number of stars"}]
Plot[P[R, 10, 2/3], {R, 0, 10}, AxesLabel -> {"Radius (Parsecs)", "Number density"}, ScalingFunctions -> "Log"]
Plot[rs[R, 10, 2/3], {R, 0, 10}, AxesLabel -> {"Radius (Parsecs)", "rs"}, ScalingFunctions -> "Log"]
Plot[{Piecewise[{{P[R, 10, 2/3], 1.25 < R < 10}, {PBW[R, 10, 2/3], 0 < R < 1.25}}]}, {R, 0, 10}, AxesLabel -> {"Radius (Parsecs)", "Number density"}, PlotRange -> 10^3, ScalingFunctions -> "Log"](*Plots final density*)

Congratulations! You just built an open cluster!

Well, we built a simplified model of an open cluster, and we made more assumptions than I'd prefer:

• All the mass of the cluster consists of stars and the central black hole, which gravitationally dominates most of the cluster. In reality, there should be gas and dust in the cluster, because the stars formed from giant molecular clouds.
• The stellar population is homogeneous, i.e. all of the stars are the same. In general, this isn't the case. Stars that formed together might be similar, but often they're quite different. Their masses should be be distributed as predicted by some sort of initial mass function (IMF). This created no problems so far, but it might when you take into account things like mass segregation that are important in globular clusters. If you're really interested in fixing this, you can use something like the MASSCLEAN packages to generate star clusters via the King model and an IMF (the default is the Salpeter IMF). It's a pretty cool tool, and can save you a lot of work.
• A King model defines the cluster. I'd be more comfortable using a King model to describe a globular cluster, and I see no reason why you couldn't just have a young globular cluster - like M15, as you suggested - but I'm okay applying our version of it to our open cluster.

None of these, however, are fatal flaws in our model.

So, is this cluster going to be stable over long timescales? That, after all, is the question we're trying to answer. Well, to start, there are several things in our favor:

• The cluster contains $\sim10^3$ stars, which is a decent amount for an open cluster. I think that should make it more likely to end up with a sizable number of stars.
• The black hole is quite massive compared to the cluster - more massive than the rest of the cluster, in fact. I did this on purpose, because it means that the black hole's sphere of influence is over half the tidal radius of the entire cluster. Again, this should make it easier to retain more stars.

Let's do a more formal assessment (based on e.g. these notes). Let the mean mass of a star in the cluster be $m_m$. By the virial theorem, the speed $v$ of a star at $R\approx. R_t$ is given by $$v^2\approx\frac{GNm_m}{R_t}\tag{7}$$ The time it takes a star to travel $R_t$, the crossing time $t_c$, is $$t_c=R/v\tag{8}$$ and the relaxation time $t_r$, the time over which encounters with other stars begin to have an effect on the motion of a star, is $$t_r\approx\frac{0.1N}{\ln N}t_c\tag{9}$$ Furthermore, if a percent $\gamma$ of stars have a velocity above a certain critical value, then the evaporation time of the system is $$t_e=\frac{t_r}{\gamma}\approx t_r\times10^2\tag{10}$$ So, what critical parameters have changed now that we added the Bahcall-Wolf cusp? Well, we can safely assume that $m_m$ and $R_t$ are the same, as we specified them. It turns out, however, than $N$ has changed. Integrating over the piecewise density distribution gives $n(R_t)=0.89$, meaning that we have fewer stars in the cluster.

At the same time, the total mass of the cluster is much, much larger, since the black hole's mass is substantially greater than the mass of the stars. Therefore, $$v^2\approx\frac{GNm_m+M_{BH}}{R_t}$$ Having fewer stars and adding a black hole will lead to a much shorter relaxation time. Even if we extend the tidal radius (as would be expected), the black hole still won't help stars at the outside of the cluster. I recall one paper noting that in the case of a black hole with a mass somewhat less than that of the cluster, most of the stars in the outer regions of the cluster will behave just as if the cluster was perfectly described by a King model, with no central black hole.

At the center of the cluster, there will indeed be some significant effects - and since the number density is much greater there, a large fraction of the stars should be impacted. I'm reluctant to draw analogies between stellar clusters and galaxies, but stars near the center should be influenced just like stars near the supermassive black hole at the center of the Milky Way, Sagittarius A*, orbit the black hole. The analogy is actually fairly good: Most stars in a galaxy don't orbit the central supermassive black hole; they orbit the total mass of the galaxy. It's a self-gravitating system, not a set of bodies orbiting one massive one. At the very center, however, stars are clearly orbiting Sagittarius A*. The same holds true for our open cluster.

So, to summarize:

• Stars in the outer regions of the open cluster will likely leave the cluster at the same time or earlier than if there was no intermediate mass black hole at the center.
• Stars within the sphere of influence of the black hole, especially within the Bahcall-Wolf cusp, will be influenced strongly by the black hole and may directly orbit it.
• The cluster will largely "evaporate" on the same timescale (though likely a shorter one) than if it had no central black hole. However, I think the central group of stars will remain bound to the black hole. Eventually, this will be all that remains of the open cluster, if the orbits are stable.
• Whether or not life can survive on planets orbiting these stars is another question entirely, and I'll deal with that in a future edit of this answer, if you want.