  We study numerically the onset of chaos and thermalization in the
Banks-Fischler-Shenker-Susskind (BFSS) matrix model with and without fermions,
considering Lyapunov exponents, entanglement generation, and quasinormal
ringing. We approximate the real-time dynamics in terms of the most general
Gaussian density matrices with parameters which obey self-consistent equations
of motion, thus extending the applicability of real-time simulations beyond the
classical limit. Initial values of these Gaussian density matrices are
optimized to be as close as possible to the thermal equilibrium state of the
system. Thus attempting to bridge between the low-energy regime with a
calculable holographic description and the classical regime at high energies,
we find that quantum corrections to classical dynamics tend to decrease the
Lyapunov exponents, which is essential for consistency with the
Maldacena-Shenker-Stanford (MSS) bound at low temperatures. The entanglement
entropy is found to exhibit an expected "scrambling" behavior - rapid initial
growth followed by saturation. At least at high temperatures the entanglement
saturation time appears to be governed by classical Lyapunov exponents. Decay
of quasinormal modes is found to be characterized by the shortest time scale of
all. We also find that while the bosonic matrix model becomes non-chaotic in
the low-temperature regime, for the full BFSS model with fermions the leading
Lyapunov exponent, entanglement saturation time, and decay rate of quasinormal
modes all remain finite and non-zero down to the lowest temperatures.
