In this work, we calculate the eigenvalues of the probe (overlap) Dirac operator on thermal gauge ensembles of $2+1$ flavor QCD generated using domain wall fermions as well as pure $SU(3)$ gauge theory on the lattice. Focusing on the infrared part of the eigenspectrum that lies within the non-perturbative magnetic scale, we propose suitable observables that allow us to categorize different regions of the eigenspectrum unambiguously. While most of these eigenmodes are completely delocalized and chaotic, i.e. their nearest-neighbor level spacing fluctuations are similar to random matrices of a Gaussian unitary ensemble (GUE), we showed that a classical non-thermal state of SU(3) gauge theory consisting of magnetic gluons is also chaotic, thus a non-trivial realization of the Bohigas-Giannoni-Schmit conjecture. This allowed us to estimate an upper bound on the thermalization time $\sim 1.44$
fm/c of magnetic gluons by matching the magnetic scales in these two regimes. Furthermore we also observe the appearance of a few eigenmodes in deep-infrared part of the spectrum near and above the chiral crossover temperature, whose fractal dimensions might carry information about the universality class of the chiral transition in QCD.