I'm interested in finding solutions a fourth order version of the standard wave equation in $d$ dimensional Minkowski spacetime $\mathcal{M}^d$. Defining $\Box := \partial_0^2 - \sum_{i = 1}^{d-1} \partial_{i}^2$, I want to find solutions to $$\Box^2\Phi(x) = 0,$$ which are not also solutions to $\Box \Phi(x) = 0$, and which transform as scalars under Lorentz transformations. I'm mostly interested in solutions for $d=4$, but I don't think that's likely to be relevant at the level of analysis I'm currently working at so I want to stick to general $d$ if possible.

I know that in solving $\Box \Phi(x) = 0$ with boundary conditions that require some kind of fast enough fall off at infinity, then the solution space can be taken to be Lorentz invariant $L^2$ integrable functions. A basis of solutions for functions with these boundary conditions can be given by plane waves $\Phi_k(x) := e^{i k\cdot x}$ where $k^2 = 0$. My understanding is that as $\Box \Phi(x) = 0$ is seperable, Sturm-Liouville theory tells us that for some choices of boundary conditions we can find a function space of solutions to the equation spanned by some basis of solutions.

(Disclaimer: I don't have a very in-depth understanding of analysis. I can see that $\Phi_k \not\in L^2(\mathcal{M}^d)$, so I understand that at some level calling this a `basis' for the space is not correct. I'm happy to call it this given that two such functions are orthogonal in the sense that $\int_{\mathcal{M}^d}dx \;\Phi_k(x)\Phi_{k'}(x) \propto \delta^d(k-k')$, and that any (I think?) function in this space can be decomposed in terms of an integral over all possible $\Phi_k$ by the Fourier transform. I appreciate that to answer my question it may be necessary to go to some more detailed level of analysis where it doesn't make sense to think of $\Phi_k$ as a basis.)

I note that, given any $f$ such that $\Box f(x) = 0,$ we can construct a solution to $\Box^2 \Phi(x) = 0$ by taking $\Phi(x) = a \cdot x f(x)$, for some vector $a$. For some choices of $a$ this new solution may also satisfy $\Box^2\left( a\cdot x f(x)\right) = 0,$ but I think it's always possible to choose $a$ such that $\Box^2\left( a\cdot x f(x)\right) \neq 0.$

Then my question is, given a function space of solutions corresponding to a choice of boundary conditions to $\Box \Phi(x) = 0$, is there a canonical way to extend this to a larger function space of soluitons to $\Box^2 \Phi(x) = 0$, which includes functions which satisfy $\Box^2 \Phi(x) = 0$ and $\Box \Phi(x) \neq 0$? If not a canonical way, is there some set of ways of doing this?

So then in the case of $L^2$ functions, is it possible to do something like extend to functions which increase like $x$ as $x\rightarrow \infty$? Then I could imagine that this new space could be spanned by $\Phi_k$, and $\chi_{k,a}(x) := a \cdot x e^{i k \cdot x}$, with some inner product under which $\Phi_k, \Phi_{k'}$, $\chi_{k,a}$ and $\chi_{k',a'}$ are orthogonal to each other? I think what I've written here is probably a little naive, but I'm hoping there may be some kind of result similar to this which makes sense.

After some searching online I've found that there exists a fourth order Sturm-Liouville theory, but I wasn't able to find anything accesible for me to read, and it appears to me that $\Box^2 \phi(x) = 0$ is anyway not seperable, at least in Cartersian coordinates. Any references to a simple introduction or review article of fourth-order Sturm-Lioville theory would be appreciated, as well as any local coordinate transformation which makes $\Box^2\phi(x) = 0$ seperable, if you think this would be relevant to my problem.

I appreciate that you could read this question and say, 'why would you expect there to be a canonical way to extend solution spaces for $\Box \Phi(x) = 0$ to $\Box^2 \Phi(x) = 0$ with boarder boundary conditions?'. So here is some evidence that I have from solving the equation in 2D.

In 2D in lightcone coordinates $u = x^0 + x^1, v = x^0-x^1$ the wave operator factorises so that $\Box = \partial_u\partial_v$. This makes it simple to write down d'Alembert's general solution to $\Box\Phi(x) = 0$ for all possible boudary conditions, $\Phi(x) = f^+(x^0 + x^1) + f^-(x^0-x^1).$

It's also simple to solve $\Box^2\Phi(x) = 0$ in this case. Writing a seperable ansatz $\Phi(u,v) = U(u)V(v)$ the equation becomes $U'' V''=0$, which is solved non-trivially by either $U' = 0$, $V' = 0$, $U''=0$ or $V''=0$. Then the general solution is given by $\Phi(u,v) = f^+_1(u) + v f^+_2(u) + f^-_1(v) + u f^-_2(v).$ Introducing null vectors $k_+ = k(1,1)$ and $k_- = q(1,-1)$, then we see that $u \propto k_+ \cdot x$, and $v \propto k_- \cdot x.$ Then the general solution can be written as a sum over solutions of the form $$\Phi(x) = f_1(k\cdot x) + a \cdot x f_2(k \cdot x),$$ where $k^2 = 0$ and $a$ is any vector. If $a \propto k_\pm$ then the second term solves $\Box\Phi(x) = 0$, otherwise it produces new solutions that satisfy $\Box^2\Phi(x) = 0$ and $\Box\Phi(x) \neq 0$. (I took a linear combination of the lightcone coordinate solutions to produce the new solution). I believe this to be the general solution to $\Box^2\Phi(x) = 0$ in 2D for all possible boundary conditions, and I've intentionally written it in a form that extends in a natural way to general dimension $d$. So this is some kind of motivation for why I think my question may have a solution; it appears to me that what I'm asking in general dimensions happens in 2D based on the form of the general solution I gave here.