Let $S(t) = {1}/{\pi} \arg \zeta \big(\hh + it \big)$ be the argument of the Riemann zeta-function at the point 1/2 + it. For n ⩾ 1 and t > 0 define its iterates
$$\begin{equation*}
S_n(t) = \int_0^t S_{n-1}(\tau) \,\d\tau\, + \delta_n\,,
\end{equation*}$$
where δn is a specific constant depending on n and S0(t) ≔ S(t). In 1924, J. E. Littlewood proved, under the Riemann hypothesis (RH), that Sn(t) = O(log t/(log log t)n + 1). The order of magnitude of this estimate was never improved up to this date. The best bounds for S(t) and S1(t) are currently due to Carneiro, Chandee and Milinovich. In this paper we establish, under RH, an explicit form of this estimate
$$\begin{equation*}
-\left( C^-_n + o(1)\right) \frac{\log t}{(\log \log t)^{n+1}} \ \leq \ S_n(t) \ \leq \ \left( C^+_n + o(1)\right) \frac{\log t}{(\log \log t)^{n+1}}\,,
\end{equation*}$$
for all n ⩾ 2, with the constants C±n decaying exponentially fast as n → ∞. This improves (for all n ⩾ 2) a result of Wakasa, who had previously obtained such bounds with constants tending to a stationary value when n → ∞. Our method uses special extremal functions of exponential type derived from the Gaussian subordination framework of Carneiro, Littmann and Vaaler for the cases when n is odd, and an optimized interpolation argument for the cases when n is even. In the final section we extend these results to a general class of L-functions.