I'll offer a few partial answers, which may eventually lead to a complete answer.

As observed in Saussol's paper (Theorem 3, Kac's lemma), if you have an ergodic invariant measure $\mu$ then the mean return time to a set $A$ is equal to $1/\mu(A)$. (Note, however, that this is the mean return time for trajectories that begin in $A$, rather than for an arbitrary initial condition, for which you'd need the hitting time).

So we need to know what the ergodic invariant measure $\mu$ is. Many topological dynamical systems have lots of invariant measures lying around, in which case it's a non-trivial task to pick the one you're really interested in. However, if we assume that your convex polygon satisfies the Veech dichotomy (see, for example, this article by John Smillie and Barak Weiss), then it's true that the flow in every direction is either periodic or uniquely ergodic. Regular polygons satisfy the Veech dichotomy, so in that setting your non-periodicity condition would be enough to guarantee that Lebesgue measure is the only invariant measure.

Of course, Lebesgue measure is always invariant for the billiard flow, so the main thing we need this for is the fact that Lebesgue measure is ergodic under these assumptions. (There may be an easier way to get ergodicity, but this is the first thing that came to my mind.)

If we write $X_\theta$ for the phase space of the billiard flow with angle $\theta$, and let $\theta_1,\dots,\theta_n$ be the angles that your $n$ particles are launched at, then the above discussion implies that the billiard flow on each $(X_{\theta_n},\text{Leb})$ is ergodic. The phase space for your entire system is $\prod_n X_{\theta_n}$, and the direct product of the $n$ different Lebesgue measures is an invariant measure for your overall system.

What we'd like to do next is say that this measure is actually ergodic. Unfortunately, it's not quite that simple, since the direct product of two ergodic systems need not be ergodic (just consider $R_\alpha\times R_\alpha$, the direct product of two circle rotations by irrational multiples of $\pi$). So to say that Lebesgue measure is ergodic for your whole system requires something more, which is where your condition that the trajectories be distinct ought to come in. I'm not sure exactly how this step should go, but you should be able to get something along the lines of "for almost every set of angles $(\theta_1,\theta_2,\dots,\theta_n)$, Lebesgue measure is ergodic for the whole system". Billiard flows have something to do with interval exchange transformations (IETs), and so this paper by Jon Chaika may well have the result that's needed at this point.

Once you know that (product) Lebesgue measure is invariant for the whole system, you're in the clear: given a ball of radius $\epsilon$, the (normalised) Lebesgue measure of that ball in each $X_\theta$ is $\pi \epsilon^2 / C$, where $C$ is the area of the polygon, and so the set in $\prod_n X_{\theta_n}$ corresponding to those configurations for which all particles lie in this ball is just $(\pi \epsilon^2 / C)^n$. The inverse of this is your expected return time, provided we have appropriate hypotheses to justify all the above steps. (Expected hitting time might be a different story, I'd need to think a bit.)

I'll offer a few partial answers, which may eventually lead to a complete answer.

As observed in Saussol's paper (Theorem 3, Kac's lemma), if you have an ergodic invariant measure $\mu$ then the mean return time to a set $A$ is equal to $1/\mu(A)$. (Note, however, that this is the mean return time for trajectories that begin in $A$, rather than for an arbitrary initial condition, for which you'd need the hitting time).

So we need to know what the ergodic invariant measure $\mu$ is. Many topological dynamical systems have lots of invariant measures lying around, in which case it's a non-trivial task to pick the one you're really interested in. However, if we assume that your convex polygon satisfies the Veech dichotomy (see, for example, this article by John Smillie and Barak Weiss), then it's true that the flow in every direction is either periodic or uniquely ergodic. Regular polygons satisfy the Veech dichotomy, so in that setting your non-periodicity condition would be enough to guarantee that Lebesgue measure is the only invariant measure.

Of course, Lebesgue measure is always invariant for the billiard flow, so the main thing we need this for is the fact that Lebesgue measure is ergodic under these assumptions. (There may be an easier way to get ergodicity, but this is the first thing that came to my mind.)

If we write $X_\theta$ for the phase space of the billiard flow with angle $\theta$, and let $\theta_1,\dots,\theta_n$ be the angles that your $n$ particles are launched at, then the above discussion implies that the billiard flow on each $(X_{\theta_n},\text{Leb})$ is ergodic. The phase space for your entire system is $\prod_n X_{\theta_n}$, and the direct product of the $n$ different Lebesgue measures is an invariant measure for your overall system.

What we'd like to do next is say that this measure is actually ergodic. Unfortunately, it's not quite that simple, since the direct product of two ergodic systems need not be ergodic (just consider $R_\alpha\times R_\alpha$, the direct product of two circle rotations by irrational multiples of $\pi$). So to say that Lebesgue measure is ergodic for your whole system requires something more, which is where your condition that the trajectories be distinct ought to come in. I'm not sure exactly how this step should go, but you should be able to get something along the lines of "for almost every set of angles $(\theta_1,\theta_2,\dots,\theta_n)$, Lebesgue measure is ergodic for the whole system". Billiard flows have something to do with interval exchange transformations (IETs), and so this paper by Jon Chaika may well have the result that's needed at this point.

Once you know that (product) Lebesgue measure is invariant for the whole system, you're in the clear: given a ball of radius $\epsilon$, the (normalised) Lebesgue measure of that ball in each $X_\theta$ is $\pi \epsilon^2 / C$, where $C$ is the area of the polygon, and so the set in $\prod_n X_{\theta_n}$ corresponding to those configurations for which all particles lie in this ball is just $(\pi \epsilon^2 / C)^n$. The inverse of this is your expected return time, provided we have appropriate hypotheses to justify all the above steps. (Expected hitting time might be a different story, I'd need to think a bit.)