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 explore types of convergence which are not encapsulated by the 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)}{e\cdot \sqrt{n \log(\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$.

Note: The last integral on the right hand side is for me hard to crack. However numerical computations suggest that:
\begin{equation}
\lim_{k\rightarrow 0} \frac{1}{\log(\log(1/k))} \cdot \int\limits_0^\infty (-\imath k)^\delta \cdot \Gamma(-\delta,-\imath e k) d \delta  = 1
\end{equation}
Indeed by using the integral representation of the Gamma function along with integration by parts we quickly establish the following identity:
\begin{eqnarray}
(-\imath k)^\delta \cdot \Gamma(-\delta,-\imath e k) = \frac{e^{-\delta}}{\delta} + (-\imath k)^\delta \cdot \Gamma(-\delta) + \sum\limits_{n=1}^\infty \frac{(\imath k)^n}{n!}\cdot \frac{e^{n-\delta}}{\delta-n}
\end{eqnarray}

Now clearly 
\begin{eqnarray}
&&\int\limits_0^\infty (-\imath k)^\delta \cdot \Gamma(-\delta,-\imath e k) d \delta =\\
&& \int\limits_0^\infty \left( \frac{e^{-\delta}}{\delta} + (-\imath k)^\delta \cdot \Gamma(-\delta) \right) d\delta + O(k)\\
&&= \int\limits_0^\infty \left( \frac{e^{-\delta}}{\delta} - \frac{(-\imath k)^\delta}{\delta}  \right) d\delta +
\int\limits_0^\infty (-\imath k)^\delta \left(\Gamma(-\delta)+\frac{1}{\delta}\right) d\delta + O(k)\\
&&= \left.\left( Ei(-\delta) - Ei(-A \delta)\right)\right|_0^\infty+
\int\limits_0^\infty (-\imath k)^\delta \left(\Gamma(-\delta)+\frac{1}{\delta}\right) d\delta + O(k)\\
&&= \log(-A) + \int\limits_0^\infty (-\imath k)^\delta \left(\Gamma(-\delta)+\frac{1}{\delta}\right) d\delta + O(k)
\end{eqnarray}
where $A=-\log(-\imath k)= \imath \pi/2 - \log(k)$. Now we have checked numericaly that the integral in the middle above decays monotonically when $k\rightarrow 0$. Since now $\log(-A) \rightarrow \log(\log(1/k))$ when $k\rightarrow 0$ the claim is established.

Now we check that our test statistic is properly normalized.


Define $c_n:=\sqrt{n\log(\log(n))}$.
Indeed we have:
\begin{eqnarray}
\log\left( \kappa_{S_n}(k)\right) &=& -\imath k \frac{3}{2} \frac{n}{c_n} + n \log\left[ \kappa_{f^{-1}(X)}(\frac{k}{e c_n})\right] \\
&=& \frac{1+4 \log(2)-4 \log(2-2 \log(k)+\log(n) + \log(\log(\log(n))))}{8\log(\log(n))} k^2 + O(\frac{k^3}{\sqrt{n \cdot \log(\log(n))}})
\end{eqnarray}
where in the second line I simply used the Series[] command in Mathematica to obtain the expansion . Now we can clearly see that:
\begin{equation}
\lim_{n\rightarrow \infty} \log\left( \kappa_{S_n}(k)\right) = -\frac{k^2}{2}
\end{equation}
as it should be.

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)/(E Sqrt[n Log[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/yJNJB.jpg
  [2]: https://i.sstatic.net/r4Zf1.jpg
  [3]: https://i.sstatic.net/gExa5.jpg