At first we notice that
$$\begin{align*} (1-\alpha)^{\frac{1}{2}}(1-\beta)^{\frac{1}{2}} & =\left(1-\frac{1}{2} \cdot \alpha-\frac{1}{8} \cdot \alpha^{2}-\cdots\right)\left(1-\frac{1}{2} \cdot \beta-\frac{1}{8} \cdot \beta^{2}-\cdots\right) \\ & =\sum_{k, \ell \geq 0} c_{k, \ell} \cdot \alpha^{k} \beta^{\ell} \quad \text{ for all } \alpha, \beta \in(0,1) \end{align*}$$
where $c_{0,0}=1$ and $c_{k, \ell}$ are certain coefficients.
For an indirect proof, we suppose that $x_{n}=\sqrt{(a^{n}-1)(b^{n}-1)} \in \mathbb{Z}$ for all positive integers $n$. Replacing $a$ by $a^{2}$ and $b$ by $b^{2}$ if necessary, we may assume that $a$ and $b$ are perfect squares, hence $\sqrt{a b}$ is an integer.
At first we shall assume that $a^{\mu} \neq b^{\nu}$ for all positive integers $\mu, \nu$. We have
$$\begin{equation*} x_{n}=(\sqrt{a b})^{n}\left(1-\frac{1}{a^{n}}\right)^{\frac{1}{2}}\left(1-\frac{1}{b^{n}}\right)^{\frac{1}{2}}=\sum_{k, \ell \geq 0} c_{k, \ell}\left(\frac{\sqrt{a b}}{a^{k} b^{\ell}}\right)^{n} . \end{equation*}$$
Choosing $k_{0}$ and $\ell_{0}$ such that $a^{k_{0}}>\sqrt{a b}, b^{\ell_{0}}>\sqrt{a b}$, we define the polynomial
$$P(x)=\prod_{k=0, \ell=0}^{k_{0}-1, \ell_{0}-1}\left(a^{k} b^{\ell} x-\sqrt{a b}\right)=: \sum_{i=0}^{k_{0} \cdot \ell_{0}} d_{i} x^{i}$$
with integer coefficients $d_{i}$. By our assumption, the zeros
$$\frac{\sqrt{a b}}{a^{k} b^{\ell}}, \quad k=0, \ldots, k_{0}-1, \quad \ell=0, \ldots, \ell_{0}-1,$$
of $P$ are pairwise distinct.
Furthermore, we consider the integer sequence
$$\begin{equation*} y_{n}=\sum_{i=0}^{k_{0} \cdot \ell_{0}} d_{i} x_{n+i}, \quad n=1,2, \ldots \end{equation*}$$
By the theory of linear recursions, we obtain
$$\begin{equation*} y_{n}=\sum_{\substack{k, \ell \geq 0 \\ k \geq k_{0} \text{ or } \ell \geq \ell_{0}}} e_{k, \ell}\left(\frac{\sqrt{a b}}{a^{k} b^{\ell}}\right)^{n}, \quad n=1,2, \ldots, \end{equation*}$$
with real numbers $e_{k, \ell}$. We have
$$|y_{n}| \leq \sum_{\substack{k, \ell \geq 0 \\ k \geq k_{0} \text{ or } \ell \geq \ell_{0}}}|e_{k, \ell}|\left(\frac{\sqrt{a b}}{a^{k} b^{\ell}}\right)^{n}=: M_{n} .$$
Because the series in (4) is obtained by a finite linear combination of the absolutely convergent series (1), we conclude that in particular $M_{1}<\infty$. Since
$$\frac{\sqrt{a b}}{a^{k} b^{\ell}} \leq \lambda:=\max \left\{\frac{\sqrt{a b}}{a^{k_{0}}}, \frac{\sqrt{a b}}{b^{\ell_{0}}}\right\} \quad \text{ for all } k, \ell \geq 0 \text{ such that } k \geq k_{0} \text{ or } \ell \geq \ell_{0},$$
we get the estimates $M_{n+1} \leq \lambda M_{n}, n=1,2, \ldots$. Our choice of $k_{0}$ and $\ell_{0}$ ensures $\lambda<1$, which implies $M_{n} \rightarrow 0$ and consequently $y_{n} \rightarrow 0$ as $n \rightarrow \infty$. It follows that $y_{n}=0$ for all sufficiently large $n$.
So, equation (3) reduces to $\sum_{i=0}^{k_{0} \cdot \ell_{0}} d_{i} x_{n+i}=0$.
Using the theory of linear recursions again, for sufficiently large $n$ we have
$$x_{n}=\sum_{k=0, \ell=0}^{k_{0}-1, \ell_{0}-1} f_{k, \ell}\left(\frac{\sqrt{a b}}{a^{k} b^{\ell}}\right)^{n}$$
for certain real numbers $f_{k, \ell}$.
Comparing with (2), we see that $f_{k, \ell}=c_{k, \ell}$ for all $k, \ell \geq 0$ with $k<k_{0}, \ell<\ell_{0}$, and $c_{k, \ell}=0$ if $k \geq k_{0}$ or $\ell \geq \ell_{0}$, since we assumed that $a^{\mu} \neq b^{\nu}$ for all positive integers $\mu, \nu$.
In view of (1), this means
$$\begin{equation*} (1-\alpha)^{\frac{1}{2}}(1-\beta)^{\frac{1}{2}}=\sum_{k=0, \ell=0}^{k_{0}-1, \ell_{0}-1} c_{k, \ell} \cdot \alpha^{k} \beta^{\ell} \end{equation*}$$
for all real numbers $\alpha, \beta \in(0,1)$. We choose $k^{*}<k_{0}$ maximal such that there is some $i$ with $c_{k^{*}, i} \neq 0$. Squaring (5) and comparing coefficients of $\alpha^{2 k^{*}} \beta^{2 i^{*}}$, where $i^{*}$ is maximal with $c_{k^{*}, i^{*}} \neq 0$, we see that $k^{*}=0$. This means that the right hand side of (5) is independent of $\alpha$, which is clearly impossible.
We are left with the case that $a^{\mu}=b^{\nu}$ for some positive integers $\mu$ and $\nu$. We may assume that $\mu$ and $\nu$ are relatively prime. Then there is some positive integer $c$ such that $a=c^{\nu}$ and $b=c^{\mu}$. Now starting with the expansion (2), i.e.,
$$x_{n}=\sum_{j \geq 0} g_{j}\left(\frac{\sqrt{c^{\mu+\nu}}}{c^{j}}\right)^{n}$$
for certain coefficients $g_{j}$, and repeating the arguments above, we see that $g_{j}=0$ for sufficiently large $j$, say $j>j_{0}$. But this means that
$$\left(1-x^{\mu}\right)^{\frac{1}{2}}\left(1-x^{\nu}\right)^{\frac{1}{2}}=\sum_{j=0}^{j_{0}} g_{j} x^{j}$$
for all real numbers $x \in(0,1)$. Squaring, we see that
$$\left(1-x^{\mu}\right)\left(1-x^{\nu}\right)$$
is the square of a polynomial in $x$. In particular, all its zeros are of order at least $2$, which implies $\mu=\nu$ by looking at roots of unity. So we obtain $\mu=\nu=1$, i.e., $a=b$, a contradiction.

We set $a^{2}=A, b^{2}=B$, and $z_{n}=\sqrt{(A^{n}-1)(B^{n}-1)}$. Let us assume that $z_{n}$ is an integer for $n=1,2, \ldots$ Without loss of generality, we may suppose that $b<a$. We determine an integer $k \geq 2$ such that $b^{k-1} \leq a<b^{k}$, and define a sequence $\gamma_{1}, \gamma_{2}, \ldots$ of rational numbers such that
$$2 \gamma_{1}=1 \quad \text{ and } \quad 2 \gamma_{n+1}=\sum_{i=1}^{n} \gamma_{i} \gamma_{n-i} \text{ for } n=1,2, \ldots$$
It could easily be shown that $\gamma_{n}=\frac{1 \cdot 1 \cdot 3 \ldots(2 n-3)}{2 \cdot 4 \cdot 6 \ldots 2 n}$, for instance by reading Vandermonde's convolution as an equation between polynomials, but we shall have no use for this fact.
Using Landau's $O$-Notation in the usual way, we have
$$\begin{aligned} & \left\{(a b)^{n}-\gamma_{1}\left(\frac{a}{b}\right)^{n}-\gamma_{2}\left(\frac{a}{b^{3}}\right)^{n}-\cdots-\gamma_{k}\left(\frac{a}{b^{2 k-1}}\right)^{n}+O\left(\frac{b}{a}\right)^{n}\right\}^{2} \\ & =A^{n} B^{n}-2 \gamma_{1} A^{n}-\sum_{i=2}^{k}\left(2 \gamma_{i}-\sum_{j=1}^{i-1} \gamma_{j} \gamma_{i-j}\right)\left(\frac{A}{B^{i-1}}\right)^{n}+O\left(\frac{A}{B^{k}}\right)^{n}+O\left(B^{n}\right) \\ & =A^{n} B^{n}-A^{n}+O\left(B^{n}\right) \end{aligned}$$
whence
$$z_{n}=(a b)^{n}-\gamma_{1}\left(\frac{a}{b}\right)^{n}-\gamma_{2}\left(\frac{a}{b^{3}}\right)^{n}-\cdots-\gamma_{k}\left(\frac{a}{b^{2 k-1}}\right)^{n}+O\left(\frac{b}{a}\right)^{n} .$$
Now choose rational numbers $r_{1}, r_{2}, \ldots, r_{k+1}$ such that
$$(x-a b) \cdot\left(x-\frac{a}{b}\right) \ldots\left(x-\frac{a}{b^{2 k-1}}\right)=x^{k+1}-r_{1} x^{k}+\cdots \pm r_{k+1},$$
and then a natural number $M$ for which $M r_{1}, M r_{2}, \ldots M r_{k+1}$ are integers. For known reasons,
$$M\left(z_{n+k+1}-r_{1} z_{n+k}+\cdots \pm r_{k+1} z_{n}\right)=O\left(\frac{b}{a}\right)^{n}$$
for all $n \in \mathbb{N}$ and thus there is a natural number $N$ which is so large, that
$$z_{n+k+1}=r_{1} z_{n+k}-r_{2} z_{n+k-1}+\cdots \mp r_{k+1} z_{n}$$
holds for all $n \geqslant N$. Now the theory of linear recursions reveals that there are some rational numbers $\delta_{0}, \delta_{1}, \delta_{2}, \ldots, \delta_{k}$ such that
$$z_{n}=\delta_{0}(a b)^{n}-\delta_{1}\left(\frac{a}{b}\right)^{n}-\delta_{2}\left(\frac{a}{b^{3}}\right)^{n}-\cdots-\delta_{k}\left(\frac{a}{b^{2 k-1}}\right)^{n}$$
for sufficiently large $n$, where $\delta_{0}>0$ as $z_{n}>0$. As before, one obtains
$$\begin{aligned} & A^{n} B^{n}-A^{n}-B^{n}+1=z_{n}^{2} \\ & =\left\{\delta_{0}(a b)^{n}-\delta_{1}\left(\frac{a}{b}\right)^{n}-\delta_{2}\left(\frac{a}{b^{3}}\right)^{n}-\cdots-\delta_{k}\left(\frac{a}{b^{2 k-1}}\right)^{n}\right\}^{2} \\ & =\delta_{0}^{2} A^{n} B^{n}-2 \delta_{0} \delta_{1} A^{n}-\sum_{i=2}^{i=k}\left(2 \delta_{0} \delta_{i}-\sum_{j=1}^{j=i-1} \delta_{j} \delta_{i-j}\right)\left(\frac{A}{B^{i-1}}\right)^{n}+O\left(\frac{A}{B^{k}}\right)^{n} . \end{aligned}$$
Easy asymptotic calculations yield $\delta_{0}=1, \delta_{1}=\frac{1}{2}, \delta_{i}=\frac{1}{2} \sum_{j=1}^{j=i-1} \delta_{j} \delta_{i-j}$ for $i=2,3, \ldots, k-2$, and then $a=b^{k-1}$. It follows that $k>2$ and there is some $P \in \mathbb{Q}[X]$ for which $(X-1)(X^{k-1}-1)= P(X)^{2}$. But this cannot occur, for instance as $X^{k-1}-1$ has no double zeros. Thus our
assumption that $z_{n}$ was an integer for $n=1,2, \ldots$ turned out to be wrong, which solves the problem.