Let $X\sim\text{Normal}(0,1)$ and let $f_X$ be its probability density function. I conducted some numerical experiments in the software Mathematica to estimate $f_X$ via a kernel method. Let $\hat{f}_X^M$ be the kernel density estimate using a sample of length $M$. Let
$$\epsilon=E\left[\|f_X-\hat{f}_X^M\|_\infty\right]$$
be the error ($E$ is the expectation). I noticed that the error decreases with $M$ until a certain length $M_0$ from which the error stabilizes. For example, in Mathematica, the built-in function SmoothKernelDistribution
employs the Gaussian kernel with Silverman's rule to determine the bandwidth by default. In the following figure in log-log scale, I show the error $\epsilon$ for different values of $M$ growing geometrically, where the expectation that defines $\epsilon$ is estimated using 20 realizations of $\|f_X-\hat{f}_X^M\|_\infty$. I also plot the estimated $90\%$ confidence interval for $\|f_X-\hat{f}_X^M\|_\infty$ (dashed lines).
Observe that the error decreases linearly in log-log scale (that is, at rate $O(M^{-r})$), up to a certain length $M$ where it starts to stabilize. Also, the confidence intervals become more narrow in the end. Is this issue due to accumulated numerical errors? Is it due to Silverman's rule?