  We report new, extremely precise, photopolarimetry of the rapidly-rotating A0
main-sequence star $\zeta$ Aql, covering the wavelength range $\sim$400--900nm,
which reveals a rotationally-induced signal. We model the polarimetry, together
with the flux distribution and line profiles, in the framework of Roche
geometry with $\omega$-model gravity darkening, to establish the stellar
parameters. An additional constraint is provided by TESS photometry, which
shows variability with a period, $P_{\rm phot}$, of 11.1 hr. Modelling based on
solid-body surface rotation gives rotation periods, $P_{\rm rot}$, that are in
only marginal agreement with this value. We compute new ESTER stellar-structure
models to predict horizontal surface velocity fields, which depart from
solid-body rotation at only the $\sim$2% level (consistent with a reasonably
strong empirical upper limit on differential rotation derived from the
line-profile analysis). These models bring the equatorial rotation period,
$P_{\rm rot,e}$, into agreement with $P_{\rm phot}$, without requiring any
'fine tuning' (for the Gaia parallax). We confirm that surface abundances are
significantly subsolar ($\mbox{[M/H]} \simeq -0.5$). The star's basic
parameters are established with reasonably good precision: ${M =
2.53\pm0.16\,\mbox{M}_\odot}$, $\log{L/\mbox{L}_\odot} = 1.72\pm0.02$, $R_{\rm
p} = 2.21\pm 0.02\,\mbox{R}_\odot$, $T_{\rm eff} = 9693 \pm 50~\mbox{K}$, $i =
85{^{+5}_{-7}}^\circ$, and $\omega/\omega_{\rm c} = 0.95\pm0.02$. Comparison
with single-star, solar-abundance stellar-evolution models incorporating
rotational effects shows excellent agreement (but somewhat poorer agreement for
models at $\mbox{[M/H]} \simeq -0.4$).
