As has been mentioned many times on this site, anytime you have an algebraic function containing the square root of a cubic or a quartic, you are bound to bump into an elliptic integral.
Usually, such things are handled by using Jacobian elliptic functions for substitutions (in a manner similar to using substitution with trigonometric or hyperbolic functions when you have the square root of a quadratic in an integral).
I'll skip the tedious details of figuring out the proper substitution, since Byrd and Friedman give a formula for handling your integral (formula 212.00 in their handbook):
$$\int_y^\infty\frac{\mathrm dt}{\sqrt{(t^2+a^2)(t^2-b^2)}}=\frac1{\sqrt{a^2+b^2}}F\left(\arcsin\left(\sqrt{\frac{a^2+b^2}{a^2+y^2}}\right) \mid\frac{a^2}{a^2+b^2}\right)$$
where $F(\phi|m)$ is the incomplete elliptic integral of the first kind.
Coming back to your integral, we let $u=2\sqrt{3}-3$ and $v=2\sqrt{3}+3$ such that
$$\int_1^\infty\frac{\mathrm dx}{\sqrt{3x^4+6x^2-1}}=\frac1{\sqrt{uv}}\int_1^\infty \frac{\mathrm dx}{\sqrt{(x^2+1/u)(x^2-1/v)}}$$
Using the quoted formula, the integral reduces to
$$\frac1{\sqrt{u+v}}F\left(\arcsin\left(\sqrt{\frac{u+v}{v+uv}}\right)\mid\frac{v}{u+v}\right)$$
Substituting the values of $u$ and $v$ into this expression and simplifying, we have the result
$$\frac1{\sqrt[4]{48}}F\left(\arcsin\left(\sqrt{\sqrt{3}-1}\right)\mid\frac{2+\sqrt{3}}{4}\right)$$
which agrees with the numerical result in the comments.
As an aside, I consider it a capital annoyance that Mathematica often returns results with complex amplitudes even for real results...
In addition to the nice set of references by Raymond Manzoni, here is my proof of the identity. Frankly, I have not seen these references yet, thus I am not sure if this already appears in one of them.
Here I refer to the following identity
$$ \binom{\alpha}{\omega} = \frac{1}{2\pi} \int_{-\pi}^{\pi} \left(1 + e^{i\theta}\right)^{\alpha} e^{-i\omega \theta} \; d\theta, \ \cdots \ (1) $$
whose proof can be found in my blog post.
Now let $x$ be a real number such that $|x| < \frac{\pi}{2}$. Then simple calculation shows that
$$ \log\left(1+e^{2ix}\right) = \log(2\cos x) + ix \quad \Longleftrightarrow \quad \Im \left( \frac{-x}{\log\left(1+e^{2ix}\right)} \right) = \frac{x^2}{x^2 + \log^2(2\cos x)},$$
hence we have
$$ \begin{align*}I
&:= \int_{0}^{\frac{\pi}{2}} \frac{x^2}{x^2 + \log^2(2\cos x)} \; dx
= -\int_{0}^{\frac{\pi}{2}} \Im \left( \frac{x}{\log\left(1+e^{2ix}\right)} \right) \; dx \\
&= -\frac{1}{8}\int_{-\pi}^{\pi} \Im \left( \frac{\theta}{\log\left(1+e^{i\theta}\right)} \right) \; d\theta
= \frac{1}{8}\Re \left( \int_{-\pi}^{\pi} \frac{i\theta}{\log\left(1+e^{i\theta}\right)} \; d\theta \right).
\end{align*}$$
Differentiating both sides of $(1)$ with respect to $\omega$ and plugging $\omega = 1$, we have
$$ \frac{1}{2\pi} \int_{-\pi}^{\pi} (-i\theta) \left(1 + e^{i\theta}\right)^{\alpha} e^{-i\theta} \; d\theta = \alpha \left(\psi_0(\alpha) - \psi_0(2)\right). $$
Now integrating both sides with respect to $\alpha$ on $[0, 1]$,
$$ \begin{align*}
-\frac{1}{2\pi} \int_{-\pi}^{\pi} \frac{i\theta}{\log \left(1 + e^{i\theta}\right)} \; d\theta
&= \int_{0}^{1} \alpha \left(\psi_0(\alpha) - \psi_0(2)\right) \; d\alpha \\
&= \left[ \alpha \log \Gamma (\alpha) \right]_{0}^{1} - \int_{0}^{1} \log \Gamma (\alpha) \; d\alpha - \frac{1}{2}\psi_0(2) \\
&= -\frac{1}{2}\left( 1 - \gamma + \log (2\pi) \right),
\end{align*}$$
where we have used the fact that
$$ \psi_0 (1+n) = -\gamma + H_n, \quad n \in \mathbb{N}$$
and
$$ \begin{align*}
\int_{0}^{1} \log \Gamma (\alpha) \; d\alpha
& = \frac{1}{2} \int_{0}^{1} \log \left[ \Gamma (\alpha) \Gamma (1-\alpha) \right] \; d\alpha \\
&= \frac{1}{2} \int_{0}^{1} \log \left( \frac{\pi}{\sin \pi \alpha} \right) \; d\alpha \\
&= \frac{1}{2} \left( \log \pi - \int_{0}^{1} \log \sin \pi \alpha \; d\alpha \right) \\
&= \frac{1}{2} \log (2\pi).
\end{align*} $$
Therefore we have the desired result.
Best Answer
As I remarked in a comment, I have spent a while trying to derive this identity by means of contour integration.
Though I have not succeeded at that,(I take great pleasure in striking that out; see my other answer below!) I have found some related identities in my search which surprised me quite a bit, and I think they're worth sharing. For example, using the integral you asked about, I have shown that $$ \begin{equation} \sum_{n = 1}^\infty \frac{H_n}{n^3} = \frac{5}{4}\cdot \frac{\pi^4}{90} = \frac{5}{4} \zeta(4). \tag{1} \end{equation} $$ Here $$ H_n = \sum_{k = 1}^n \frac{1}{k} $$ is the $n$th harmonic number. Here are a couple more examples. I will use the first to derive $(1)$.More details about the derivations of 1. and 2. are contained in the answer I referenced above (near the bottom).
Let me now turn to deriving $(1)$. Taking $x = - e^{-2y}$ in the series expansion $$ \log{(1+x)} = \sum_{n = 1}^\infty \frac{(-1)^{n+1}x^n}{n} $$ and inserting the result into the integral evaluated in 1. gives $$ \begin{align} \frac{1}{8}\cdot\frac{\pi^4}{90} &= \int_0^\infty y \log^2(1 - e^{-2y})\,dy \\ & = \int_0^\infty y \left(\sum_{n = 1}^\infty \frac{e^{-2ny}}{n}\right)^2\,dy \\ & = \sum_{m = 1}^\infty\sum_{n = 1}^\infty\frac{1}{nm} \int_0^\infty y e^{-2(n + m)y}\,dy \\ & = \frac{1}{4}\sum_{m = 1}^\infty\sum_{n = 1}^\infty \frac{1}{nm(n+m)^2}. \tag{2} \end{align} $$ Now put $r = m+n$ and $s=n$ in order to write $$ \begin{align} \sum_{m = 1}^\infty\sum_{n = 1}^\infty \frac{1}{nm(n+m)^2} & = \sum_{r = 2}^\infty \frac{1}{r^2}\sum_{s = 1}^{r-1} \frac{1}{s(r-s)} = 2 \sum_{r=2}^\infty \frac{1}{r^3} \sum_{s=1}^{r-1}\frac{1}{s}, \tag{3} \end{align} $$ with the last equation a consequence of the identity $\frac{1}{s(r-s)} = \frac{1}{r}\left(\frac{1}{s}+ \frac{1}{r-s}\right)$. Insert $(3)$ into $(2)$ and multiply both sides by $2$ to get $$ \frac{1}{4} \cdot \frac{\pi^4}{90} = \sum_{r=2}^\infty \frac{1}{r^3} \sum_{s=1}^{r-1}\frac{1}{s} = \sum_{r=2}^\infty \frac{H_{r-1}}{r^3}. $$ Since $H_r - H_{r-1} = 1/r$, the identity $(1)$ now follows from the equations $$ \frac{\pi^4}{90} = \sum_{r=1}^\infty \frac{1}{r^4} = \sum_{r = 1}^\infty \frac{H_r}{r^3} - \sum_{r=2}^\infty \frac{H_{r-1}}{r^3} = \sum_{r = 1}^\infty \frac{H_r}{r^3} - \frac{1}{4}\cdot\frac{\pi^4}{90}. $$