Difference of two non-central chi squared random variables

I have a problem where I want to find the distribution of the difference of two non-central chi squared random variables (RV), both independent.
Given
$$
X=a+A\
Y=b+B
$$
where $a=a_r+ja_i$, $b=b_r+jb_i$, $a_r,a_i,b_r,b_iinmathbb{R}$ are constants, $A=A_r+jA_i$, $B=B_r+jB_i$, $A_r,A_i,B_r,B_isimmathcal{N}(0,1)$, all independent. The problem is calculating
$$
P(Z>z)=int_z^{infty}f_Z(x)dx
$$
where $Z=|X|^2-|Y|^2$. I tried the following approaches:

1. $|X|^2=X_r^2+X_i^2$, where $X_r$ and $X_i$ are the real and imaginary parts of $X$ and the same for $Y$. $X_r^2$ and $X_i^2$ are non-central chi squared with non-centrality parameter $lambda_X=|a|^2$ and $2$ degrees of freedom (DoF). I know that adding two non-central chi squared is also a non-central chi squared, but it is not clear to me what is the difference of two non-central chi squared.
2. Via Moment Generating Function (MGF). The MGF of the difference of two RV is the product of individual MGF, one with negative argument:
   $$
   M_{U-V}(t)=M_U(t)M_V(-t)
   $$
   for arbitrary $U$, $V$ RV. Thus,
   $$
   M_Z(t)=M_{|X|^2}(t)M_{|Y|^2}(-t)=frac{e^{frac{lambda_X t}{1-2t}-frac{lambda_Y t}{1+2t}}}{1-4t^2}
   $$
   The pdf can be found via the inverse of Laplace transform of the MGF, but I cannot see if the $M_Z$ has a closed-form inverse.
3. Convolving.
   $$
   f_Z(z)=int_{-infty}^{infty}f_{|X|^2}(x)f_{|Y|^2}(x+z)dx.
   $$
   I tried to solve this integral using:
   • The pdf expressed as the modified Bessel function of first kind $I_0$ :
     $$
     f_{|X|^2}(x)=frac{1}{2}e^{-frac{x+lambda_X}{2}}I_0left(sqrt{lambda_X x}right)
     $$
     In this case, the convolution is an integral of an exponential and product of two Bessel functions. I used the book of “Tables of Integrals” of Gradshteyn and Ryzhik but I could not find any solution.
   • The expression of Bessel $I_0$ as infinite series. At the end I have an integral of two infinite series. I tried to express this product as a Cauchy product but arguments are different (one is $x$ and the other is $x+z$).
   • The pdf expressed as a mixture. But at the end I have a convolution of two infinite series of functions.
4. Developing the squares. $Z$ can be rewritten as
   $$
   Z=|a|^2-|b|^2+2left(a_rA_r+a_iA_i-b_rB_r-b_iB_iright)+left(A_r^2+A_i^2-B_r^2-B_i^2right)=N+L
   $$
   where $Nsimmathcal{N}left(|a|^2-|b|^2,|a|^2+|b|^2right)$ and $Lsimmathcal{L}(0,2)$ is a Laplace RV. The sum of $N+L$, though, is not clear what distribution is, since $N$ and $L$ are dependent. The MGF of $N+L$ is the previous one, as expected.

I would appreciate any hint you could give me. Thank you.