Recurrence relation involving derivative

Question

I have the following recurrence equation: \begin{align*} Z_{n+1}(x, y) = -2i\frac{\partial Z_n}{\partial x}(x, y) - Z_{n-1}(x, y), \end{align*} with the following conditions: \begin{equation} \begin{cases} Z_0(x, y) = 2\pi J_0\left(\sqrt{x^2 + y^2}\right)e^{id}, \\ Z_1(x, y) = \frac{1}{i}\frac{\partial Z_0}{\partial x}(x, y) + \frac{\partial Z_0}{\partial y}(x, y), \end{cases} \end{equation} where $J_0$ is the Bessel function of the first kind of order $0$ and $d$ is a real number.

I would like to get a closed form of $Z_n$ involving Bessel functions, but I really don't know if there is a way to do it.

Remark: This question is a follow-up of another one, for more context please refer to this post where you can see the development to obtain the recurrence equation.

Edit: following @Sal advice

I am not the best user of Fourier transform so bare with me if there is any big mistakes. Defining the fourier transform of $Z_n$ as: \begin{equation*} \tilde{Z}_n(u, v) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} Z_n(x,y) e^{-i (xu+yv)} dx dy \end{equation*} and taking the Fourier transform of my recurrence relation, holds: \begin{equation*} \tilde{Z}_{n+1}(u, v) = 2u\tilde{Z}_n(u, v) - \tilde{Z}_{n-1}(u, v). \end{equation*}

Using the following characteristic equation: \begin{equation*} \lambda^2 - 2u\lambda + 1 = 0, \end{equation*} we can compute two solutions: \begin{equation} \begin{cases} \lambda_1 = u - \sqrt{u^2 - 1}, \\ \lambda_2 = u + \sqrt{u^2 - 1}. \end{cases} \end{equation} Using all that, we have: \begin{equation*} \tilde{Z}_n(u, v) = a\lambda_1^n + b\lambda_2^n, \end{equation*} where the coefficients $a$ and $b$ are given by the conditions: \begin{equation} \begin{cases} \tilde{Z}_0 = a+b, \\ \tilde{Z}_1 = a\lambda_1 + b\lambda_2. \end{cases} \Rightarrow \begin{cases} a=\frac{\lambda_2\tilde{Z}_0 - \tilde{Z}_1}{\lambda_2-\lambda_1}, \\ b=\frac{\tilde{Z}_1 - \lambda_1\tilde{Z}_0}{\lambda_2-\lambda_1}. \end{cases} \end{equation} Finally, we can write that: \begin{align*} \tilde{Z}_n(u, v) &= \frac{1}{2\sqrt{u^2 - 1}}\left[(\lambda_2\tilde{Z}_0 - \tilde{Z}_1)\lambda_1^n + (\tilde{Z}_1 - \lambda_1\tilde{Z}_0)\lambda_2^n\right] \end{align*}

Edit 2:

We can do some more simplifications by remarking that: \begin{equation*} \tilde{Z}_1 = u\tilde{Z}_0 + iv\tilde{Z}_0, \end{equation*} such that: \begin{equation} \begin{cases} \lambda_2\tilde{Z}_0 - \tilde{Z}_1 = \left(\sqrt{u^2-1} - iv\right)\tilde{Z}_0 \\ \tilde{Z}_1 - \lambda_1\tilde{Z}_0 = \left(iv + \sqrt{u^2-1}\right)\tilde{Z}_0 \end{cases} \end{equation}

Finally, we can write: \begin{equation*} \tilde{Z}_n(u, v) = \tilde{Z}_0\frac{\lambda_1^n + \lambda_2^n}{2} - \tilde{Z}_0\frac{iv}{\sqrt{u^2 - 1}}\frac{\lambda_1^n - \lambda_2^n}{2} \end{equation*} which, by using the Binomial theorem, can be put into the form: \begin{equation*} \tilde{Z}_n(u, v) = \tilde{Z}_0\sum_{k=0}^{n} \begin{pmatrix} n\\ k \end{pmatrix} u^{n-k}\left(u^2 - 1\right)^\frac{k}{2}\left(\frac{1+(-1)^k}{2}+\frac{1-(-1)^k}{2}\frac{iv}{\sqrt{u^2-1}}\right) \end{equation*}

Edit 3:

Following the computation done in this post, we can find that: \begin{equation} \tilde{Z}_0(u, v) = 2\pi e^{id} \delta(\sqrt{u^2 + v^2}-1), \\ \end{equation}

Now, my problem is to compute the inverse Fourier transform to get $Z_n$.

Answer 1

Doing some brute force computation using Sympy, we get the following for the first few terms (the term $2\pi e^{id}$ is omitted): \begin{align*} Z_0(x, y) &= J_0(r)\\ Z_1(x, y) &= \left(i\frac{x+iy}{x^2+y^2}\right)rJ_1(r)\\ Z_2(x, y) &= \left(\frac{x+iy}{x^2 + y^2}\right)^2 r^2J_0(r) + 2\left(i\frac{x+iy}{x^2+y^2}\right)^2 rJ_1(r)\\ Z_3(x, y) &= 4i\left(\frac{x+iy}{x^2 + y^2}\right)^3 r^2J_0(r) + \left(i\frac{x+iy}{x^2+y^2}\right)^3\left[8 - \left(x^2 + y^2\right)\right] rJ_1(r)\\ Z_4(x, y) &= \left(\frac{x+iy}{x^2 + y^2}\right)^4\left[-24 + \left(x^2+y^2\right)\right] r^2J_0(r) + \left(i\frac{x+iy}{x^2+y^2}\right)^4\left[48 - 8\left(x^2 + y^2\right)\right] rJ_1(r)\\ Z_5(x, y) &= \left(\frac{x+iy}{x^2 + y^2}\right)^5\left[-192i + 12i\left(x^2+y^2\right)\right] r^2J_0(r) + \left(i\frac{x+iy}{x^2+y^2}\right)^5\left[384 - 72\left(x^2 + y^2\right) + \left(x^2 + y^2\right)^2\right]rJ_1(r)\\ Z_6(x, y) &= \left(\frac{x+iy}{x^2 + y^2}\right)^6\left[1920 - 144\left(x^2+y^2\right) + \left(x^2+y^2\right)^2\right] r^2J_0(r) + \left(i\frac{x+iy}{x^2+y^2}\right)^6\left[3840 - 768\left(x^2 + y^2\right) + 18\left(x^2 + y^2\right)^2\right]rJ_1(r)\\ Z_7(x, y) &= \left(\frac{x+iy}{x^2 + y^2}\right)^7\left[23040i-1920i\left(x^2+y^2\right) + 24i\left(x^2+y^2\right)^2\right] r^2J_0(r) + \left(i\frac{x+iy}{x^2+y^2}\right)^7\left[46080 - 9600\left(x^2 + y^2\right) + 288\left(x^2 + y^2\right)^2 - \left(x^2 + y^2\right)^3\right]rJ_1(r)\\ Z_8(x, y) &= ... \end{align*}

where $r=\sqrt{x^2 + y^2}$. Clearly, we can see a pattern.

Looking at the brute force computation, we can say that: \begin{align*} \forall n\geq 0, \quad Z_{n} = \left(\frac{x+iy}{x^2+y^2}\right)^{n}\left[r^2 J_0(r) \sum_{k=1}^{u_n}a_k^{n} \left(x^2 + y^2\right)^{k-1} + i^n r J_1(r) \sum_{k=1}^{u_{n+1}}b_k^{n} \left(x^2 + y^2\right)^{k-1}\right] \end{align*} where $u_n$ is defined as $u_n = \left(2n+(-1)^n -1\right)/4$ and $n\geq 1$. Using the brute force computations, we can write the first few terms of $a_k^n$ and $b_k^n$ \begin{array}{|c|c|c|c|c|c|c|c|} \hline a_1^2 = 1 & a_1^3 = 4i & a_1^4=-24 & a_1^5=-192i & a_1^6=1920 & a_1^7=23040i & a_1^8=-322560 & a_1^9=-5160960i\\ \hline && a_2^4 = 1 & a_2^5 = 12i & a_2^6=-144 & a_2^7=-1920i & a_2^8=28800& a_2^9=483840i\\ \hline &&&& a_3^6 = 1 & a_3^7=24i & a_3^8=-480& a_3^9=-9600i\\ \hline &&&&&& a_4^8=1& a_4^9=40i\\ \hline \end{array} and \begin{array}{|c|c|c|c|c|c|c|c|c|} \hline b_1^1 = 1 & b_1^2 = 2 & b_1^3=8 & b_1^4=48 & b_1^5=384 & b_1^6=3840 & b_1^7=46080 & b_1^8=645120 & b_1^9=10321920\\ \hline && b_2^3 = -1 & b_2^4 = -8 & b_2^5 = -72 & b_2^6=-768 & b_2^7=-9600 & b_2^8=-138240& b_2^9=-2257920\\ \hline &&&& b_3^5 = 1 & b_3^6=18 & b_3^7=288 & b_3^8=4800& b_3^9=86400\\ \hline &&&&&& b_4^7=-1 & b_4^8=-32& b_4^9=-800\\ \hline &&&&&&&& b_5^9=1\\ \hline \end{array} from which we can compute relations for the two first coefficients: \begin{align*} \begin{pmatrix} b_1^{n+1} = 2n b_1^n\\ b_1^1= 1 \end{pmatrix}&\Rightarrow b_1^{n} = 2^{n-1}(n-1)! \quad \forall n > 0 \\ \begin{pmatrix} a_1^{n+1} = 2in a_1^n\\ a_1^2= 1 \end{pmatrix}&\Rightarrow a_1^{n} = (2i)^{n-2}(n-1)! \quad \forall n > 1 \end{align*}

Doing the same for higher order coefficients, I obtained: \begin{align*} b_1^{n}&=2^{n-1}(n-1)! \quad \forall n>0\\ b_2^{n}&=-2^{n-3}(n-3)!(n-2)^2=-2^{n-3}(n-2)!(n-2) \quad \forall n>2\\ b_3^{n}&=2^{n-5}(n-5)!\left(\frac{(n-4)(n-3)}{2}\right)^2=2^{n-5}(n-3)!\frac{(n-4)(n-3)}{2^2} \quad \forall n>4\\ b_4^{n}&=-2^{n-7}(n-7)!\left(\frac{(n-6)(n-5)(n-4)}{6}\right)^2=-2^{n-7}(n-4)!\frac{(n-6)(n-5)(n-4)}{6^2}\quad \forall n>6 \end{align*} which can be put into the general form of: \begin{equation*} \forall k\geq 1, n\geq 2k-1,\quad b_k^{n} = (-1)^{k+1}2^{n-(2k-1)}\frac{(n-k)!(n-k)!}{(k-1)!(k-1)!(n-(2k-1))!)} \end{equation*}

Doing the same thing for the other coefficients, I obtained: \begin{align*} a_1^{n} &= (2i)^{n-2}(n-1)! \quad \forall n > 1\\ a_2^{n} &= (2i)^{n-4}(n-3)!\frac{(n-3)(n-2)}{2} = (2i)^{n-4}(n-2)!\frac{(n-3)}{2} \quad \forall n>3\\ a_3^{n} &= (2i)^{n-6}(n-5)!\frac{(n-5)(n-4)^2(n-3)}{12} = (2i)^{n-6}(n-3)!\frac{(n-5)(n-4)}{12} \quad \forall n>5\\ a_4^{n} &= (2i)^{n-8}(n-4)!\frac{(n-5)(n-6)(n-7)}{144} \quad \forall n>7 \end{align*} which can be put into: \begin{equation*} \forall k\geq 1, n\geq 2k,\quad a_k^{n} = (2i)^{n-2k}\frac{(n-k)!(n-(k+1))!}{(k-1)!k!(n-2k)!} \end{equation*}

The same kind of computation can be done for negative value of n, it appears that: \begin{align*} \forall n\geq 0, \quad Z_{-n} = \left(\frac{x-iy}{x^2+y^2}\right)^{n}\left[r^2 J_0(r) \sum_{k=1}^{u_n}a_k^{n} \left(x^2 + y^2\right)^{k-1} + i^n r J_1(r) \sum_{k=1}^{u_{n+1}}b_k^{n} \left(x^2 + y^2\right)^{k-1}\right] \end{align*} where $a_k^n$ and $b_k^n$ are defined by the same relations as before.