Magnitude of Lennard-Jones force evaluated with molecular dynamics simulation
The following plots the force for reasonable distances

Intuitively, if you are implementing the verlet algorithm then the step size will be $\delta x \propto a\delta t^2$ (ignoring velocities). If you want a step size of order ~$f \times \sigma$ (where f is a scaling factor) then you'd like to satisfy a relation such as $\delta t &lt; \sqrt{\frac{f\sigma}{a}}$ (if ignoring velocities). The timestep cutoff $\delta t$ computed according to this equation with a force of $\pu{10^-12 N}$ and f=0.1 is 0.3 psec. But 0.3 psec is still a relatively long timestep. Clearly the problem is not with acceleration (provided the timestep is shorter than the cutoff just calculated), but maybe there is one with the velocities. You can proceed in the same way but solving the complete quadratic verlet equation for the timestep in terms of the step size, and inserting a value for the rms speed based on the desired T and step size, to arrive at a timestep cutoff. Using this reasoning I came up with a cutoff $\delta t &lt; \pu{ 15 fsec}$
One problem with a large step size (and discretization of the trajectory in general) is that the energy is not conserved and can easily blow up (along with the temperature): if you take a step of order $\sigma$ from the potential well minimum toward the maximum in the potential then the atoms will suddenly find themselves in the strongly repulsive high energy region.
You may want to check Ref. 1 which discusses this to some extent, and provides this advice:

One wants to find the largest time step that will maintain the conservation of energy.

Reference

Sangrak Kim. Issues on the Choice of a Proper Time Step in Molecular Dynamics. Physics  Procedia 53 (2014) 60–62.
doi: 10.1016/j.phpro.2014.06.027