We performed numerical simulations with the PLUTO code in order to analyze the non-linear behavior of the Kelvin-Helmholtz instability in non-magnetized relativistic fluids. The relevance of the instability at the cosmological QCD phase transition was explored using an equation of state based on lattice QCD results with the addition of leptons. The results of the simulations were compared with the theoretical predictions of the linearized theory. For small Mach numbers up to $M_s \sim 0.1$ we find that both results are in good agreement. However, for higher Mach numbers, non-linear effects are significant. In particular, many initial conditions that look stable according to the linear analysis are shown to be unstable according to the full calculation. Since according to lattice calculations the cosmological QCD transition is a smooth crossover, violent fluid motions are not expected. Thus, in order to assess the role of the Kelvin-Helmholtz instability at the QCD epoch, we focus on simulations with low shear velocity and use monochromatic as well as random perturbations to trigger the instability. We find that the Kelvin-Helmholtz instability can strongly amplify turbulence in the primordial plasma and as a consequence it may increase the amount of primordial gravitational radiation. Such turbulence may be relevant for the evolution of the Universe at later stages and may have an impact in the stochastic gravitational wave background.

We study the evolution of turbulence in the early universe at the QCD epoch using a state-of-the-art equation of state derived from lattice QCD simulations. Since the transition is a crossover we assume that temperature and velocity fluctuations were generated by some event in the previous history of the Universe and survive until the QCD epoch due to the extremely large Reynolds number of the primordial fluid. The fluid at the QCD epoch is assumed to be non-viscous, based on the fact that the viscosity per entropy density of the quark gluon plasma obtained from heavy-ion collision experiments at the RHIC and the LHC is extremely small. Our hydrodynamic simulations show that the velocity spectrum is very different from the Kolmogorov power law considered in studies of primordial turbulence that focus on first order phase transitions. This is due to the fact that there is no continuous injection of energy in the system and the viscosity of the fluid is negligible. Thus, as kinetic energy cascades from the larger to the smaller scales, a large amount of kinetic energy is accumulated at the smallest scales due to the lack of dissipation. We have obtained the spectrum of the gravitational radiation emitted by the motion of the fluid finding that, if typical velocity and temperature fluctuations have an amplitude $(\Delta v) /c \gtrsim 10^{-2}$ and/or $\Delta T/T_c \gtrsim 10^{-3}$, they would be detected by eLISA at frequencies larger than $\sim 10^{-4}$ Hz.