real analysis - Dependence of error on mesh for Riemann sums

If $f(x) = x$ on $[0,1]$, then a mesh less than $delta$ means the Riemann sum differs from the integral by about $delta/2$, so $delta(epsilon)$ is about $2epsilon$. Thus, I presume that you mean $epsilon$ instead of $epsilon^2$.

A function with bounded variation satisfies a linear upper bound for the error in terms of the mesh size and total variation, hence a linear lower bound for the required $delta(epsilon)$.

A function without bounded variation may require a smaller mesh. For example, consider the continuous extension of
$f(x) = x cos(pi/x^n)$ to $[0,1]$. The difference in $f$ between sampling $f$ at $x=(2k+1)^{-1/n}$ and at
$x=(2k)^{-1/n}$ is about $2x$ at locations differing by about $1/n~x^{n+1}$. That means with mesh $delta$, you can sample $f$ to be either $1$ or $-1$ up to about $x=c delta^{1/{(n+1)}}$, for an error of $cdelta^{2/(n+1)}$.

Therefore, $delta(epsilon)$ can't be larger than $c epsilon^{(n+1)/2}$.