Does this recurrence have a closed form limit $x_{n+1}=x_n-\frac{a}{3^{2n+1}x_n}$?
I would be astounded if the limit has a nice closed form. Here are some facts I have investigated.
Let $f(z) = 3z - \frac{1}{z}$. Then we easily check that $y_n = \frac{3^n}{\sqrt{a}} x_n$ satisfies $y_{n+1} = f(y_n)$. Thus $x_n$ has the following general form
$$ x_n = \frac{\sqrt{a}}{3^n} f^{\circ n}(y_0) = \sqrt{a} F_n\left(\frac{x_0}{\sqrt{a}}\right), \qquad F_n(z) := z \prod_{k=0}^{n-1} \left(1 - \frac{1}{3f^{\circ k}(z)^2} \right). $$
So it suffices to consider the case where $a = 1$ and then investigate the limit of $F_n(z)$ if exists. To this end, we claim the following:
Proposition. Let $I = [-\frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}}]$. Then the limit $F(z) = \lim_{n\to\infty} F_n(z)$ converges and defines a holomorphic function on $\Bbb{C}\setminus I$. Moreover,
There exists a finite Borel measure $\mu$ on $\Bbb{R}$ which is supported on $I$ and satisfies \begin{align*} F(z) &= z - \int_I \frac{\mu(d\lambda)}{z-\lambda} \\ &= z - \frac{3}{8z} - \frac{9}{640z^3} - \frac{2241}{465920z^5} - \cdots \quad \text{as } z \to \infty. \end{align*}
$F(f(z)) = 3F(z)$.
The proof of this fact is at the end of this answer. Here are some comments:
Identifying $F$ amounts to identifying the measure $\mu$. Although not entirely sure, it seems to me that $\mu$ is discrete and $\operatorname{supp}(\mu)$ is a Cantor-like set. That is, poles aggregate at uncountably many points and they do not average out. I can hardly imagine of an elementary function which exhibits this kind of behavior, so I suspect that $F$ is not an elementary function.
On $(\frac{1}{\sqrt{2}}, \infty)$ the function $f$ has inverse. If we denote by $g = f^{-1}$ the inverse, then for each $\epsilon > 0$ we can prove that
$$ g^{\circ n}\left(\frac{1}{\sqrt{2}} + \epsilon\right) = \frac{1}{\sqrt{2}} + \frac{\ell + o(1)}{5^n} \qquad \text{as } n\to\infty \tag{*}$$
for some constant $\ell = \ell(\epsilon) \geq$ which depends on $\epsilon$. Then in view of the identity
$$ F\left(\frac{1}{\sqrt{2}} + \epsilon\right) = 3^n F\left(g^{\circ n}\left(\frac{1}{\sqrt{2}} + \epsilon\right) \right), $$
if $F(\frac{1}{\sqrt{2}} + \epsilon) \sim c \epsilon^{\alpha}$ for some constants $c > 0$ and $\alpha$, then we must have $\alpha = \frac{\log 3}{\log 5}$. Although I have not formally written down my idea, using a detailed version of $\text{(*)}$ I checked that $F(\frac{1}{\sqrt{2}} + \epsilon) = \Theta(\epsilon^{\alpha})$ for this $\alpha$, where $\Theta$ is a Landau asymptotic notation.
Proof of Proposition. By the extra assumption $a =1$, we have $x_n = F_n(x_0)$ and $y_n = f^{\circ n}(x_0)$. From this, we have the following relations.
$$ F_{n+1}(z) = F_n(z) - \frac{1}{3^{2n+1}F_n(z)}, \qquad F_n(z) = 3^{-n} f^{\circ n}(z). \tag{1}$$
Using the first relation in $\text{(1)}$, we have
$$ \operatorname{Im}(F_{n+1}(z)) = \operatorname{Im}(F_n(z))\left( 1 + \frac{1}{3^{2n+1}|F_n(z)|^2} \right). $$
In particular, $F_n$ is a Nevanlinna function and thus we can write write $F_n(z)$ as
$$F_n(z) = z - \int_{\Bbb{R}} \frac{\mu_n(d\lambda)}{z-\lambda} $$
for some point mass $\mu_n$ on $\Bbb{R}$. (This can also be proved directly by some tedious algebra.) This representation shows that $F_n(z) = z - \mu_n(\Bbb{R})z^{-1} + \mathcal{O}(z^{-2})$ as $z \to \infty$. Plugging this to the recurrence relation in $\text{(1)}$, we find that
$$\mu_{n+1}(\Bbb{R}) = \mu_n(\Bbb{R}) + \frac{1}{3^{2n+1}}. $$
In particular, the total mass of $(\mu_n)$ is bounded. On the other hand, if $x \in \Bbb{R}\setminus I$, then we easily check that $|f(x)| > |x|$. This implies that
$\text{(2)}$ $F_n$ is finite on $\Bbb{R}\setminus I$, and thus $\mu_n$ is supported on $I$.
$\text{(3)}$ $F_n$ converges on $\Bbb{R}\setminus I$, since $F_n(x)$ is monotone and bounded for each $x \in \Bbb{R}\setminus I$.
So it follows from $\text{(2)}$ that $(\mu_n)$ is weakly compact and hence $F_n$ is a normal family on $\Bbb{C}\setminus I$. Then the first claim follows from the identity theorem together with $\text{(3)}$.
Once we know that $F_n$ converges on $\Bbb{C}\setminus I$, then the second claim follows by taking the limit to $F_n(f(z)) = 3F_{n+1}(z)$. Using this, we can compute the Laurent expansion of $F$ near $\infty$. Since $F$ is an odd function, with $a_n = \int_{\Bbb{R}} \lambda^{2n} \, \mu(d\lambda)$ we have
$$ F(z) = z - \sum_{n=0} \frac{a_n}{z^{2n+1}}. $$
Plugging this to the identity $F(f(z)) = 3F(z)$ produces an system of equations for $(a_n)$, from which we can compute $(a_n)$ at least theoretically.