Space of solutions to a fourth order wave equation 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)$, where $a$ is a vector. For some choices of $a$ this new solution may also satisfy $\Box\left( a\cdot x f(x)\right) = 0,$ but I think it's always possible to choose $a$ such that $\Box\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$, which can be expanded in terms of a basis of functions which are orthogonal under some inner product? 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$ 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.
 A: You talk about the non-separability of the $\Box^2 \phi = 0$ equation, which I don't understand. Each plane wave $e^{ik\cdot x} = \prod_{j=0}^{d-1} e^{i k_j x^j}$ is already in separated form with the components of the null vector $k=(k_j)$ playing the role of the separation constants.
But rather than get into technicalities about what does or does not constitute a separable equation or a basis of separated solutions, why not just solve the equation via the Fourier transform?
Recall first the ordinary wave equation. If $\phi(x)$ is smooth and grows no faster than polynomially in any direction, its Fourier transform $\hat{\phi}(k)$ exists as a distribution (a subclass of tempered distributions in this case). The equation $\Box \phi(x) = 0$ translates to $k^2 \hat{\phi}(k) = 0$, where of course $k^2 = k_0^2 - \sum_{j=1}^{d-1} k_j^2$. No non-vanishing continuous, let alone smooth function of $k$ can satisfy that condition. Hence, $\hat{\phi}(k)$ must be a distribution. A natural candidate is $\hat{\phi}(k) = f(k) \delta(k^2)$, with $f(k)$ an at least continuous function, because of the simple identity $u \delta(u) = 0$. The function $f(k)$ is of course not unique, any functions $f_1,f_2$ such that the difference $f_1(k)-f_2(k)$ vanishes smoothly on the locus of $k^2$ define the same $\hat{\phi}(k)$. The level of regularity (including behavior at infinity) of $f(k)$ translates in a certain way to the level of regularity of $\phi(x)$ via the Fourier transform. The above parametrization of $\hat{\phi}(k)$ is actually exhaustive when these regularity classes are appropriately fixed. To see how the data in $f(k)$ translates to initial data for $\phi(x)$, try the $d=1$ example.
Now, on to the squared wave equation. Under the Fourier transform, $\Box^2 \phi(x) = 0$ translates to $(k^2) \hat{\phi}(k) = 0$. Taking derivatives of the $\delta$-function identity, we can get $u \delta'(u) + \delta(u) = 0$ or more importantly $u^2 \delta'(u)=0$. Hence, a natural candidate for a solution is
\begin{align*}
  \hat{\phi}(k) &= f(k) \delta(k^2) + g(k) \delta'(k^2) , \\
    &= F(k) \delta(k^2)
      + ia\cdot\partial_k [G(k) \delta(k^2)]
\end{align*}
for functions $f(k), g(k)$ of appropriate regularity and $G(k) = g(k)/(2ia\cdot k)$, $F(k) = f(k) - (a\cdot\partial_k) G(k)$. To make it easier to go between $g(k)$ and $G(k)$, the vector $a$ should not be null, so that $a\cdot k$ only vanishes at $k=0$. The corresponding formula in real space is
$$
  \phi(x) = \phi_1(x) + (a\cdot x) \phi_2(x) ,
$$
where $\phi_1(x), \phi_2(x)$ are independent solutions of the wave equation, $\Box \phi_{1,2}(x)=0$. This is I think the parametrization of solutions that you were looking for. As before, it is a matter of figuring out the right regularity classes for $F(k), G(k)$ and the corresponding initial data for $\phi(x)$ to make sure that the above parametrization is exhaustive.
