Define the subset $D\subseteq\mathbb{R}_{>0}^{2}$ of the first quadrant by the inequalities,
$$D:=\{\left(x,y\right)\in\mathbb{R}_{>0}^{2}\mid1\le xy\le 3\land1\le x^{2}-y^{2}\le3\}.$$
Given a density function $\sigma:\mathbb{R}^{2}\rightarrow\mathbb{R}_{\ge0}$ defined via the polynomial
$$\sigma{\left(x,y\right)}:=\left(x^{2}+y^{2}\right)\left(x^{4}+y^{4}\right),$$
the total mass of the region $D$ will be equal to the double integral
$$\operatorname{m}{\left(D\right)}=\iint_{D}\mathrm{d}x\,\mathrm{d}y\,\sigma{\left(x,y\right)}.$$
As you discovered for yourself, changing to polar coordinates doesn't seem to make the region of integration any less complicated. If we want the new integration region to be as simple as possible after changing variables, then you need a substitution that maps the region $D$ to a rectangle $R=\left[a,b\right]\times\left[c,d\right]\subseteq\mathbb{R}^{2}$, for some $a,b,c,d\in\mathbb{R}\land a<b\land c<d$.
Consider the following change of variables:
$$\begin{cases}
&xy=u\\
&\frac{x^{2}-y^{2}}{2}=v.\\
\end{cases}$$
For each ordered pair $\left(u,v\right)\in\mathbb{R}_{>0}^{2}$, there is a (unique) ordered pair $\left(x,y\right)\in\{\left(a,b\right)\in\mathbb{R}_{>0}^{2}\mid0<b<a\}$ satisfying the change-of-variable relations above, i.e., the transformation is a bijection from $\{\left(a,b\right)\in\mathbb{R}_{>0}^{2}\mid0<b<a\}$ to $\mathbb{R}_{>0}^{2}$ with inverse substitution
$$\begin{cases}
&x=\sqrt{v+\sqrt{u^{2}+v^{2}}}\\
&y=\frac{u}{\sqrt{v+\sqrt{u^{2}+v^{2}}}}.\\
\end{cases}$$
The determinant of the Jacobian matrix is found to be (after some perhaps tedious calculations),
$$\begin{align}
\frac{\partial\left(x,y\right)}{\partial\left(u,v\right)}
&=\det\pmatrix{\frac{u}{2\sqrt{u^{2}+v^{2}}\sqrt{v+\sqrt{u^{2}+v^{2}}}}
&\frac{\sqrt{v+\sqrt{u^{2}+v^{2}}}}{2\sqrt{u^{2}+v^{2}}}\\\
\frac{\sqrt{v+\sqrt{u^{2}+v^{2}}}}{2\sqrt{u^{2}+v^{2}}}
&-\frac{u}{2\sqrt{u^{2}+v^{2}}\sqrt{v+\sqrt{u^{2}+v^{2}}}}}\\
&=-\frac{1}{2\sqrt{u^{2}+v^{2}}}.\\
\end{align}$$
Finally, our change of variables yields the following result:
$$\begin{align}
\operatorname{m}{\left(D\right)}
&=\iint_{D}\mathrm{d}x\,\mathrm{d}y\,\sigma{\left(x,y\right)}\\
&=\iint_{R}\mathrm{d}u\,\mathrm{d}v\,\left|\frac{\partial\left(x,y\right)}{\partial\left(u,v\right)}\right|\sigma{\left(x\left(u,v\right),y\left(u,v\right)\right)}\\
&=\iint_{\left[1,3\right]\times\left[\frac12,\frac32\right]}\mathrm{d}u\,\mathrm{d}v\,\frac{2\left(u^{2}+2v^{2}\right)\cdot2\sqrt{u^{2}+v^{2}}}{2\sqrt{u^{2}+v^{2}}}\\
&=2\int_{1}^{3}\mathrm{d}u\int_{\frac12}^{\frac32}\mathrm{d}v\,\left(u^{2}+2v^{2}\right)\\
&=26.\blacksquare\\
\end{align}$$
The simplicity of the last integral is our payoff for all the tedious intermediate algebraic steps (many of which I've omitted here for the sake of space).
Ending Note: The coordinate transformation used in the solution above actually does have a name. In addition to Cartesian and polar coordinates, there are also parabolic coordinates, and if you're really curious there is an entire zoo of exotic coordinate systems you can read about there as well, but I digress... ;)