$\newcommand{\ep}{\varepsilon}$
Each of the two proposed inequalities can be written in the form $H(a)>0$ for $a\in(1,\infty)$, where $H$ is a smooth function. All such inequalities -- but only for $a$ in a finite closed interval $I$ -- can be proved, at least in principle, by the interval method; that is, by partitioning $I$ into many small subintervals and bounding, tightly enough, terms in the expression of $H(a)$ on each small subinterval. A difficulty of using the interval method here is that $x_1$ and $x_2$ as functions of $a$ are known only implicitly. Another difficulty is that the interval $I$ is neither finite nor closed, and so, we shall need to consider separately ("small") values of $a$ close to $1+$ and large values of $a$, close enough to $\infty$. For intermediate values of $a$, we shall indeed use the interval method.
Therefore, the solution presented below is long and technical. Since the problem does not seem to belong to a general class of problems, I'd be surprised if there is an "official" solution to it.
All and, in particular, large values of $a$. Let
\begin{equation*}
g(x):=x^2-x-\ln x
\end{equation*}
for $x>0$. Then $g(1)=g'(1)=0$ and $g$ is strictly convex, with $g(0+)=g(\infty-)=\infty$. So, $g$ is strictly decreasing from $\infty$ to $0$ on $(0,1]$ and strictly increasing from $0$ to $\infty$ on $[1,\infty)$. So, for each real $a>0$ there are uniquely determined $x_1=x_1(a)\in(0,1)$ and $x_2=x_2(a)\in(1,\infty)$ such that
\begin{equation*}
g(x_1)=g(x_2)=\ln a;
\end{equation*}
of course, these $x_1$ and $x_2$ are the same as in the OP.
One may note that, if $a\to\infty$, then the conditions $g(x_1)=\ln a$ and $x_1\in(0,1)$ imply, successively, the following: $\ln x_1\to-\infty$, $x_1\to0$, $-\ln x_1=\ln a+o(1)$, and $x_1\sim1/a$.
Also, again if $a\to\infty$, then the conditions $g(x_2)=\ln a$ and $x_2>1$ imply $x_2^2\sim\ln a$ and hence $x_2\sim\sqrt{\ln a}$. So, if $a\to\infty$, then $x_1x_2\sim\frac1a\,\sqrt{\ln a}$, while the proposed lower and upper bounds on $x_1x_2$ are $\frac3{2a+1}\sim\frac3{2a}$ and $\frac{\ln a}{a-1}\sim\frac{\ln a}a$.
By tweaking the asymptotic expressions $1/a$ and $\sqrt{\ln a}$ for $x_1$ and $x_2$ as $a\to\infty$, we obtain lower and upper bounds on $x_1$ and $x_2$, useful for large enough $a$, as follows.
The conditions $g(x_1)=\ln a$ and $x_1\in(0,1)$ yield $x_1^2-x_1<0$ and hence $\ln a=g(x_1)<-\ln x_1$, so that we get an upper bound on $x_1$:
\begin{equation*}
x_1<X_1^{up}:=X_1^{up}(a):=\tfrac1a.
\end{equation*}
Therefore and because $x^2-x$ is decreasing in $x\in(0,1/2]$, we have $x_1^2-x_1>\frac1{a^2}-\frac1a\ge\frac1{2^2}-\frac12=-\frac14$ if $a\ge2$, whence $\ln a=g(x_1)>-\frac14-\ln x_1$, so that we get a lower bound on $x_1$:
\begin{equation*}
x_1>X_1^{lo}:=X_1^{lo}(a):=\tfrac1{e^{1/4}a}\quad\text{for}\quad a\ge2.
\end{equation*}
The conditions $g(x_2)=\ln a$ and $x_2>1$ yield $x_2^2=x_2+\ln x_2+\ln a\ge 1+\ln a$, so that we get a lower bound on $x_2$:
\begin{equation*}
x_2>X_2^{lo}:=X_2^{lo}(a):=\sqrt{1+\ln a}.
\end{equation*}
Therefore and because $h(x):=\frac{x+\ln x}{x^2}$ is decreasing in $x\ge1$ (from $1$ to $0$), we have $h(x_2)<h(\sqrt{1+\ln a})\le h(\sqrt{1+\ln 17})$ if $a\ge17$. So, $x_2^2=h(x_2)x_2^2+\ln a<h(\sqrt{1+\ln 17})x_2^2+\ln a$, and
we get an upper bound on $x_2$:
\begin{equation*}
x_2<X_2^{up}:=X_2^{up}(a):=\frac{\sqrt{\ln a}}c \quad\text{for}\quad a\ge17,
\end{equation*}
where $c:=\sqrt{1-h(\sqrt{1+\ln 17})}\approx0.560$.
Consider now the ratio
\begin{equation*}
R^{up}(a):=X_1^{up}X_2^{up}\Big/\frac{\ln a}{a-1}=\frac{a-1}{ca\sqrt{\ln a}}.
\end{equation*}
It is easy to see that $R^{up}(a)$ decreases in $a>4$, and $R^{up}(17)\approx0.9979<1$. So,
\begin{equation*}
x_1x_2<X_1^{up}X_2^{up}<\frac{\ln a}{a-1} \quad\text{for}\quad a\ge17.
\end{equation*}
Similarly, the ratio
\begin{equation*}
R^{lo}(a):=X_1^{lo}X_2^{lo}\Big/\frac3{2a+1}=\frac{(2a+1)\sqrt{1+\ln a}}{3e^{1/4}a}
\end{equation*}
increases in $a>1$, and $R^{lo}(11)\approx1.00056>1$. So,
\begin{equation*}
x_1x_2>X_1^{lo}X_2^{lo}>\frac3{2a+1} \quad\text{for}\quad a\ge11.
\end{equation*}
Small values of $a$. Here we consider somewhat small values of
\begin{equation*}
\ep:=a-1>0.
\end{equation*}
Then Maclaurin expansions of the roots in powers of $\sqrt\ep$ are as follows:
\begin{equation*}
x_i\approx1+\sum_{j=1}^4c_j\,(-\sqrt\ep)^{ij},
\end{equation*}
where $i\in\{1,2\}$ and
\begin{equation*}
(c_1,\dots,c_4):=\Big(\sqrt{\frac{2}{3}},\frac{2}{27},-\frac{277}{486 \sqrt{6}},-\frac{893}{32805}\Big).
\end{equation*}
One does not have to check this, as we are just going to tweak these expansions a bit to obtain true lower and upper bounds on $x_1$ and $x_2$ for small enough $\ep=a-1>0$.
Indeed, let
\begin{align}
x_1^{lo}&:=x_1^{lo}(\ep):=1+\sum_{j=1}^3c_j\,(-\sqrt\ep)^{j}-\frac8{100}\ep^2, \\
x_2^{lo}&:=x_2^{lo}(\ep):=1+\sum_{j=1}^4c_j\,(\sqrt\ep)^{j}, \\
x_1^{up}&:=x_1^{up}(\ep):=1+\sum_{j=1}^4c_j\,(-\sqrt\ep)^{j}, \\
x_2^{up}&:=x_2^{up}(\ep):=1+\sum_{j=1}^3c_j\,(\sqrt\ep)^{j}+\frac1{20}\ep^2.
\end{align}
The derivatives of $d_i^{lo}(\ep):=g(x_i^{lo})-\ln a=g(x_i^{lo})-\ln(1+\ep)$ and $d_i^{up}(\ep):=g(x_i^{up})-\ln a=g(x_i^{up})-\ln(1+\ep)$ are rational functions of $\sqrt\ep$, and monotonicity patterns of these functions can be quickly determined algorithmically, using a computer algebra program. E.g., $(d_1^{lo})'(\ep)>0$ for $\ep\in(0,\frac{175}{1000}]$ and hence $d_1^{lo}(\ep)$ increases in $\ep\in(0,\frac{175}{1000}]$, from $d_1^{lo}(0)=0$. So, on $(0,\frac{175}{1000}]$ we have $d_1^{lo}>0$, that is, $g(x_1^{lo})>\ln a=g(x_1)$; since $g$ is decreasing on $(0,1]$, it follows that
\begin{equation*}
x_1>x_1^{lo}\quad\text{for}\quad \ep\in(0,\tfrac{175}{1000}].
\end{equation*}
Similarly we get
\begin{align}
x_2&>x_2^{lo}\quad\text{for}\quad \ep\in(0,\tfrac{175}{1000}], \\
x_1&<x_1^{up}\quad\text{for}\quad \ep\in(0,\tfrac{127}{1000}], \\
x_2&<x_2^{up}\quad\text{for}\quad \ep\in(0,\tfrac{127}{1000}].
\end{align}
Further, $x_1^{lo}x_2^{lo}-\frac3{2a+1}$ and the derivative of $x_1^{lo}x_2^{lo}-\frac{\ln a}{a-1}$ are rational functions of $\sqrt\ep$, which allows one to similarly use a computer algebra program to quickly show that
\begin{equation*}
x_1x_2<x_1^{up}x_2^{up}<\frac{\ln a}{a-1} \quad\text{for}\quad a\in(1,1+\tfrac{175}{1000}],
\end{equation*}
\begin{equation*}
x_1x_2>x_1^{lo}x_2^{lo}>\frac3{2a+1} \quad\text{for}\quad a\in(1,1+\tfrac{127}{1000}].
\end{equation*}
It remains to consider
Intermediate values of a. In view of the above considerations of large and small values of $a$, it remains to prove the lower bound $\frac3{2a+1}$ on $x_1x_2$ for $a\in[1+\tfrac{175}{1000},11]$ and the lower bound $\frac{\ln a}{a-1}$ on $x_1x_2$ for $a\in[1+\tfrac{127}{1000},17]$.
This is done by the mentioned interval method.
Note that the monotonicity pattern of the function $g$ allows one to approximate the values of $x_1$ and $x_2$ with any prescribed degree of accuracy, for each $a>1$. However, here it is slightly more convenient to use the explicit (albeit very complicated) equivalent form of the inequalities in question, given in the last display in the OP:
\begin{equation*}
\frac3{2e^{g(F(t))}+1}< t F(t)^2<\frac{g(F(t))}{e^{g(F(t))}-1}, \tag{1}
\end{equation*}
where
\begin{equation*}
F(t):=\frac{t-1+\sqrt{(t-1)^2+4\ln{t}\cdot(t^2-1)}}{2(t^2-1)}=x_1\in(0,1)
\end{equation*}
and $t:=t(a):=x_2/x_1>1$. (The OP used the symbol $f$ to denote both $F$ and $g(x)-\ln a$.)
Note that $x_1$ is decreasing and $x_2$ is increasing in $a$, so that $t$ is increasing in $a$. So, for $a\in[1+\tfrac{127}{1000},17]$ we have $t\in[t(1+\tfrac{127}{1000}),t(17)]\subseteq[t^{small},t^{large}]$,
where $t^{small}:=x_2^{up}(\tfrac{127}{1000})/x_1^{lo}(\tfrac{127}{1000})\approx1.7699$ and $t^{large}:=X_2^{up}(17)/X_1^{lo}(17)=53.974\ldots$.
Since $t$ is increasing and $x_1=F(t)$ is decreasing in $a$, it follows that $F(t)\in(0,1)$ is decreasing in $t>1$. Therefore and because $g$ is decreasing on $(0,1)$, $g(F(t))$ is increasing in $t$. Also, $\frac b{e^b-1}$ is decreasing in $b>0$. So,
\begin{equation*}
\frac{g(F(t))}{e^{g(F(t))}-1}
\end{equation*}
is decreasing in $t>1$. So, to prove the second inequality in (1) for $a\in[1+\tfrac{127}{1000},17]$, it suffices to show that
\begin{equation*}
t_{j+1} F(t_j)^2<\frac{g(F(t_{j+1}))}{e^{g(F(t_{j+1}))}-1} \tag{2}
\end{equation*}
for all integers $j\in\{0,\dots,53n\}$, where $n=300$,
$t_j:=\frac{176}{100}+\frac jn$, so that $t_{53n}=53+\frac{176}{100}\ge53.974\ldots$.
Indeed, then it will follow that
\begin{equation*}
t F(t)^2\le t_{j+1} F(t_j)^2<\frac{g(F(t_{j+1}))}{e^{g(F(t_{j+1}))}-1}
\le\frac{g(F(t))}{e^{g(F(t))}-1}
\end{equation*}
for all $j\in\{0,\dots,53n\}$ and all $t\in[\frac{176}{100}+\frac jn,\frac{176}{100}+\frac{j+1}n]$, so that we will indeed have the second inequality in (1) for all $t\in[t(1+\tfrac{127}{1000}),t(17)]$.
The proof of (2) is done by straightforward calculation (taking about 14 sec with Mathematica).
The proof of the first inequality in (1) (for $t\in[t(1+\tfrac{175}{1000}),t(11)]$) is quite similar (and a bit easier); it takes about 6 sec with Mathematica.
The proof is now complete.
