I am doing a gravitational N-body simulation of a star cluster. I set up my initial conditions using the Plummer model, which has potential $$\Phi_P = -\frac{G M}{\sqrt{r^2 + a^2}},\tag{1}$$ where $M$ is the total cluster mass, and $a$ is a parameter which I have seen explained in 2 different ways: some sources describe it as a scaling parameter for the size of the cluster, whereas other sources motivate the inclusion of $a$ as a numerical trick to prevent divergence of the gravitational potential at small separations.
For the force between stars, I am using a softened gravitational potential $$\Phi_G = -\frac{G m_1 m_2}{\sqrt{r^2 + \epsilon^2}},\tag{2}$$ where $\epsilon$ is motivated as a numerical trick for preventing divergence of the gravitational potential at small separations (note the similarity with the second definition of $a$).
My question boils down to this: should $a$ and $\epsilon$ have the same value? In case the answer is context dependent, I am interested in the studying the relaxation time of clusters with fewer than 1000 stars.
So far, I have been treating them as different values, because I understand $a$ as just specifying the size of the system and $\epsilon$ as something I needed to set appropriately in order to suppress close encounters that lead to the ejection of stars via unphysically large forces. A problem with this is that when I generate the Plummer model, it is meant to be born in Virial equilibrium, but if I then use a softened form of gravity, this modifies the potential (without modifying the kinetic energy), so that the system is not in fact in Virial equilibrium. Also, I am finding it difficult to decide on an appropriate value of $\epsilon$, and thought that perhaps the appropriate value is $a$ since the second explanation given for $a$ is identical to the explanation of $\epsilon$.
