  We investigate the Dirac eigenvalue spectrum ($\rho(\lambda,m_l)$) to study
the microscopic origin of axial anomaly in high temperature phase of QCD. We
propose novel relations between the derivatives ($\partial^n
\rho(\lambda,m_l)/\partial m_l^n$) of the Dirac eigenvalue spectrum with
respect to the quark mass ($m_l$) and the $(n+1)$-point correlations among the
eigenvalues ($\lambda$) of the massless Dirac operator. Based on these
relations, we present lattice QCD results for $\partial^n
\rho(\lambda,m_l)/\partial m_l^n$ ($n=1, 2, 3$) with $m_l$ corresponding to
pion masses $m_\pi=160-55$ MeV, and at a temperature of about 1.6 times the
chiral phase transition temperature. Calculations were carried out using
(2+1)-flavors of highly improved staggered quarks and the tree-level Symanzik
gauge action with the physical strange quark mass, three lattice spacings
$a=0.12, 0.08, 0.06$ fm, and lattices having aspect ratios $4-9$. We find that
$\rho(\lambda\to0,m_l)$ develops a peaked structure. This peaked structure,
which arises due to non-Poisson correlations within the infrared part of the
Dirac eigenvalue spectrum, becomes sharper as $a\to0$, and its amplitude is
proportional to $m_l^2$. After continuum and chiral extrapolations, we find
that the axial anomaly remains manifested in two-point correlation functions of
scalar and pseudo-scalar mesons in the chiral limit. We demonstrate that the
behavior of $\rho(\lambda\to0,m_l)$ is responsible for it.
