Show that $\lim_{N\to\infty}\int_0^1 \left|\frac1N\sum_{n=1}^N f(x+\xi_n) \right|^{\,2}\, dx = 0$ if $\int_0^1 f(x)\, dx = 0$ and $f$ is periodic

I'll do what I think is the harder part, and leave some details to you. Note that Riemann-integrability is actually a red-herring; I think they added that hypothesis in because up to this point in the text, they only talked about the Riemann integral. Anyways, since you said you know some measure theory, let us prove a slightly more general claim (but not more difficult).

Let $\mathcal{F}$ be the set of all Lebesgue measurable functions $f:\Bbb{R}\to\Bbb{C}$ which are periodic of period $1$, and such that $\int_0^1|f(x)|^2\,dx<\infty$ and $\int_0^1f=0$ (this latter requirement makes sense since $f\in L^2([0,1])$ implies $f\in L^1([0,1])$ by Cauchy-Schwarz/Holder). Now, for each $N\in\Bbb{N}$, define the operator $T_N:\mathcal{F}\to\mathcal{F}$ by \begin{align} (T_Nf)(x):=\frac{1}{N}\sum_{n=1}^Nf(x+\xi_n) \end{align} You can clearly see that $T_Nf$ is also periodic with period $1$, and has finite squared integral over $[0,1]$, which means $T_N$ does indeed have image lying inside $\mathcal{F}$, so it's a well-defined linear mapping. Now, I claim that $T_N$ is actually bounded, with respect to the $L^2([0,1])$ norm. We have \begin{align} \|T_Nf\|_2^2 &:=\int_0^1\left|\frac{1}{N}\sum_{n=1}^Nf(x+\xi_n)\right|^2\,dx\\ &\leq \int_0^1\frac{1}{N}\sum_{n=1}^N\left|f(x+\xi_n)\right|^2\,dx&\tag{Jensen's inequality}\\ &=\frac{1}{N}\sum_{n=1}^N\int_{\xi_n}^{1+\xi_n}|f(y)|^2\,dy\\ &\leq \frac{1}{N}\sum_{n=1}^N\int_0^2|f(y)|^2\,dy\\ &=2\int_0^1|f(y)|^2\,dy\tag{periodicity}\\ &=2\|f\|_2^2 \end{align} Hence, we have $\|T_Nf\|_2\leq \sqrt{2}\|f\|_2$; the explicit constant $\sqrt{2}$ isn't really important though. Note that in my application of Jensen's inequality, what I used is that $\phi:\Bbb{R}\to\Bbb{R}$, $\phi(t)=t^2$ is a convex function, and $\frac{1}{N}\sum_{n=1}^Na_n$ can be viewed as integration with respect to a probability measure (namely the counting measure on $\{1,\dots, N\}$ divided by $N$). Notice also that so far we didn't use the fact that for functions $f\in\mathcal{F}$, we have $\int_0^1f=0$; this will be useful in the next step.

Next, what happens is a typical analysis idea. Fix a collection of functions $\mathcal{G}\subset \mathcal{F}$. Then, for any $g\in \mathcal{G}$, any $f\in\mathcal{F}$, and any $N\in\Bbb{N}$, we have \begin{align} \|T_Nf\|_2&=\|T_Nf-T_Ng+T_Ng\|_2\\ &\leq \|T_N(f-g)\|_2+\|T_Ng\|_2\\ &\leq \sqrt{2}\|f-g\|_2+\|T_Ng\|_2,\tag{$*$} \end{align} where in the middle we used linearity of $T_N$ and the fact that $\|\cdot\|_2$ is a norm, hence satisfies the triangle inequality. From this estimate, it is almost clear how to finish. We should choose a collection $\mathcal{G}$ of "good" functions in $\mathcal{F}$, for which we already know the conclusion, i.e that for every $g\in\mathcal{G}$, $\lim\limits_{N\to \infty}\|T_Ng\|_2=0$; and we should also ensure that this collection of functions is big enough, i.e that it is dense in $\mathcal{F}$, so that we can conclude that for all $f\in\mathcal{F}$, $\lim\limits_{N\to\infty}\|T_Nf\|_2=0$ as well.

Such density arguments are very standard and important in analysis, so I leave it to you to carefully justify how the inequality $(*)$ allows us to conclude the result for all $f\in\mathcal{F}$.

All that remains is choosing a good class $\mathcal{G}$ of "good functions". My first instinct would have been to choose the collection of trigonometric polynomials. But since you said you've done part (a), you can take $\mathcal{G}$ to be the set of continuous $1$-periodic functions $g:\Bbb{R}\to\Bbb{C}$ such that $\int_0^1g=0$. For such functions, you already know that $\lim\limits_{N\to\infty}\|T_Ng\|_2=0$, so all you have to prove is that they are dense in $\mathcal{F}$ (with respect to the $L^2([0,1])$ norm). I leave this bit to you.

I hope you realize that at no point would Riemann-integrability simplify things (ok maybe using Riemann-integrability, one could explicitly show that $\mathcal{G}$ is dense in $\mathcal{F}\cap \text{set of Riemann-integrable functions on $[0,1]$}$, as opposed to invoking some measure-theoretic results).

You say

In the referenced post, $f$ was continuous - so we could approximate it by trigonometric polynomials and use Weyl's criterion of equidistribution. However, now f is merely Riemann integrable, so we cannot do the same thing.

Well, yes and no. Continuity of $f$ allowed you to get a uniform approximation by trigonometric polynomials, and hence use Weyl's criterion on each of the exponentials there. But for this part of the question, we don't need a uniform approximation (it would be nice if we could, but you're right that we can't). But that doesn't matter because the question isn't asking about uniformity; all we need to do is approximate in the $L^2$ integral norm. Note that more generally, the assertion holds for any function belonging to the closure of $\mathcal{G}$ with respect to any of the $L^p([0,1])$ norms.