This is a follow-up question to https://mathoverflow.net/questions/296717/rate-of-convergence-of-frac1-sqrtn-ln-n-sum-k-1n-1-sqrtx-k-2n/296718#296718 . My motivation is to construct a statistic whose rate convergence to a Gaussian will be very slow and as such to formulate a  counterexample to Berry-Esseen' theorem. We therefore define the following statistic:
\begin{equation}
S_n := \frac{\left(\sum\limits_{k=1}^n f^{-1}(X_k) - \frac{3}{2} e \cdot n\right)}{\sqrt{n \log(n)}}
\end{equation}
where $X_k$ are i.i.d. uniformly distributed in $(0,1)$ and the function $f()$ is defined as follows:
\begin{equation}
f(x) := \frac{e^2}{2} \cdot \frac{1+\log(x)}{x^2 \log(x)^2} 1_{x \ge e}
\end{equation}
Now, the probability density of $f(X)$ is as follows:
\begin{eqnarray}
\rho_{f(X)}(z) &=& \int\limits_0^1 \delta(z - f^{-1}(x)) dx
=-\int\limits_{e}^\infty \delta(z-u) f^{'}(u) du= -f^{'}(z) 1_{z \ge e}\\
&=& \frac{e^2}{2} \cdot \frac{2+3 \log(z)+2 \log(z)^2}{z^3 \log(z)^3}1_{z \ge e}
\end{eqnarray}
From this we readily get the moments:
\begin{eqnarray}
E\left[ f^{-1}(X) \right] = \frac{3}{2} e\\
E\left[ (f^{-1}(X))^2 \right] = \infty
\end{eqnarray}
We also get the characteristic function. It reads:
\begin{eqnarray}
\kappa_{f^{-1}(X)}(k) = e^{\imath k e}+ \imath k \frac{e}{2} e^{\imath k e}-k^2 \frac{e^2}{2}\cdot \int\limits_0^\infty (-\imath k)^\delta \cdot \Gamma(-\delta,-\imath e k) d \delta
\end{eqnarray}
for $0<k<1$.

Now, I carried out a Monte Carlo simulation and computed the sample Cumulative Distribution Function (CDF) of our statistic and plotted it along with the CDF of a standardized Gaussian distribution with the former and the later being plotted in Blue and Purple respectively. Here I took $n=5,10,15$ and in each case I used $m=1000$ realizations. The figures are below:

[CDFs at $n=5$][1]

[CDFs at $n=10$][2]

[CDFs at $n=15$][3]

I have used the following Mathematica code to produce those figures:

    m = 1000; n = 15; delta = 1/10;
    bins = Table[-5 + delta/2 + j delta, {j, 1, (10 - delta)/delta}];
    limD = CDF[NormalDistribution[0, 1], bins];
    X = RandomReal[{0, 1}, {m, n}];
    x =.; {t0, Y} = 
     Timing[(x /. 
        Map[First[
           NSolve[(E^2 (1 + Log[x]))/(2 x^2 Log[x]^2) == # && x > E, x, 
            Reals]] &, X, {2}])];
    ll = (Total[#] & /@ Y - 3/2 E n)/Sqrt[n Log[n]];
    emp = EmpiricalDistribution[ll];
    DD = CDF[emp, bins];
    pl = ListPlot[Transpose[{bins, #}] & /@ {DD, limD}, ImageSize -> 800, 
       LabelStyle -> {15, FontFamily -> "Arial"}, 
       BaseStyle -> {15, FontFamily -> "Bold"}, 
       PlotLabel -> "n=" <> ToString[n]];
    Export["LimitBehavior1_n_" <> ToString[n] <> ".jpg", pl, "JPEG"];
    Import["LimitBehavior1_n_" <> ToString[n] <> ".jpg"]

Having said all this my question is the following. What is the rate of convergence of our statistic towards a Gaussian. To be specific we are asking  about the behavior of the supremum norm of the difference in CDFs for large values of $n$.

  [1]: https://i.sstatic.net/lMsKd.jpg
  [2]: https://i.sstatic.net/XBktf.jpg
  [3]: https://i.sstatic.net/18wQB.jpg