One heuristic is to replace the $n^{th}$ roots of unity by $n$ iid elements $\zeta_1,\dots,\zeta_n$ of the unit circle, drawn uniformly at random. For any sum $\zeta_{i_1} + \dots + \zeta_{i_k}$ of $k$ of these with $i_1 < \dots < i_k$, a standard Fourier analytic calculation (Esseen concentration inequality) shows that if $k \geq 5$ (to make the $k$-fold convolution power of the Fourier transform of the unit circle have absolutely integrable Fourier transform), then
$$ {\bf P}( |\zeta_{i_1} + \dots + \zeta_{i_k}| < r ) \ll_k r^2 $$
for $0 < r < 1$. Taking the union bound over all $\binom{n}{k}$ different sums, we see that the probability that $k$ distinct elements of the $\zeta_1,\dots,\zeta_n$ sum to something of magnitude less than $r$ is $O_k( \binom{n}{k} r^2 )$. This shows that with probability $\gg 1$, we should have the lower bound
$$ |\zeta_{i_1} + \dots + \zeta_{i_k}| \gg_k \binom{n}{k}^{-1/2} \asymp_k n^{-k/2}$$
for all $i_1 < \dots < i_k$. This is for a fixed $n$ (which after averaging corresponds to predicting a lower bound $\gg_k n^{-k/2}$ for a positive density fraction of $n$, rather than for all $n$). To get a bound that is true almost surely for all but finitely many $n$ (in the spirit of the "strong" law of large numbers, as opposed to the "weak" law), the Borel-Cantelli lemma suggests that we need to get the failure probability for a fixed $n$ down to something like $1/n^{1+\varepsilon}$, and so the lower bound now worsens from $\gg_k n^{-k/2}$ to $\gg_{k} n^{-(k+1+\varepsilon)/2}$.
Some remarks:
- The gap between the $k/2$ type exponents for the lower bound and the $k/4$ type exponent for the upper bound is related to the square root present in the birthday paradox.
- These bounds are fairly reversible in this model due to all the independence, and with a little more effort one should be able to show that the bounds of $\gg_k n^{-k/2}$ (for a positive density set of $n$) and $\gg_k n^{-(k+1+\varepsilon)/2}$ (for all but finitely many $n$) cannot be significantly improved except possibly for the $\varepsilon$ factor, but I have not attempted to work this out rigorously.
- To complete the heuristic analysis one should also treat the cases where there are collisions amongst the $i_1,\dots,i_k$, which can be handled by the same methods as long as there are at least $5$ distinct values of the $i_1,\dots,i_k$; the cases of only four or fewer values need to be treated by more direct computations but I think they should give lower order contributions in the $k \geq 5$ regime.
EDIT: while the analysis of the above random model remains accurate within the scope of that model, I have just realised that the rotation invariance symmetry of the $n^{th}$ roots of unity leads to a larger lower bound prediction than the above naive model, as it indicates that the roots of unity have some relevant structure to them that prevents them from behaving like completely random elements of the unit circle. The key point is that for the original $n^{th}$ roots of unity problem, one can always rotate one of the roots to equal $1$. So one should really be looking at sums of the form $1 + \zeta_{i_2} + \dots + \zeta_{i_k}$ with $1 < i_2 < \dots < i_k$. This effectively lowers $\binom{n}{k}$ to $\binom{n-1}{k-1}$ and raises all the exponents in the previous analysis by $1/2$, thus the revised prediction would be a lower bound of $\gg_k n^{-(k-1)/2}$ for a positive fraction of $n$ and $\gg_k n^{-(k+\varepsilon)/2}$ for all but finitely many $n$. There are some additional symmetries coming from cyclic permutation of the roots and conjugation symmetry but these only affect the $k$-dependent implied constants and not the exponents. In principle, the further Galois group symmetries of the problem could lead to additional refined corrections to the model, but I tentatively am of the opinion that this will not change the predicted exponents further.