Thank you for your replies. Unfortunately, as so it often happens, the problem was between the keyboard and the chair. While writing the MWE for a single k I went over the equations once more and noticed that I had forgotten a single a^2 term from the equation for \rho_2 (which should’ve been \rho_2 = \frac{1}{2\rho_0}\left(0.25+\omega^2\rho_0^2+\frac{a^2\dot\rho_0^2}{4}\right)). With this correction the “artifact” disappeared and @stevengj was indeed correct, it was an instability of the system.
Notwithstanding, I edited the calculations slightly based on @apo383 's suggestions, mainly changing m_{Eff}^2 for \frac{m_{Eff}^2}{a^2} = -R\left(\xi - \frac{1}{6}\right) and \omega^2 for \frac{\omega^2}{a^2} = \frac{k^2}{a^2} + \frac{m_{Eff}^2}{a^2} and merging the equations for \ddot\rho_0 and \rho_2 to \ddot\rho_0 = \frac{1}{2\rho_0}\left(\dot\rho_0^2 + \frac{1}{a^2}\right) - H\dot\rho_0 -2\rho_0\frac{\omega^2}{a^2}. Hopefully they’ll help with the stability with larges values of t.
