On encounter rates in star clusters

Close encounters between stars in star forming regions are important as they can perturb or destroy protoplanetary discs, young planetary systems, and stellar multiple systems. We simulate simple, viralised, equal-mass $N$-body star clusters and find that both the rate and total number of encounters between stars varies by factors of several in statistically identical clusters due to the stochastic/chaotic details of orbits and stellar dynamics. Encounters tend to rapidly `saturate' in the core of a cluster, with stars there each having many encounters, while more distant stars have none. However, we find that the fraction of stars that have had at least one encounter within a particular distance grows in the same way (scaling with crossing time and half-mass radius) in all clusters, and we present a new (empirical) way of estimating the fraction of stars that have had at least one encounter at a particular distance.


INTRODUCTION
Young stars are commonly found with circumstellar discs (e.g. Hillenbrand et al. 1998;Lada et al. 2000;Haisch et al. 2000Haisch et al. , 2001, and these discs are thought to be where planet formation occurs. Since most stars are formed in relatively dense environments (e.g. Lada & Lada 2003), it is possible for the discs, and the on going planet formation process within, to be affected by close encounters between stars.
Simulations have shown that the effect of tidal perturbation from a stellar fly-by can range from slightly changing the density distribution in the disc to truncating or even destroying it (e.g. Clarke & Pringle 1993;Cuello et al. 2022), depending on how close the encounter is. This dynamical truncation, as well as photoevaporation (e.g. Concha-Ramírez et al. 2022), and face-on accretion (Wijnen et al. 2017), can significantly affect the population of young stars with discs in the early stages of star formation. Perturbations can also trigger disc instabilities (e.g. Thies et al. 2005Thies et al. , 2010 and may determine the population of planets forming in the disc (Ndugu et al. 2022). Another interesting effect of encounters on the disc is the misalignment between the rotational planes of the disc and the host star due to a non-coplanar encounter (e.g. Heller 1993;Larwood 1997). Encounters may alter already formed planetary systems, changing orbits (e.g. Breslau & Pfalzner 2019), or disrupting them (e.g. Parker & Quanz 2012). And encounters can similarly alter or destroy multiple stellar systems (e.g. Goodwin 2010;Reipurth et al. 2014).
Young stars with masses 1 M typically have discs with radii of a few hundreds of au (e.g. Andrews & Williams 2007). For the discs of those stars to be significantly perturbed in encounters, the periastron distance between the encountering stars needs to be less than ∼ 1000 au. Therefore, to understand how important encounters are in affecting discs/planet formation, a key question is how many young stars have encounters within 1000 au?
There are two approaches one might take to finding the rates and numbers of encounters in some star cluster of interest. The first is to perform N -body simulations of a variety of similar systems which is time consuming and computationally expensive (e.g. Parker & Quanz 2012;Craig & Krumholz 2013), the second would be to have some (ideally analytic) estimate to quickly get at least a 'feel' for the expected values.
In this paper, we examine the encounter rate in a number of N -body simulations of bound star clusters. We show that encounter rates can vary by up to an order of magnitude between statistically identical clusters, but the fraction of stars that have had an encounter remains statistically the same. We present an empirical way of estimating the fraction of stars in a cluster that have had at least one encounter within a particular distance.

THE ENCOUNTER RATE
The number of encounters per unit time (ε) for a star seems like it should depend on several factors: the encounter distance of interest, some average number density of stars, and the typical velocity of stars. The velocity of stars will affect the encounter rate by both changing how often stars encounter other stars, and also changing how effective gravitational focusing is.

The standard method
The most common method of calculating encounter times is based on the fundamental assumptions that a star is travelling through an effectively infinite, uniform density medium at a constant speed (see e.g. the derivation in Binney & Tremaine 2008, note that these assumptions are perfectly adequate if one is interested in e.g. the Galactic disc).
The encounter rate for any individual star is typically given by where n is the number density of the stars, σ is the velocity dispersion, r e is the closest distance during the encounter, G is the gravitational constant, and m is the typical mass of the stars (Binney & Tremaine 2008). The second term in the brackets is associated with the gravitational focusing effect which deflects the trajectories and decreases the distance of closest approach for slow encounters or encounters between more massive stars. For an ensemble of N stars (e.g. a cluster), it seems reasonable to assume that the total encounter rate (total number of encounters per unit time, E) scales with the total number of stars. Since each encounter involves two stars, the encounter rate should scale with N/2, i.e. E N ε/2.
In convenient units where r e is in au, σ in km s −1 , n in pc −3 , and m in M , the encounter rate in a cluster of N stars is then E = 8.5 × 10 −11 N nσ r 2 e + 886 m σ 2 r e Myr −1 .
(2) Therefore, it would seem that to calculate the rate of encounters, E, at a particular distance of interest, r e , in an ensemble of N stars, the correct values of (a) number density, n, and (b) velocity dispersion, σ, are required. If there is a distribution of stellar masses, an appropriate value for m must be taken.
While various assumptions that go into this simple calculation are clearly wrong for star clusters (e.g. moving through an effectively infinite uniform density medium at a constant speed) one might think that some simple variation on this approach might work (e.g. taking some appropriate average speed and density). However, we will show that this approach in star clusters gives an often wrong, and an always misleading, 'answer'.

What do we want to know?
It is important to clarify what we want to know about encounters in a cluster. In most cases, what we would like is an estimate of what fraction of stars have had a close encounter as this tells us the relative levels of disc/planetary system/multiple system perturbation/destruction. It is important to remember that this is not what an estimate of an encounter rate gives without a further assumption of how the encounters are distributed between stars.
As an example let us take a cluster that we shall examine in detail later: an N = 300, M = 300 M , equalmass (so m = 1 M ) virialised Plummer sphere cluster with half-mass radius r h = 0.5 pc. If we want to know the number of stars that have had an encounter within, e.g. r e = 1000 au, we can calculate that E ∼ 25 Myr −1 by taking the values for the velocity dispersion and half-mass density of this cluster and putting them into equation (2).
If we assume this encounter rate estimate is correct, to calculate how many stars have had an encounter after some time we need to make a further assumption that encounters are random so that after 2 Myr there will be 50 stars that have had an encounter, and after 10 Myr 250 stars (ie. > 80 per cent) will have had an encounter. (One could be somewhat more sophisticated and estimate as the encounter fraction starts to approach unity how many stars have had zero, one, two etc. encounters.) In assuming encounters are random, this calculation ignores that encounters are much more likely to occur in the core, and that after a few crossing times some stars in the core are likely to have had multiple encounters, Table 1. The properties of cluster ensembles. Each ensemble contains 10 clusters with different random number seeds (labelled A to J) with the same number of stars (N ), stellar masses (m), half-mass radii (r h ), and crossing time (tcr). The ID of a simulation contains information on the initial conditions: N3 or N6 are N = 300 and N = 600 respectively, SM stands for single-mass, and after R is the initial half mass radius of the cluster.

Cluster IDs
while those in the halo may have had none. Indeed, when stated like this, this approach does seem extremely naive and it would be surprising if it gave the correct answer.

SIMULATIONS
We investigate encounter rates in clusters by performing N -body simulations in which we can record individual encounters, the stars involved, and their distances of closest approach. The simulations we report here are of the simplest bound systems: virialised Plummer spheres of equal-mass stars.
We simulate N = 300 and N = 600 virialised Plummer spheres (Plummer 1911) initialised by the method described in Aarseth et al. (1974) with initial half-mass radii of 0.5, 0.75 and 1 pc. Simulations are run only with equal-mass stars so that we can ignore any complicating effects of mass spectra.
The number of stars (N ), the stellar mass (m), the initial half-mass radius (r h ), and the label of simulated clusters are shown in Table 1. Clusters with r h = 0.5 pc are truncated at 3 pc, those with r h = 0.75 pc are truncated at 5 pc, and r h = 1 pc clusters are truncated at 7 pc so that they have approximately the same relative sizes.
All simulations are run for 10 Myr using our own Nbody code. The code uses a forth-order Hermite scheme (Makino & Aarseth 1992) as the integrator. We keep the energy error well below 10 −4 by employing an adaptive timestep, i.e. using equation (7) from Makino & Aarseth (1992) with parameter η = 4×10 −4 for N = 300 clusters and η = 1 × 10 −4 for N = 600 clusters. We also use block timestepping for N = 600 runs to speed up the calculations.
The separation between any pair of stars is monitored at every timestep. Once two stars are closer to each other than 1000 au, whether they are bound or unbound, they are considered as having close encounter. During this period, the closest separation is recorded and taken as the encounter distance once the stars move away from each other beyond 1000 au. In binaries multiple 'encounters' will occur, if the separation stays below 1000 au this is only counted as one 'encounter'. This is to prevent hard binaries inflating the close encounter rate, however as we discuss below hard binaries form extremely rarely in our simulations.
Simulations are run with a gravitational softening length of 0.01 au to avoid collisions or computationally expensive very close encounters. This is only of importance for the details of extremely close encounters at 10 au which is much closer than the vast majority of encounters and at a distance that would completely disrupt discs or planetary systems.

Number density and velocity dispersion
It would seem reasonable to assume that encounter rates and fractions should depend in some way on both number density and velocity dispersion. Below we go into some detail on how the number density and velocity dispersion of a cluster might be quantified. This is important to show that the encounter rates we measure in simulations disagree with simple calculations because the assumptions that underlie them are wrong for clusters, rather than that we are using the 'wrong' values for number density or velocity dispersion, or not accounting for how they change with time. A reader happy to take our word for this can skip the details below.

Number density
The derivation of equation (2) assumes that stars are uniformly distributed and so n is constant in time and space. This is a reasonable assumption for e.g. encounters in the Galactic disc, but not for encounters in a cluster (or any region where the number density varies on short length-scales).
There are a number of ways one could quantify some average number density, and they can result in very different values for the estimated encounter rate.
The first average number density we use is simply calculated from the half-mass radius of the cluster (r h ): The second average number density is the mean number density defined by where n and f are the number density and the probability density function of the distance of stars from the centre of mass of the cluster (r). In practice, we can approximate the integral by where ∆r and N bin are the size and the number of bins of stellar distances from the centre of mass of the cluster. Theoretically, the probability density function is defined as f = dP/dr, where dP is the probability of finding stars at distance between r and r+dr from the centre of mass of the cluster. But practically, the value of f i in equation (5) may be obtained from where N i is the number of stars in the i th -bin and N is the number of stars in the cluster. The number density n i in equation (5) is related to f i via where ∆V is a spherical volume element containing ∆N stars. Substituting equations (6) and (7) in (5) gives The half-mass number density is often used as it is simple to calculate, the more complex mean number density has the advantage of including information on the full density distribution of the cluster. We also note here that the half-mass radius is often used as a characteristic radius as it remains roughly constant over the long-term evolution of a cluster (Aarseth et al. 1974), however the half-mass radius does fluctuate, especially at early times and so even the half-mass density changes (sometimes by factors of several).

Velocity dispersion
The velocity dispersion (σ) in equation (1) and (2) comes from the assumption that the velocity distribution of the stars in the cluster is Maxwellian. From the Maxwell-Boltzmann distribution, the velocity dispersion is related to the mode of the distribution (v m ) by σ = v m / √ 2. The mode can simply be determined by constructing the velocity histogram and then fitting it with a polynomial regression to find the position of the peak.
It should be noted that this way of finding the velocity dispersion is only possible when all (3D) velocities are well known. In any observation the value of σ is either 'guessed' by assuming virial equilibrium, or observed in either 1D (radial velocities) or 2D (proper motions) with usually quite significant errors and biases (e.g. binary inflation, see Cottaar et al. (2012)).

Time-averaging
There are two ways to calculate the number density and velocity dispersion to use in equation (2). One is to take the values of n or σ calculated instantaneously at the end of the simulation, the other is to take a time average. In a simulation it is possible to calculate full 3D timeaveraged values for various quantities. However, to estimate encounter rates in an observed region or for a single snapshot only the current values for any quantity can be calculated, and even they might be uncertain/guesstimated (e.g. no velocity data is available and only 2D positions).
For later reference, Table 2 shows the initial, timeaveraged and final (i.e. those at 10 Myr) values of σ, n h , and n m for all simulations.

RESULTS
In our simulations we follow all encounters at distances of < 1000 au. We record when the encounter occurred, which two stars were involved, and the distance of closest approach. This allows us to find the encounter rate within a particular distance (i.e. what equation (2) attempts to estimate), and the number of stars that have had such an encounter -something equation (2) cannot tell us without further assumptions, but is often what we wish to know. 4.1. Comparing an N = 300 and an N = 600 cluster We start by comparing encounter rates at various distances in two clusters with N = 300 and N = 600 equal-mass stars. Both are initially virialised Plummer Spheres with r h = 0.5 pc, with stars each of mass 1 M . Figure 1 shows the cumulative number of encounters over 10 Myr in the N = 300 (top panel) and N = 600 (bottom panel) clusters. In each panel the lines from bottom-to-top are the cumulative numbers of encounters at r e < 50 (orange), 100 (green), 500 (magenta), and 1000 (blue) au respectively. At the top left of each sub-figure are three numbers for each of the encounter distances: the first is the actual encounter rate (Myr −1 ) as found in the simulation, the next two are time-averaged estimates that we will return to later, but note for now that all three numbers are often quite differ-ent. The two simulations shown are N3SMR050-A (top) and N6SMR050-A (bottom).  1. The total number of encounters grows roughly linearly with time (in these two clusters at least).
2. The number of encounters at different distances scales very roughly with r 2 e (e.g. for N = 300 after 10 Myr there have been 317 encounters at r e < 500 au, and 27 at < 50 au).
3. Increasing both N (and therefore also n) by a factor of two results in about 4 times more encounters (e.g. 317 when N = 300, and 1165 when N = 600 at r e < 500 au).
The second and third numbers in the top left are the estimates of encounter rate as calculated from equation (2) using the time-averaged values of the half-mass number density (n h ) and the mean number density (n m ) respectively.
In both cases using n h under-estimates the number of encounters by a factor of ∼ 2. Using n m seems better, often giving a reasonable estimate (but sometimes being off by a factor of ∼ 2). So at a first glance at just these two simulations one might consider that using n m often provides a reasonable estimate of the encounter rate in a cluster.
However, we have only compared two simulations which just happened to be those labelled A in our ensembles. As we show below, when we look at the whole ensemble the picture becomes much more complicated, and this emphasises the importance of looking at ensembles of simulations when dealing with N -body systems.

An ensemble of statistically identical clusters
We now consider all ten clusters in our N = 300 equalmass stars, and a half-mass radius of r h = 0.5 pc ensemble. The only difference between these clusters is the random number seed used to generate the initial positions and velocities. Therefore one would expect that the encounter rates in each would be similar -and ideally be close to an analytic estimate.
The top panel of Fig. 2 shows the final encounter rates (E sim ) for each of our identical clusters measured after 10 Myr from r e = 0 to 1000 au. 1 The simulation shown in the top panel of Fig. 1 is A which is the black line towards the middle-bottom of all the lines.
The most important thing to note about the top panel of Fig. 2 is the total encounter rates after 10 Myr varies very significantly between clusters with a difference of almost an order of magnitude. There appears to be no 'typical' clusters and some outliers -just a seemingly Table 2.
The first three triple columns are the initial (Ini.), time averaged (Avg.) and final (End) values of the velocity dispersion (σ), the half-mass number density (n h ) and the mean number density (nm) of clusters N3SMR050-A..J (initial conditions in Table 1). In the last triple column are the analytic estimates of the encounter rates at re < 1000 au, using the average half-mass number density (Eest(n h )), and the average mean number density (Eest(nm)), compared with the actual values of the encounter rate measured in the simulations (Esim). Note that the encounter rate can also be calculated from the initial or final values of σ, n h and nm.
Cluster random spread in encounter rates between clusters that are initially statistically identical. Interestingly, all of the curves in the top panel of Fig.  2 follow a distribution that goes roughly as r 2 e suggesting that the distribution of encounter distances is what would be expected for unbound encounters. However, the distribution can slightly deviate from r 2 e , often due to three-body encounters between a single star and a binary (which has formed during the simulation). These encounters can cause an increase in the encounter rate at r e ∼ 1000 au, as can be most obviously seen in cluster I (dark red line at the top of the figure).
The other panels of Fig. 2 show the difference between the analytically estimated encounter rates and the actual encounter rates for each cluster. In the middle panel the half-mass density (n h ) is used for the estimate, and in the bottom panel it is the mean number density (n m ). A good match to the analytic estimate is the black dashed line at a ratio of unity.
Using the half-mass density never gives a good estimate, and can be wrong by an order of magnitude. Using the mean density is slightly better -three or four clusters stay reasonably close to unity, but most clusters are always wrong by factors of several.
One obvious explanation would be that different clusters have changed significantly (some expanding and some contracting?).
In Table 2 we give the initial, time averaged and final values of σ, n h , and n m for all ten clusters, as well as the final cumulative encounter rate. The averages and variance of each quantity over all the clusters are given in the bottom line.
The velocity dispersions (σ) are very similar between all clusters, but measures of density can change significantly over time with quite different time averaged and final values, and between different measures (half mass or mean).
However, these variations do not seem to correlate with the vastly different encounter rates. The last three columns show the encounter rates estimated with n h and n m and then the actual encounter rates from the simulations. Only once (cluster F) does n h get close to predicting the actual encounter rate. Estimates using n m are reasonable for 5/10 of the clusters, but very wrong for 5/10. We think this was seen by Craig & Krumholz (2013) who note that their simulations had far more encounters than one might expect, but substructure complicated their analysis.
What is particularly interesting is that there is no systematic change in encounters rates with any measure of density or how they evolve. Cluster F has the lowest final mean density (484 pc −1 ) and the lowest encounter rate (30 Myr −1 ), but cluster I has the second lowest final mean density (531 pc −1 ), but the highest encounter rate (208 Myr −1 ). Clusters B and E have almost the same encounter rate, but final mean densities that are different by a factor of over two (735 and 1813 pc −1 ).
One might think that maybe there was some extreme deviation in density at some point in time that the time averaged densities do not properly include, however this is not the case. In Fig. 3 we show the numbers of encounters with r e < 1000 au (top panel), half-mass density (middle panel), and mean density (lower panel) for clusters D and I (cf. Fig. 1).
Cluster D (black line) shows a relatively linear increase in encounter numbers to end with ∼ 1400 encounters within 10 Myr. Cluster I (red line) is similar to cluster D until a period between 6-7 Myr when the encounter rate increases significantly.
There is no reason to think that the increased encounter rate in cluster I is due to density variations however. The middle panel shows that both cluster's half-mass densities are very similar, and both fairly constant and declining slightly. The bottom panel shows more variation in the mean densities with short-lived fluctuations of factors of a few, but both clusters show this behaviour, and, if anything, cluster D has higher densities. There are fluctuations in the mean density of cluster I around when the encounter rate increases significantly, but there are others when it does not.
Examination of the data shows that the large numbers of encounters in cluster I at 6-7 Myr is due to a few pairs of stars having multiple self-encounters in weakly bound pairs (cf. Moeckel & Clarke 2011).
In all the ensembles of statistically identical clusters we have run we find no systematic relationship between any measure of cluster density and encounter rates (apart from occasionally in just a few of the clusters, but these could be chance given that many fluctuations do not correlate).

Binaries
All our simulations start with no binaries. An interesting question is how many binaries can form, and how they might alter the evolution.
Soft binaries are extremely easy to form (see Moeckel & Clarke 2011), and can inflate the encounter rate. Any wide binary with periastron below 1000au and apastron above 1000au will be included as multiple encountershowever such binaries are extremely soft and short lived in our simulations (this was seen in cluster I).
Hard binaries, however, are much more difficult to form. The soft binaries we find are very weakly bound and can appear due to fluctuations in the global potential. However, to form a hard, long-lived, binary system requires a three-body encounter as the third body is needed to carry-away the excess energy (see Goodman & Hut 1993). A back-of-the-envelope calculation suggests hard binary formation should be rare in our clusters, and an examination of the simulations finds only a few hard binaries have managed to form (one every few simulations, and never more than one in a simulation).

The number of stars having had an encounter
We clearly see that the total number of encounters between stars at any particular distance can be different by maybe an order of magnitude in initially statistically identical clusters, in a way that cannot be explained by density fluctuations.
The encounter rates we have shown in Fig. 2 and Table 2 are the number of times two stars come closer together than a particular distance. However, this measure does not include information on if a particular star, or a particular pair of stars, have had multiple encounters.
Star clusters have a density distribution with a high density core and increasingly lower density as one moves outwards, and a Plummer profile is a reasonable approximation to young, bound clusters. Within this density distribution stars can have a variety of orbits (which can change after encounters). Some stars will spend a significant amount of their time in the high density core, some will spend most of their time in the low density halo, and various combinations in between (orbits can be radial or circular etc. and can change over time).
This means that each individual star will have a unique encounter history. Those that spend a lot of time in the core may have multiple encounters, while those in the halo may have none. In addition (as seen above), some stars may get into loosely bound multiples (cf. Moeckel & Clarke 2011) and potentially have numerous self-encounters which can inflate the encounter rate significantly (see above).
Despite the large variation in encounter rates, When we measure the encounter fraction in each of our ten clusters we find that this value is statistically the same.
In Table 3 we show the number (N s ) and fraction (f s ) of stars in each of the ten statistically identical clusters from Fig. 2 and Table 2 that have had an encounter within 1000 au after 10 Myr. This number is between 143 and 170 (157±9) -statistically consistent with being the same number, and a little over half the stars in the cluster at f s = 0.58 ± 0.03. This tells us that in all clusters in this ensemble the same fraction of stars are having very different numbers of encounters.
We can also examine other ensembles of statistically identical clusters and we find that the encounter fraction in different clusters in the same ensemble is statistically the same.
The mean and variances of encounter fractions for each ensemble are given in Table 4, but to summarise for r e < 1000 au: for N = 300 clusters with r h = 0.75pc, the encounter fraction is 0.38 ± 0.02; for N = 300 clusters with r h = 1pc, the encounter fraction is 0.28 ± 0.03; and for N = 600 clusters with r h = 0.5pc, the encounter fraction is 0.60 ± 0.02.
We do not present the data here in detail, but the same is true for different encounter distances: ie. the encounter fraction is lower when the distance is e.g. 500au, but the encounter fraction is statistically the same within each ensemble. (It is difficult to say anything about extremely close encounters as we are into small-N statistics.) That the encounter fraction is constant is extremely interesting, as the most useful measure of encounters is often how many stars have had at least one encounter closer than a particular distance over a particular timescale.

Encounter fractions in different clusters
As we saw above, within an ensemble of statistically identical clusters the fraction of stars that have at least one encounter within 1000 au within 10 Myr is statistically the same, but it is different between different ensembles.
The top panel of Fig. 4 shows how the encounter fraction, f s increases with (absolute) time for all of our ensembles. The blue line and shaded region which shows the variance at the top are for the N = 600 clusters with r h = 0.5 pc. The red line and shaded region are for the N = 300 clusters with r h = 0.5 pc. The purple line and shaded region are for the N = 300 clusters with r h = 0.75 pc. And at the bottom the green line and shaded region are for the N = 300 clusters with r h = 1 pc.
As can be seen, in each case the encounter fraction within ensembles evolves in the same general way -ris- Table 3. For each cluster in the N = 300, r h = 0.5 pc ensemble (column 1) with N stars (column 2), we show the half-mass radius r h (column 3), crossing time tcr (column 4), and the total number Ns (column 5) and fraction fs (column 6) of stars that have had an encounter within 1000 au in 10 Myr. Half mass radii and crossing times are calculated explicitly for each cluster. The final row shows the means and variances of each quantity over the ensemble. ing rapidly and then flattening -but different ensembles seem to evolve at different rates.
That encounters occur at different rates in these different ensembles should not be a surprise as each of the clusters have a different internal dynamical timescale set by their crossing time. In the middle panel of Fig. 4 we show the encounter fractions by crossing time, rather than by physical time, and the differences between the different ensembles becomes less pronounced. That the (green) N = 300 clusters with r h = 1 pc have had the fewest encounters is clearly to a large extent because these clusters are dynamically much younger.
However, it is clearly not just dynamical age that is important as the lines are still somewhat different. The reason for this is that the cluster size plays a role. In all of these simulations we are counting encounters within 1000 au which is a more significant fraction of the distances between stars in an r h = 0.5 pc cluster than in an r h = 1 pc cluster. So, we would expect the encounter timescale to also be sensitive to the relative impact parameter (r h /r e ) 2 .
In the bottom panel of Fig. 4, we plot encounter fraction against an 'encounter crossing time', t * cr , defined as t * cr = t cr r h /pc r e /1000 au 2 . ( Now in the bottom panel we appear to have found a timescale on which all clusters show extremely similar behaviour. For encounters within 1000 au there is a sharp rise in f s in the first 5 t * cr to a point where roughly a third of all stars have had an encounter. It then takes another ∼ 50 t * cr for the next third of stars to have an encounter. To test the scaling with (r h /r e ) 2 , in Fig. 5 we compare encounters within 500 au in an r h = 0.5 pc cluster (red) with encounters within 1000 au in an r h = 1 pc cluster (green). Here (r h /r e ) 2 is the same in both clusters (half the encounter distance, but half the half-mass radius), therefore we would expect the growth of the different encounter fractions with just the crossing time to be the same, which they are as is clear in the figure.
Note that there are various subtleties at play such as different velocity dispersions causing the effect of gravitational focusing to be different, and we are dealing with small-N stochastic systems. But overall, the overlap between evolution between different clusters within ensembles and very different ensembles is impressive when scaled by crossing time and relative encounter cross section. Figure 5. The encounter fractions against crossing time with re < 1000 au for N = 300, r h = 1 pc (green shaded region), and with re < 500 au for N = 300, r h = 0.5 pc (red shaded region).

An empirical relationship
The bottom panel of Fig. 4 provides a way of getting a rough estimate of the encounter fraction of stars, f s , at a particular encounter distance, r e , after some time, t, in any cluster if one knows the crossing time, t cr and half-mass radius, r h . Then from equation (9) one can calculate t * cr , and so t/t * cr . It is worth pointing-out that the curve in the bottom panel looks like it should have a fairly simple function that would provide a fit. However, we have struggled to find a simple (2 or 3 parameter) function that fits (the problem is that the initial rise is much steeper than e.g. an exponential will fit). Therefore we would suggest simply reading-off the value of f s from the bottom panel of Fig. 4 for whatever value of t/t * cr . We stress that this is a rough estimate, however, as rough-and-ready as this may be, it will still almost certainly provide a much better feeling for how many stars have had an encounter than any attempt to use equation (2), and then to extrapolate to an encounter fraction.

An example
We can take a roughly Orion Nebula-like cluster with M = 1000 M , N = 2500, a half-mass radius r h = 0.7 pc, and age 3 Myr and attempt to estimate what fraction of stars have had an encounter at < 1000 au. For such a cluster, n = 1800 pc −3 , and assuming virial equilibrium σ = 2.6 km s −1 , and so t cr = 0.27 Myr.
From equation (9) if r e = 1000 au, then t * cr = 0.5t cr = 0.14 Myr. Therefore this cluster has an 'encounter age' of t/t * cr ∼ 20. From the bottom panel of Fig. 4 this would suggest around 40 per cent of stars will have had an encounter within 1000 au.
If we use equation (2), we find E ∼ 1000 Myr −1 . For an age of 3 Myr this suggests 3000 encounters among the N = 2500 stars. An extremely naive extrapolation might suggest that therefore all stars had an encounter within 1000 au. One can be a little more sophisticated and assume that encounters are random which finds that ∼ 90 per cent of stars will have had an encounter (as each one of the 3000 encounters involves two stars, even if they are random most stars will have been involved in one). Even if the value of 3000 encounters in 3 Myr happened, by luck, to be right, the extension of this encounter number to the number of stars involved in encounters is certainly not random.

CONCLUSIONS
We have performed N -body simulations of small star clusters to investigate stellar encounters with separations r e < 1000 au. This is the regime in which discs, planetary systems, and multiple stellar systems can be significantly perturbed or destroyed.
We find that the encounter rates vary by up to an order of magnitude between statistically identical clusters. However, we find that the fraction of stars that have had an encounter is statistically the same within statistically identical clusters.
The fraction of stars that have had an encounter increases rapidly at early dynamical times before flattening significantly once stars in orbits particularly susceptible to encounters have had at least one encounter. This depends on both the dynamical timescale of the cluster (t cr ), and the relative impact parameter (r h /r e ) 2 .
We find a consistent, and reasonably tight, relationship between the fraction of stars that have had an encounter and a modified crossing time t * cr ∝ t cr (r h /r e ) 2 . The relationship we have found has a seemingly solid physical basis, but no detailed theoretical underpinning (we are working on this). However, it provides a simple way of extracting an estimate of the encounter fraction for a particular cluster of a particular age from a figure. While this is empirical, it almost certainly provides a much better estimate of the true encounter fraction than any attempt to apply standard theory.