  For a full rank lattice $\Lambda \subset \mathbb{R}^d$ and $\mathbf{A} \in
\mathbb{R}^d$, consider $N_{d,0;\Lambda,\mathbf{A}}(\Sigma) = \#
([\Lambda+\mathbf{A}] \cap \Sigma \mathbb{B}^d) = \# \{\mathbf{k}\in \Lambda :
|\mathbf{k}+\mathbf{A}| \leq \Sigma \}$. Consider the iterated integrals \[
  N_{d,k+1;\Lambda,\mathbf{A}}(\Sigma) = \int_0^\Sigma
N_{d,k;\Lambda,\mathbf{A}}(\sigma) \,\mathrm{d} \sigma, \] for $k\in
\mathbb{N}$. After an elementary derivation via the Poisson summation formula
of the sharp large-$\Sigma$ asymptotics of $N_{3,k;\Lambda,\mathbf{A}}(\Sigma)$
for $k\geq 2$ (these having an $O(\Sigma)$ error term), we discuss how they are
encoded in the structure of the Fourier transform
$\mathcal{F}N_{3;\Lambda,\mathbf{A}}(\tau)$. The analysis is related to
H\"ormander's analysis of spectral Riesz means, as the iterated integrals above
are weighted spectral Riesz means for the simplest magnetic Schr\"odinger
operator on the flat $3$-torus. That the $N_{3,k;\Lambda,\mathbf{A}}(\Sigma)$
obey an asymptotic expansion to $O(\Sigma^2)$ is a special case of a general
result holding for all magnetic Schr\"odinger operators on all manifolds, and
the subleading polynomial corrections can be identified in terms of the Laurent
series of the half-wave trace at $\tau=0$. The improvement to $O(\Sigma)$ for
$k\geq 2$ follows from a bound on the growth rate of the half-wave trace at
late times.
