Let $U=(u_1, u_2)$ and $V=(v_1, v_2)$ be two randomly distributed points on the Euclidean plane assuming bivariate normal distributions $U \sim N(\mu_u, \Sigma_u)$ and $V \sim N(\mu_v, \Sigma_v)$ with unequal means $\mu_u = (\mu_{u1}, \mu_{u2})$, $\mu_v =( \mu_{v1}, \mu_{v2} )$ and unequal covariance matrices
$$ \Sigma_u = \begin{bmatrix} \sigma_{u1}^2 & 0 \\ 0 & \sigma_{u2}^2 \end{bmatrix} % , \Sigma_v= \begin{bmatrix} \sigma_{v1}^2 & 0 \\ 0 & \sigma_{v2}^2 \end{bmatrix}. $$
What is the expected value $\mathbb{E}(r)$ of the Euclidean distance $r=\sqrt{(u_1-v_1)^2+(u_2-v_2)^2}$ between $U$ and $V$? Is there a closed-form expression for the moment generating function of $r$? Note that I am specifically seeking moments of $r$ and not $r^2$.
This question is equivalent to deriving the moments of the Beckmann distribution as follows. Since $U$ and $V$ are both uncorrelated and independently distributed bivariate normal distributions, it follows that $X=U_1-V_1$ and $Y=U_2-V_2$ are also independently distributed normal variables with means $\mu_x = \mu_{u1}-\mu_{v1}$, $\mu_y = \mu_{u2}-\mu_{v2}$ and variances $\sigma_x^2 = \sigma_{u1}^2+\sigma_{v1}^2$, $\sigma_y^2 = \sigma_{u2}^2+\sigma_{v2}^2$. Then $r=\sqrt{X^2+Y^2}$ can be described by a Beckmann distribution [1]:
$$ f(r) = \frac{r}{2\pi\sigma_x\sigma_y} \int_0^{2\pi} \exp{ \left( -\frac{(r \cos(\theta)-\mu_x)^2}{2\sigma_x^2} - \frac{(r \sin(\theta)-\mu_y)^2}{2\sigma_y^2} \right)}d\theta, \quad r\geq 0, $$ where $X$ and $Y$ are transformed into polar coordinates using $X=r \cos \theta$ and $Y=r \sin \theta$. Its special cases comprise the Hoyt ($\mu_x=\mu_y=0$), Rice ($\sigma_x=\sigma_y$) and Rayleigh ($\mu_x=\mu_y=0$, $\sigma_x=\sigma_y$) distributions (also see [2, 3]). Unlike some of the pdfs of these special cases, the Beckmann distribution does not have a closed-form expression.
Tried so far
The moment generating function of the squared Beckmann-distributed variable $r^2$ is widely used in communication theory and given by [4]: $$ M_{r^2}(t) = \frac{1}{\sqrt{1-2\sigma_x^2t}\sqrt{1-2\sigma_y^2t}} \exp \left( \frac{\mu_x^2t}{1-2\sigma_x^2t} + \frac{\mu_y^2t}{1-2\sigma_y^2t} \right). $$ Its closed-form is a result of the property that $r^{2n}=(x^2+y^2 )^n$ can be rewritten as $r^{2n}=\sum_{k=0}^n \binom{n}{k} X^{2n} Y^{2(n-k)}$. Therefore, the integral expressions in $M_{r^2}(t)=\mathbb{E}[r^{2n} e^{tr^2} ]$ can be separated into distinct terms that either rely on $X$ or on $Y$, but not on both [5]. Since the 'unsquared' variable $r=\sqrt{x^2+y^2 }$ cannot be separated in the same way, I soon got stuck while trying to compute $M_{r}(t)=\mathbb{E}[r^{n} e^{tr} ]$ directly using the same approach.
An alternative approach could be to transform the result for $M_{r^2}(t)$ above directly into $M_{r}(t)$ using the same steps as in this related question. Unfortunately I have not yet been able to rewrite the result into a closed-form for the expected value of $r$.
