It may be instructive to demonstrate this result from first principles and basic results, exploiting properties of cumulant generating functions (exactly as in standard proofs of the Central Limit Theorem). It requires us to understand the rate of growth of generalized harmonic numbers $$H(n,s)=\sum_{k=1}^n k^{-s}$$ for $s=1, 2, \ldots.$ These growth rates are well-known and easily obtained by comparison to the integrals $\int_1^n x^{-s}dx$: they converge for $s \gt 1$ and otherwise diverge logarithmically for $s=1$.
Let $n \ge 2$ and $1 \le k \le n$. By definition, the cumulant generating function (cgf) of $(X_k - 1/k)/B_n$ is
$$\psi_{k,n}(t) = \log\mathbb{E}\left(\exp\left(\frac{X_k - 1/k}{B_n}t\right)\right) = -\frac{t}{k B_n} + \log\left(1 + \frac{-1 + \exp(t/B_n)}{k}\right).$$
The series expansion of the right hand side, obtained from the expansion of $\log(1+z)$ around $z=0$, takes the form
$$\psi_{k,n}(t) = \frac{(k-1)}{2 k^2 B_n^2}t^2 + \frac{k^2 - 3k + 2}{6 k^3 B_n^3} t^3 + \cdots + \frac{k^{j-1} - \cdots \pm (j-1)!}{j! k^j B_n^j}t^j + \cdots.$$
The numerators of the fractions are polynomials in $k$ with leading term $k^{j-1}$. Because the log expansion converges absolutely for $\left|\frac{-1 + \exp(t/B_n)}{k}\right| \lt 1$, this expansion converges absolutely when
$$\left|\exp(t/B_n) - 1\right| \lt k.$$
(In case $k=1$ it converges everywhere.) For fixed $k$ and increasing values of $n$, the (obvious) divergence of $B_n$ implies the domain of absolute convergence grows arbitrarily large. Thus, for any fixed $t$ and sufficiently large $n$, this expansion converges absolutely.
For sufficiently large $n$, then, we may therefore sum the individual $\psi_{k,n}$ over $k$ term by term in powers of $t$ to obtain the cgf of $S_n/B_n$,
$$\psi_n(t) = \sum_{k=1}^n \psi_{k,n}(t) = \frac{1}{2}t^2 + \cdots + \frac{1}{B_n^j}\left(\sum_{k=1}^n \left(k^{-1} - \cdots \pm (j-1)!k^{-j}\right)\right)\frac{t^j}{j} + \cdots.$$
Taking the terms in the sums over $k$ one at a time requires us to evaluate expressions proportional to
$$b(s,j) = \frac{1}{B_n^j}\sum_{k=1}^n k^{-s}$$
for $j \ge 3$ and $s=1, 2, \ldots, j$. Using the asymptotics of generalized harmonic numbers mentioned in the introduction, it follows easily from
$$B_n^2 = H(n,1) - H(n,2) \sim \log(n)$$
that
$$b(1,j) \sim (\log(n))^{1-j/2}\to 0$$
and (for $s \gt 1$)
$$b(s,j) \sim (\log(n))^{-j/2}\to 0$$
as $n$ grows large. Consequently all terms in the expansion of $\psi_n(t)$ beyond $t^2$ converge to zero, whence $\psi_n(t)$ converges to $t^2/2$ for any value of $t$. Since convergence of the cgf implies convergence of the characteristic function, we conclude from the Levy Continuity Theorem that $S_n/B_n$ approaches a random variable whose cgf is $t^2/2$: that is the standard Normal variable, QED.
This analysis uncovers just how delicate the convergence is: whereas in many versions of the Central Limit Theorem the coefficient of $t^j$ is $O(n^{1-j/2})$ (for $j \ge 3$), here the coefficient is only $O(((\log(n))^{1-j/2})$: the convergence is much slower. In this sense the sequence of standardized variables "just barely" becomes Normal.
We can see this slow convergence in a series of simulations. The histograms display $10^5$ independent iterations for four values of $n$. The red curves are graphs of standard normal density functions for visual reference. Although there is evidently a gradual tendency towards normality, even at $n=1000$ (where $(\log(n))^{-1/2} \approx 0.38$ is still sizable) there remains appreciable non-normality, as evidenced in the skewness (equal to $0.35$ in this sample). (It is no surprise the skewness of this histogram is close to $(\log(n))^{-1/2}$, because that's precisely what the $t^3$ term in the cgf is.)

Here is the R
code for those who would like to experiment further.
set.seed(17)
par(mfrow=c(1,4))
n.iter <- 1e5
for(n in c(30, 100, 300, 1000)) {
B.n <- sqrt(sum(rev((((1:n)-1) / (1:n)^2))))
x <- matrix(rbinom(n*n.iter, 1, 1/(1:n)), nrow=n, byrow=FALSE)
z <- colSums(x - 1/(1:n)) / B.n
hist(z, main=paste("n =", n), freq=FALSE, ylim=c(0, 1/2))
curve(dnorm(x), add=TRUE, col="Red", lwd=2)
}