We study the number $P(n)$ of partitions of an integer $n$ into sums of distinct squares and derive an integral representation of the function $P(n)$. Using semi-classical and quantum statistical methods, we determine its asymptotic average part $P_{as}(n)$, deriving higher-order contributions to the known leading-order expression [M. Tran {it et al.}, Ann. Phys. (N.Y.) {bf 311}, 204 (2004)], which yield a faster convergence to the average values of the exact $P(n)$. From the Fourier spectrum of $P(n)$ we obtain hints that integer-valued frequencies belonging to the smallest Pythagorean triples $(m,p,q)$ of integers with $m^2+p^2=q^2$ play an important role in the oscillations of $P(n)$. Finally we analyze the oscillating part $delta P(n)=P(n)-P_{as}(n)$ in the spirit of semi-classical periodic orbit theory [M. Brack and R. K. Bhaduri: {it Semiclassical Physics} (Bolder, Westview Press, 2003)]. A semi-classical trace formula is derived which accurately reproduces the exact $delta P(n)$ for $n > sim 500$ using 10 pairs of `orbits. For $n > sim 4000$ only two pairs of orbits with the frequencies 4 and 5 -- belonging to the lowest Pythagorean triple (3,4,5) -- are relevant and create the prominent beating pattern in the oscillations. For $n > sim 100,000$ the beat fades away and the oscillations are given by just one pair of orbits with frequency 4.