  We prove a variant of the multidimensional polynomial Szemer\'edi theorem of
Bergelson and Leibman where one replaces polynomial sequences with other sparse
sequences defined by functions that belong to some Hardy field and satisfy
certain growth conditions. We do this by studying the limiting behavior of the
corresponding multiple ergodic averages and obtaining a simple limit formula. A
consequence of this formula in topological dynamics shows denseness of certain
orbits when the iterates are restricted to suitably chosen sparse subsequences.
Another consequence is that every syndetic set of integers contains certain
non-shift invariant patterns, and every finite coloring of $\N$, with each
color class a syndetic set, contains certain polychromatic patterns, results
very particular to our non-polynomial setup.
