# Symmetry operators for the conformal wave equation in rotating black hole
spacetimes

Finnian Gray<EMAIL_ADDRESS>Perimeter Institute, 31 Caroline
Street North, Waterloo, ON, N2L 2Y5, Canada  Department of Physics and
Astronomy, University of Waterloo, Waterloo, Ontario, Canada, N2L 3G1
Tsuyoshi Houri<EMAIL_ADDRESS>National Institute of Technology,
Maizuru College, Kyoto 625-8511, Japan  David Kubizňák
<EMAIL_ADDRESS>Perimeter Institute, 31 Caroline Street North,
Waterloo, ON, N2L 2Y5, Canada  Department of Physics and Astronomy, University
of Waterloo, Waterloo, Ontario, Canada, N2L 3G1  Yukinori Yasui
<EMAIL_ADDRESS>Institute for Fundamental Sciences, Setsunan
University, Osaka 572-8508, Japan

(April 29, 2021)

###### Abstract

We present covariant symmetry operators for the conformal wave equation in the
(off-shell) Kerr–NUT–AdS spacetimes. These operators, that are constructed
from the principal Killing–Yano tensor, its ‘symmetry descendants’, and the
curvature tensor, guarantee separability of the conformal wave equation in
these spacetimes. We next discuss how these operators give rise to a full set
of conformally invariant mutually commuting operators for the conformally
rescaled spacetimes and underlie the $R$-separability of the conformal wave
equation therein. Finally, by employing the WKB approximation we derive the
associated Hamilton–Jacobi equation with a scalar curvature potential term and
show its separability in the Kerr–NUT–AdS spacetimes.

##  I Introduction

Symmetries, both explicit and hidden, play an important role in general
relativity – in their presence one may be able to explicitly integrate the
Einstein equations and/or significantly simplify the study of matter fields in
a given curved spacetime. Perhaps one of the most remarkable symmetries is a
hidden symmetry of the principal Killing–Yano tensor Frolov:2017kze . Such a
symmetry appears for the Kerr family of spacetimes in all dimensions, or more
precisely for all the so called (off-shell) Kerr–NUT–AdS metrics Houri:2007xz
; Krtous:2008tb ; Houri:2008th , and underlies many of their remarkable
properties. In particular, it stands behind the separability of the massless
and massive scalar, spinor, and vector field equations in the Kerr–NUT–AdS
backgrounds Frolov:2006pe ; Oota:2007vx ; Lunin:2017drx ; Frolov:2018ezx  (see
also Lunin:2019pwz  for a separability of $p$-form fields).

Most recently, it has been demonstrated Gray:2020rtr  that also the
conformally coupled scalar wave equation

$\bigl{(}\Box-\eta R\bigr{)}\Phi=0\,,\quad\eta=\frac{1}{4}\frac{D-2}{D-1}\,,$
(1)

separates in the general off-shell Kerr–NUT–AdS spacetimes. Here, $D$ stands
for the number of spacetime dimensions, $R$ is the Ricci scalar of the
background metric ${\boldsymbol{g}}$, and prefactor $\eta$ is chosen so that
the equation enjoys conformal symmetry, (see, e.g., appendix D of
wald1984general  ). Namely, a solution to this equation remains a solution in
a conformally scaled spacetime

$\widetilde{{\boldsymbol{g}}}=\Omega^{2}{\boldsymbol{g}}\,,$ (2)

provided it also scales as $\widetilde{\Phi}=\Omega^{w}\Phi$ , with the
conformal weight $w=1-D/2$ . The wave equation equation (1) is of fundamental
importance and has a number of applications, see e.g. recent study of the
asymptotic structure of Kerr spacetime via conformal compactification
Hennig:2020rns .

The purpose of the present paper is to further our understanding of the
conformal wave equation (1) in the Kerr–NUT–AdS spacetime —filling some
important gaps in the previous analysis. In particular, we want to
‘intrinsically characterize’ the obtained separability by finding an explicit
covariant form of the corresponding symmetry operators that were found in
Gray:2020rtr  in a given coordinate basis. As we shall see, such operators can
be written in terms of the principal Killing–Yano tensor, its symmetry
descendants, and the curvature tensor. Moreover, following Michel:2013dfa ,
such operators can be ‘lifted up’ to conformal operators and guarantee
$R$-separability of the conformal wave equation in any conformally related
spacetime (2).

Finally, by applying the WKB approximation we derive an associated with (1)
Hamilton–Jacobi equation with a scalar curvature potential,

$g^{ab}\partial_{a}S\,\partial_{b}S+\eta R=0\,,$ (3)

and demonstrate its separability in the Kerr–NUT–AdS spacetimes. The equation
(3) has a long history, going back at least to a paper by DeWitt DeWitt:1952js
which considers quantum Hamiltonians arising from classical systems. Therein,
couplings to the geometrical objects can naturally arise. In a similar vein,
the extra term we find in the Hamiltonian can arise due to ambiguities in
operator ordering when quantizing non-linear systems Omote:1976fx . It has
also found use when considering the quantum mechanics of the motion of a free
particle constrained to a Riemannian surface Destri:1992sg ; Lian:2017qoh .
Here we understand it as a purely classical equation that describes certain
modification of the free particle motion in a curved space.

Our plan for the remainder of the paper is as follows. In the next section we
review the Kerr–NUT–AdS spacetimes, their hidden symmetry of the principal
Killing–Yano tensor, and its ‘symmetry descendants’. In Sec. III we construct
the covariant form of the symmetry operators for the conformal wave equation
in these spacetimes. The associated operators for the conformally rescaled
metrics are studied in Sec. IV. In Sec. V we derive the Hamilton–Jacobi
equation (3) and demonstrate its separability in Kerr–NUT–AdS spacetimes. Sec.
VI is devoted to the final discussion. Technical results are summarized in
Appendices A and B.

##  II Principal Killing–Yano tensor and Kerr–NUT–AdS spacetimes

The principal Killing–Yano tensor ${\boldsymbol{h}}$ is a non-degenerate
closed conformal Killing–Yano 2-form ${\boldsymbol{h}}$ obeying the following
equation:

$\nabla_{a}h_{bc}=g_{ab}\,\xi_{c}-g_{ac}\,\xi_{b}\,,$ (4)

where

$\xi^{a}=\frac{1}{D-1}\nabla_{b}h^{ba}\,$ (5)

is the associated primary Killing vector field Krtous:2008tb .

Starting with a single principal Killing–Yano tensor ${\boldsymbol{h}}$, one
can generate the whole towers of explicit and hidden symmetries – the
‘symmetry descendants’ of ${\boldsymbol{h}}$. In brief, we can construct the
following tower of closed conformal Killing–Yano tensors:

${\boldsymbol{h}}^{(j)}=\frac{1}{j!}\underbrace{{\boldsymbol{h}}\wedge\dots\wedge{\boldsymbol{h}}}_{j\
\mbox{\tiny times}}\,.$ (6)

Their Hodge duals ${\boldsymbol{f}}^{(j)}=\star{\boldsymbol{h}}^{(j)}$ are
Killing–Yano tensors, and their square gives rise to a tower of rank-2 Killing
tensors:

$k^{ab}_{(j)}=\frac{1}{(D-2j-1)!}f^{(j)a}{}_{c_{1}\dots
d_{d-2j-1}}f^{(j)bc_{1}\dots c_{d-2j-1}}\,$ (7)

for $j\in(0,\dots,n-1)$. In turn, these tensors give rise to the tower of
Killing vectors:

$\boldsymbol{l}_{(j)}=\boldsymbol{k}_{(j)}\cdot\boldsymbol{\xi}^{\flat}\,.$
(8)

Note that the $j=0$ Killing tensor is just the inverse metric and the zeroth
Killing vector is the primary Killing vector,
$\boldsymbol{l}_{(0)}={\boldsymbol{\xi}}$. We also have in odd dimensions an
extra redundant Killing tensor
${\boldsymbol{k}}_{(n)}={\boldsymbol{l}}_{(n)}\otimes{\boldsymbol{l}}_{(n)}$.

All of the above constructed symmetries mutually Schouten–Nijenhuis commute

$\displaystyle\left[\boldsymbol{l}_{(i)},\boldsymbol{k}_{(j)}\right]_{\mbox{\tiny
SN}}$ $\displaystyle=$ $\displaystyle
0\;,\;\left[\boldsymbol{l}_{(i)},\boldsymbol{l}_{(j)}\right]_{\mbox{\tiny
SN}}=0\;,$
$\displaystyle\left[\boldsymbol{k}^{(i)},\boldsymbol{k}^{(j)}\right]_{\mbox{\tiny
SN}}$ $\displaystyle=$ $\displaystyle
k^{(i)}_{e(a}\nabla^{e}k^{(j)}_{bc)}-k^{(j)}_{e(a}\nabla^{e}k^{(i)}_{bc)}=0\;\,.$
(9)

In addition, the Killing tensors obey the following algebraic identity (i.e.
they commute as matrices):

$k^{a}_{(i)\,b}k^{b}_{(j)\,c}-k^{a}_{(j)\,b}k^{b}_{(i)\,c}=0\,,$ (10)

see Frolov:2017kze  for all the details and proofs of the above statements.

The most general spacetime admitting the principal Killing–Yano tensor is the
(off-shell) Kerr–NUT–AdS spacetime Houri:2007xz ; Krtous:2008tb  (see also
Houri:2008th ). Denoting by ${D=2n+\varepsilon}$ the total number of spacetime
dimensions (with $\varepsilon=0$ in even and $\varepsilon=1$ in odd
dimensions), the metric takes the following explicit form:

$\displaystyle{\boldsymbol{g}}$ $\displaystyle=$
$\displaystyle\sum_{\mu=1}^{n}\;\biggl{[}\;\frac{U_{\mu}}{X_{\mu}}\,{{{\boldsymbol{d}}}x_{\mu}^{2}}+\,\frac{X_{\mu}}{U_{\mu}}\,\Bigl{(}\,\sum_{j=0}^{n-1}A^{\\!(j)}_{\mu}{{\boldsymbol{d}}}\psi_{j}\Bigr{)}^{\\!2}\;\biggr{]}$
(11) $\displaystyle\qquad+\frac{\varepsilon
c}{A^{\\!(n)}}\Bigl{(}\sum_{k=0}^{n}A^{\\!(k)}{{\boldsymbol{d}}}\psi_{k}\\!\Bigr{)}^{\\!2}\;,$

while the principal Killing–Yano tensor reads

${\boldsymbol{h}}=\sum_{\mu=1}^{n}\,x_{\mu}dx^{\mu}\wedge\left(\sum_{k=0}^{{{n}}-1}A^{\\!(k)}_{\mu}{{\boldsymbol{d}}}\psi_{k}\right)\;.$
(12)

The employed coordinates $\\{x_{\mu},\psi_{k}\\}$ have a natural geometrical
meaning associated with the principal Killing–Yano tensor. They split into
(time and azimuthal angle) _Killing coordinates_ ${\psi_{k}}$
(${k}={0,\,\dots,\,{{n}}{-}1{+}\varepsilon}$) that correspond to the Killing
vectors (8),

$\boldsymbol{l}_{(k)}={{\boldsymbol{\partial}}}_{\psi_{k}}\,,$ (13)

and the non-trivial (radial and longitudinal angle) coordinates ${x_{\mu}}$
($\mu=1,\,\dots,\,{{n}}$) that represent the ‘eigenvalues’ of
${\boldsymbol{h}}$, see Frolov:2017kze .

In the above, the functions ${A^{\\!(k)}}$, ${A^{\\!(j)}_{\mu}}$, and
${U_{\mu}}$ are ‘symmetric polynomials’ of the coordinates ${x_{\mu}}$, and
are defined by:

$\displaystyle A^{\\!(k)}$
$\displaystyle=\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}\nu_{1},\dots,\nu_{k}=1\\\
\nu_{1}<\dots<\nu_{k}\end{subarray}}^{{n}}\\!\\!\\!\\!\\!x^{2}_{\nu_{1}}\dots
x^{2}_{\nu_{k}}\;,\>\>\>A^{\\!(j)}_{\mu}=\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}\nu_{1},\dots,\nu_{j}=1\\\
\nu_{1}<\dots<\nu_{j}\\\
\nu_{i}\neq\mu\end{subarray}}^{{n}}\\!\\!\\!\\!\\!x^{2}_{\nu_{1}}\dots
x^{2}_{\nu_{j}}\;,$ $\displaystyle U_{\mu}$
$\displaystyle=\prod_{\begin{subarray}{c}\nu=1\\\
\nu\neq\mu\end{subarray}}^{{n}}(x_{\nu}^{2}-x_{\mu}^{2})\;,\;\>\>U=\prod_{\begin{subarray}{c}\mu,\nu=1\\\
\mu<\nu\end{subarray}}^{n}(x^{2}_{\mu}-x^{2}_{\nu})=\det\bigl(A^{\\!(j)}_{\mu}\bigr{missing})\;,$
(14)

where we have fixed $A^{\\!(0)}=1=A^{\\!(0)}_{\mu}$. Each metric function
${X_{\mu}}$ is an unspecified function of a single coordinate ${x_{\mu}}$:

$X_{\mu}=X_{\mu}(x_{\mu})\,.$ (15)

Lastly, the constant $c$ only appearing in odd dimensions is a free parameter.

Despite the fact that the metric is rather complex, its Ricci scalar takes a
fairly simple form Hamamoto:2006zf

$R=\sum_{\mu=1}^{n}\frac{r_{\mu}}{U_{\mu}}\,,$ (16)

where each function $r_{\mu}$ depends only on a single variable $x_{\mu}$:

$r_{\mu}=-X_{\mu}^{\prime\prime}-\frac{2\varepsilon
X_{\mu}^{\prime}}{x_{\mu}}-\frac{2\varepsilon c}{x_{\mu}^{4}}\,.$ (17)

The determinant of the metric reads

$\sqrt{\absolutevalue{g}}=\bigl{(}cA^{\\!(n)}\bigr{)}^{\frac{\varepsilon}{2}}\,U\,.$
(18)

Importantly for our purposes the Killing tensors ${\boldsymbol{k}}_{(j)}$ take
the following coordinate form:

$\displaystyle{\boldsymbol{k}}_{(j)}\\!$ $\displaystyle=$
$\displaystyle\\!\sum_{\mu=1}^{n}\\!A^{\\!(j)}_{\mu}\\!\\!\left[\\!\frac{X_{\mu}}{U_{\mu}}\,{{{\boldsymbol{\partial}}}_{x_{\mu}}^{2}}\\!+\\!\frac{U_{\mu}}{X_{\mu}}\\!\left(\\!\sum_{k=0}^{n-1+\varepsilon}\\!{\frac{(-x_{\mu}^{2})^{n-1-k}}{U_{\mu}}}{{\boldsymbol{\partial}}}_{\psi_{k}}\right)^{\\!\\!2}\\!\right]$
(19)
$\displaystyle+\varepsilon\,\frac{A^{\\!(j)}}{cA^{\\!(n)}}{{\boldsymbol{\partial}}}_{\psi_{n}}^{2}\,,$

where $j=0$ corresponds to the inverse metric,
${\boldsymbol{g}}^{-1}={\boldsymbol{k}}_{(0)}$.

##  III Separability of the conformal wave equation and its intrinsic
characterization

Recently it was shown Gray:2020rtr , that a solution to the conformal wave
equation (1) in the background (11) can be found in the multiplicative
separated form,

$\Phi=\prod_{\mu=1}^{n}Z_{\mu}(x_{\mu})\prod_{k=0}^{n-1+\varepsilon}e^{i\Psi_{k}\psi_{k}}\,,$
(20)

where $\Psi_{k}$ are the Killing vector separation constants, and each of the
$Z_{\mu}$, which is a function of the single corresponding variable $x_{\mu}$,
obeys the following ordinary differential equation:

$\displaystyle{Z_{\mu}^{\prime\prime}}$
$\displaystyle+{Z_{\mu}^{\prime}}\Bigl{(}\frac{X_{\mu}^{\prime}}{X_{\mu}}+\frac{\varepsilon}{x_{\mu}}\Bigr{)}-\frac{Z_{\mu}}{X_{\mu}^{2}}\Bigl{(}\sum_{k=0}^{n-1+\varepsilon}{(-x_{\mu}^{2})^{n-1-k}}\Psi_{k}\Bigr{)}^{2}$
(21) $\displaystyle-\frac{Z_{\mu}}{X_{\mu}}\Bigl{(}\eta
r_{\mu}+\frac{\varepsilon}{cx_{\mu}^{2}}\Psi_{n}^{2}+\sum_{k=0}^{n-1}C_{k}(-x_{\mu}^{2})^{n-1-k}\Bigr{)}=0\,,\quad\
\ $

where $C_{k}$ ($k=0,\dots,n-1)$ are the (non-trivial) separation constants and
we have set $C_{0}=0$.

As also shown in Gray:2020rtr , underlying this separability is a complete set
of symmetry operators $\\{{\cal K}_{(j)},{\cal L}_{(j)}\\}$,

$\displaystyle{\cal K}_{(j)}$ $\displaystyle=$
$\displaystyle\nabla_{a}k_{(j)}^{ab}\nabla_{b}-\eta R_{(j)}\,,$ (22)
$\displaystyle{\cal L}_{(j)}$ $\displaystyle=$
$\displaystyle-i\,l_{(j)}^{a}\nabla_{\\!a}\,,$ (23)

that all mutually commute one of another,

$\bigl{[}{\cal K}_{(k)},{\cal L}_{(l)}\bigr{]}=0\;,\;\bigl{[}{\cal
L}_{(k)},{\cal L}_{(l)}\bigr{]}=0\;,\;\bigl{[}{\cal K}_{(k)},{\cal
K}_{(l)}\bigr{]}=0\;,$ (24)

and one of which is the wave conformal operator. Namely,

${\cal K}_{(0)}\Psi=0$ (25)

is the conformal wave equation (1). The fact that these commuting operators
exist means that, there exists a common eigenfunction of these operators
$\Phi$ obeying

$\displaystyle{\cal K}_{(j)}\Phi$ $\displaystyle=$ $\displaystyle
C_{j}\Phi\,,$ (26) $\displaystyle{\cal L}_{(j)}\Phi$ $\displaystyle=$
$\displaystyle\Psi_{j}\Phi\,.$ (27)

It is precisely this eigenfunction which is the separated solution (20).

The operators ${\cal L}_{(j)}$ are the standard scalar operators that are
generated from Killing vectors $\boldsymbol{l}_{(j)}$. On the other hand, the
Killing tensor operators ${\cal K}_{(j)}$ pick up, in addition to the standard
Killing tensor part $\nabla_{a}k_{(j)}^{ab}\nabla_{b}$, also an ‘anomalous
conformal term’ $R_{(j)}$ which ensures the commutation with the conformal
wave operator ${\cal K}_{(0)}$. In Gray:2020rtr  an explicit coordinate
expression for this term has been found, it reads

$R_{(j)}=\sum_{\mu=1}^{n}\frac{A_{\mu}^{(j)}}{U_{\mu}}r_{\mu}\,,$ (28)

where $r_{\mu}$ are the ‘Ricci scalar functions’ (17). However, no covariant
expression for $R_{(j)}$ has been given in Gray:2020rtr .

Here we amend this situation. That is to say, we show in the appendix A that
$R_{(j)}$ are given in terms of the principal Killing–Yano tensor, its
symmetry descendants, and the curvature tensor by the following covariant
formula:

$\displaystyle R_{(j)}$
$\displaystyle=k^{ab}_{(j)}R_{ab}+\frac{D-4}{2(D-2)}\Box\Tr({\boldsymbol{k}}_{(j)})$
$\displaystyle\qquad+\alpha_{j}k_{(j-1)}^{ac}h_{c}{}^{b}\,(d\xi)_{ab}-\beta_{j}l_{(j-1)}^{a}\,\xi_{a}$
$\displaystyle=k^{ab}_{(j)}R_{ab}+\frac{D-4}{2(D-2)}\Box\Tr({\boldsymbol{k}}_{(j)})$
$\displaystyle\qquad-
k_{(j-1)}^{ab}\Bigl{(}\alpha_{j}h_{a}{}^{c}\,(d\xi)_{cb}+\beta_{j}\xi_{a}\,\xi_{b}\Bigr{)}\;,$
(29)

where ${\boldsymbol{\xi}}={\boldsymbol{l}}_{(0)}$ is the primary Killing
vector (5), for $j=0$ we defined ${\boldsymbol{k}}_{(-1)}\equiv
0\equiv{\boldsymbol{l}}_{(-1)}$, and the constants $\alpha_{j}$ and
$\beta_{j}$ are given by

$\displaystyle\alpha_{j}$
$\displaystyle=\frac{(n-j+\frac{\epsilon}{2})}{(n-1+\frac{\epsilon}{2})}\,\quad\beta_{j}=2\frac{(n-j+\frac{\epsilon}{2})}{(n-1+\frac{\epsilon}{2})}(2j-3)\,.$
(30)

Interestingly these objects can be understood as follows. Let us define the
following 1-forms ${\boldsymbol{\kappa}}^{(j)}$:

$\kappa^{(j)}_{a}=k_{(j)\,a}^{\;\;\;\;\;\;\;b}\nabla_{b}R\;.$ (31)

Then, quantities $R_{(j)}$ can be understood as ‘potentials’ for the above
1-forms:

${\boldsymbol{\kappa}}_{(j)}={\boldsymbol{d}}R_{(j)}\;,$ (32)

see appendix A for the proof. In fact, it is this property which underlies the
commutation of the operators (22). Given that
$[\nabla_{a}k_{(i)}^{ab}\nabla_{b},\nabla_{c}k_{(j)}^{cd}\nabla_{d}]f=0$
Sergyeyev:2007gf ; Gray:2020rtr  for any scalar function $f$, we have

$\displaystyle\left[\mathcal{K}_{(i)},\mathcal{K}_{(j)}\right]f\\!\\!\\!\\!\\!\\!\\!\\!$
$\displaystyle=-\eta$
$\displaystyle\left(\left[\nabla_{a}k_{(i)}^{ab}\nabla_{b},R_{(j)}\right]f+\left[R_{(i)},\nabla_{a}k_{(j)}^{ab}\nabla_{b}\right]f\right)$
$\displaystyle=-\eta$
$\displaystyle\left\\{\nabla_{a}(k_{(i)}^{ab}\nabla_{b}(R_{(j)}f))-R_{(j)}\nabla_{a}(k_{(i)}^{ab}\nabla_{b}(f))\right.$
$\displaystyle\left.+R_{(i)}\nabla_{a}(k_{(i)}^{ab}\nabla_{b}(f))-\nabla_{a}(k_{(j)}^{ab}\nabla_{b}(R_{(i)}f))\right\\}$
$\displaystyle=-\eta$
$\displaystyle\left\\{\nabla_{a}\Big{(}f\,[k^{ab}_{(i)}\nabla_{b}R_{(j)}-k^{ab}_{(j)}\nabla_{b}R_{(i)}]\Big{)}\right.$
$\displaystyle+$
$\displaystyle\left.(\nabla_{a}f)\Big{(}k^{ab}_{(i)}\nabla_{b}R_{(j)}-k^{ab}_{(j)}\nabla_{b}R_{(i)}\Big{)}\right\\}$
$\displaystyle=-\eta$
$\displaystyle\left\\{\nabla_{a}\Big{(}f\,[k^{a}_{(i)\,b}k^{b}_{(j)\,c}-k^{a}_{(j)\,b}k^{b}_{(i)\,c}]\nabla^{c}R\Big{)}\right.$
$\displaystyle+$
$\displaystyle\left.(\nabla_{a}f)(k^{a}_{(i)\,b}k^{b}_{(j)\,c}-k^{a}_{(j)\,b}k^{b}_{(i)\,c})\nabla^{c}R\right\\}$
$\displaystyle=0\ \ $ $\displaystyle\,,$ (33)

where in the final step we have used the algebraic identity (10).

We note that, this is a special case of the result presented in benenti2002 .
Therein, it is shown that the commutation of any operators, $\Box+g$ and
$\nabla_{a}K^{ab}\nabla_{b}+f$, where $f,g\in{\cal C}^{\infty}({\cal M})$ and
$K^{ab}$ is a Killing tensor is guaranteed provided

$\nabla_{a}f=\tensor{K}{{}_{a}^{b}}\nabla_{a}g-\frac{1}{3}\nabla_{b}(\tensor{K}{{}_{a}^{c}}\tensor{R}{{}_{c}^{b}}-\tensor{R}{{}_{a}^{c}}\tensor{K}{{}_{c}^{b}})\,.$
(34)

In the case of the off-shell Kerr–NUT–AdS metrics the final term on the right
hand side vanishes as the Killing and Ricci tensors are diagonal in the same
basis Frolov:2017kze ; Hamamoto:2006zf  (See (82) and (84) in Appendix A).
Thus, this equation reduces to the relationship between (31) and (32).

##  IV Symmetry operators in conformally related spacetimes

As mentioned in the introduction, the conformal wave equation (1) enjoys the
conformal symmetry. That is, provided we have a solution $\Phi$ in the
spacetime ${\boldsymbol{g}}$, then

$\widetilde{\Phi}=\Omega^{w}\Phi\,,\quad w=1-D/2\,$ (35)

is a solution of the same equation in the conformally rescaled spacetime

$\widetilde{{\boldsymbol{g}}}=\Omega^{2}{\boldsymbol{g}}\,.$ (36)

In particular, this means that (35) with $\Phi$ given by (20) yields an
$R$-separated solution of the conformal wave equation in any spacetime related
to the off-shell Kerr–NUT–AdS metric by the conformal transformation (36).

It is interesting to ask if also such $R$-separability can be intrinsically
characterized by some complete set of mutually commuting operators. In what
follows we show that this is indeed the case – we explicitly construct such
operators and discuss their properties. First, starting from the special
conformal frame with $\Omega=1$, we scale the operators $\\{{\cal
K}_{(j)},{\cal L}_{(j)}\\}$, to construct a complete set of mutually commuting
operators for the metric $\widetilde{{\boldsymbol{g}}}$, (36). Second,
following Michel:2013dfa , we show that such operators can in fact be lifted
to conformally invariant operators, providing thus a complete set of
conformally invariant mutually commuting operators for the conformal wave
equation (1) in any spacetime related to the Kerr–NUT–AdS metric by a
conformal transformation.

### Mutually commuting operators

Starting from the mutually commuting operators $\\{{\cal K}_{(j)},{\cal
L}_{(j)}\\}$ in the special frame with $\Omega=1$, let us define new operators
$\\{\widetilde{\mathcal{O}}_{(j)},\widetilde{\mathcal{P}}_{(j)}\\}$ for
general $\Omega$ by:

$\displaystyle\widetilde{\mathcal{O}}_{(j)}$ $\displaystyle\equiv$
$\displaystyle\Omega^{w}{\cal K}_{(j)}\Omega^{-w}\,,$
$\displaystyle\widetilde{\mathcal{P}}_{(j)}$ $\displaystyle\equiv$
$\displaystyle\Omega^{w}{\cal L}_{(j)}\Omega^{-w}\,.$ (37)

By construction such operators mutually commute, as we have

$\displaystyle\left[\widetilde{\mathcal{O}}_{(i)},\widetilde{\mathcal{O}}_{(j)}\right]$
$\displaystyle=\Omega^{w}\left[\mathcal{K}_{(i)},\mathcal{K}_{(j)}\right]\Omega^{-w}=0\;,$
(38)
$\displaystyle\left[\widetilde{\mathcal{O}}_{(i)},\widetilde{\mathcal{{P}}}_{(j)}\right]$
$\displaystyle=\Omega^{w}\left[\mathcal{K}_{(i)},\mathcal{L}_{(j)}\right]\Omega^{-w}=0\,,$
(39)
$\displaystyle\left[\widetilde{\mathcal{P}}_{(i)},\widetilde{\mathcal{P}}_{(j)}\right]$
$\displaystyle=\Omega^{w}\left[\mathcal{L}_{(i)},\mathcal{L}_{(j)}\right]\Omega^{-w}=0\,.$
(40)

Moreover, it follows that when $\Phi$ satisfies the eigenvalue problem (26) in
the spacetime ${\boldsymbol{g}}$, $\widetilde{\Phi}=\Omega^{w}\Phi$ given by
(35) obeys the ‘associated’ eigenvalue problem:

$\displaystyle\widetilde{\mathcal{O}}_{(j)}\widetilde{\Phi}$ $\displaystyle=$
$\displaystyle C_{j}\widetilde{\Phi}\,,$
$\displaystyle\widetilde{\mathcal{P}}_{(j)}\widetilde{\Phi}$ $\displaystyle=$
$\displaystyle\Psi_{j}\widetilde{\Phi}\,,$ (41)

in the conformal spacetime $\widetilde{{\boldsymbol{g}}}$. In other words, the
operators $\\{\widetilde{\mathcal{O}}_{(j)},\widetilde{\mathcal{P}}_{(j)}\\}$,
(IV), intrinsically characterize the separability of the conformal wave
equation in the conformal spacetime (36).

The only ‘problem’ with (IV) is that the new operators
$\\{\widetilde{\mathcal{O}}_{(j)},\widetilde{\mathcal{P}}_{(j)}\\}$ remain
expressed in terms of the ‘old’ connection $\nabla_{a}$, the old Ricci tensor
$R_{ab}$, and other objects associated with the metric ${\boldsymbol{g}}$
rather than the conformally rescaled metric $\widetilde{{\boldsymbol{g}}}$.
However, using the well known transformation properties of the connection and
curvature tensor, one can straightforwardly amend this situation.

For example, let us define the following tilded objects:111We stress that
these objects are not the conformal symmetries of the spacetime
$\widetilde{{\boldsymbol{g}}}$, although it is possible to define such
symmetries. Namely, the following objects:
$k^{ab}_{(j>0)}\,,\quad\Omega^{3}h_{ab}\,,\quad l^{a}_{(j\geq 0)}\,,$ are the
conformal Killing tensors, conformal Killing–Yano 2-form, and conformal
Killing vectors of the spacetime $\widetilde{{\boldsymbol{g}}}$. Notice that
in doing this, necessarily ${\boldsymbol{k}}_{(0)}={\boldsymbol{g}}$
transforms differently to the rest of the Killing tensors. One could, of
course, use these objects to define the transformed operators, leading to
different (seemingly more complex) expressions. We will adopt this strategy
for the Killing tensors at least in the next section (IV).

$\widetilde{k}^{ab}_{(j)}=\Omega^{-2}k^{ab}_{(j)}\,,\quad\widetilde{h}_{ab}=\Omega^{2}h_{ab}\,,\quad\widetilde{l}^{a}_{(j)}=\Omega^{-2}l^{a}_{(j)}\,,$
(42)

and raise or lower their indices with the metric
$\widetilde{{\boldsymbol{g}}}$ and its inverse. We further denote by
$\widetilde{\nabla}_{a}$ the covariant derivative in the spacetime
$\widetilde{{\boldsymbol{g}}}$ and by $\widetilde{R}_{ab}$ its Ricci tensor.
With these at hand, the operators (IV) can be expressed as follows (see
appendix B for details):

$\displaystyle\widetilde{\mathcal{O}}_{(j)}:=$
$\displaystyle\Omega^{2}\left(\widetilde{\mathcal{K}}_{(j)}+\eta\left[\Bigl{(}\widetilde{\nabla}_{a}\widetilde{\nabla}_{b}\bigl{(}\widetilde{k}_{(j)}^{ab}+\frac{1}{2}\widetilde{k}^{c}_{(j)\,c}\widetilde{g}^{ab}\bigr{)}\Bigr{)}\right]\right)\,,$
(43) $\displaystyle\widetilde{\mathcal{P}}_{(j)}:=$
$\displaystyle\Omega^{2}\left(\widetilde{\mathcal{L}}_{(j)}-\frac{w}{D-2}\widetilde{\nabla}_{a}\widetilde{l}^{a}_{(j)}\right)\;,$
(44)

where $\widetilde{\mathcal{K}}_{(j)}$ and $\widetilde{\mathcal{L}}_{(j)}$ are
given by expressions (22), (23), and (III), with all the objects replaced by
the tilded ones. Note that the quantities
$\widetilde{\nabla}_{b}\left[\widetilde{k}_{(j)}^{ab}+\frac{1}{2}\widetilde{k}^{c}_{(j)\,c}\widetilde{g}^{ab}\right]$
and $\widetilde{\nabla}_{a}\widetilde{l}^{a}_{(j)}$ vanishes identically when
$\Omega=1$ due to the Killing tensor and Killing vector equations

$\nabla\tensor{[}_{(a}]{k}{{}^{(j)}_{b}{}_{c)}}=0\,,\quad\nabla_{\\!(a}l^{(j)}_{b)}=0\,,$
(45)

respectively.

Moreover, $\widetilde{\mathcal{O}}_{(0)}$ is just a conformally rescaled
$\widetilde{\mathcal{K}}_{(0)}$,

$\widetilde{\mathcal{K}}_{(0)}=\Omega^{-2}\widetilde{\mathcal{O}}_{(0)}=\Omega^{w-2}\mathcal{K}_{(0)}\Omega^{-w}\,,$
(46)

highlighting the conformal invariance of this operator. The other operators,
however, take a more complicated form, as is to be expected from the
privileged role of the conformal frame with $\Omega=1.$222This is the only
frame where the spacetime admits full (not only conformal) Killing tensors and
the Ricci tensor is diagonal in the natural orthonormal frame Kubiznak:2007kh
. We shall return to this issue in the next subsection where we discuss the
conformal form of these operators.

### Conformal symmetry operators

Conformal symmetry operators for the conformal wave equation have been studied
for many years, see e.g. carter1977killing ; boyer1976symmetry ;
kalnins1982intrinsic ; kamran1985separation2 ; duval1999conformally ;
benenti2002 ; eastwood2005higher ; eastwood2008higher ; gover2012higher ;
andersson2014second . This work culminated in ref. Michel:2013dfa  where a
complete and constructive theory was finally formulated. Our goal for the
remainder of this section is to review this theory in a more physics community
oriented language, and briefly discuss how it applies to the problem at hand.

To start with, we define a conformally invariant operator as an operator that
preserves its form under a conformal transformation. More specifically, a
conformally invariant operator of weights $s_{1}$ and $s_{2}$ obeys the
following equality:

$\widetilde{Q}_{s_{1},s_{2}}=\Omega^{s_{2}}Q_{s_{1},s_{2}}\Omega^{-s_{1}}\,,$
(47)

under the conformal transformation (36). That is,
$\widetilde{Q}_{s_{1},s_{2}}$ has exactly the ‘same form’ as
${Q}_{s_{1},s_{2}}$ but is constructed out of conformally scaled (tilded)
tensors associated with the metric $\widetilde{{\boldsymbol{g}}}$ rather than
${\boldsymbol{g}}$. To give an example, the conformal wave operator
$\mathcal{K}_{(0)}$ obeys the equation (46) and thence is a conformal operator
with weights $s_{1}=w$ and $s_{2}=w-2$.

In what follows, we are going to concentrate on conformal operators of equal
weights, $s_{1}=s_{2}=s$. In particular, as shown in Michel:2013dfa  the most
general second-order conformal operator with weight $s$ that is built out of a
symmetric tensor $K^{ab}$ is given by

$\displaystyle Q_{s}(K)=$
$\displaystyle\nabla_{a}K^{ab}\nabla_{b}+\Bigl{(}\gamma_{1}[\nabla_{a}K^{ab}]+\gamma_{2}[\nabla^{b}\Tr
K]\Bigr{)}\nabla_{b}$
$\displaystyle+\gamma_{3}(\nabla_{a}\nabla_{b}K^{ab})+\gamma_{4}(\Box\Tr
K)+\gamma_{5}\,R_{ab}K^{ab}$ $\displaystyle+\gamma_{6}\,R\,\Tr K+f\;.$ (48)

Here $f$ is a function which does not scale under conformal transformation, we
assume $\widetilde{K}^{ab}=K^{ab}$, and the coefficients are

$\displaystyle\gamma_{1}$
$\displaystyle=2\gamma_{2}=-\frac{(2s+D)}{D+2}\,,\,\gamma_{3}=\frac{(s-1)s}{(D+1)(D+2)}\,,$
$\displaystyle\gamma_{4}$
$\displaystyle=\frac{s(D+2s-1)}{2(D+1)(D+2)}\,,\,\gamma_{5}=\frac{s(D+s)}{(D-2)(D+1)}\,,$
$\displaystyle\gamma_{6}$
$\displaystyle=-\frac{2s(D+s)}{(D-2)(D-1)(D+1)(D+2)}\,.$ (49)

Similarly, having a vector $L^{a}$, the corresponding conformal operator is
given by

$Q_{s}(L)=L^{a}\nabla_{\\!a}-\frac{s}{D}(\nabla_{\\!a}\,L^{a})\;.$ (50)

In particular, we consider conformal operators of weight $w=1-D/2$, c.f. (35),

$\widetilde{Q}_{w}=\Omega^{w}Q_{w}\Omega^{-w}\,,$ (51)

that are symmetry operators of the conformal wave operator
$\mathcal{K}_{(0)}$, that is, they satisfy the following relation:

$\mathcal{K}_{(0)}\circ Q_{w}={\cal D}\circ\mathcal{K}_{(0)}\,,$ (52)

for some operator ${\cal D}$; in fact, it is easy to see that the conformal
invariance implies ${\cal D}\equiv{\cal D}_{-2+w}$. Note that the equation
(52) obviously preserves the kernel of $\mathcal{K}_{(0)}$.

To find such symmetry operators we can use the following theorem
Michel:2013dfa :
Theorem 1. _Let $K^{ab}$ be a (special) Killing tensor of the metric
${\boldsymbol{g}}$, so that the following conformally invariant ‘geometric
obstruction’ built from the Weyl tensor $C_{abcd}$:_

$\text{Obs}(K)_{a}=\frac{2(D-2)}{3(D+1)}\Bigl{(}\nabla_{b}K^{cd}\tensor{C}{{}^{b}_{c}{}_{d}{}_{a}}-\frac{3}{D-3}K^{cd}\nabla_{b}\tensor{C}{{}^{b}_{c}{}_{d}{}_{a}}\Bigr{)}$
(53)

_is exact, that is,_

${\bf Obs}(K)=-2{\boldsymbol{d}}f\,.$ (54)

_Then ( IV) with $f$ given by (54) (up to a constant) is a symmetry operator
for the conformal wave operator and in fact satisfies_

$\mathcal{K}_{(0)}\circ Q_{w}(K)=Q_{-2+w}(K)\circ\mathcal{K}_{(0)}\,.$ (55)

When $K^{ab}$ is a Killing tensor we can simplify the operator (IV) via the
Killing equation,

$\nabla\tensor{[}_{(a}]{K}{{}_{b}{}_{c)}}=0\,,$ (56)

however this will only hold for a particular metric of the conformal class.
For this particular metric, we then have

$\displaystyle Q_{w}(K)$ $\displaystyle=$ $\displaystyle Q_{w-2}(K)$ (57)
$\displaystyle=$
$\displaystyle\nabla_{a}K^{ab}\nabla_{b}-\frac{(D-2)}{8(D+1)}\left[\Box\Tr
K\right]$ $\displaystyle-$
$\displaystyle\frac{(D+2)}{4(D+1)}R_{ab}K^{ab}+\frac{R\,\Tr
K}{2(D+1)(D-1)}+f\;.\quad\quad$

In this case, therefore the corresponding symmetry operator (52) actually
commutes with the conformal wave equation

$\left[Q_{w},\mathcal{K}_{(0)}\right]=0\,,$ (58)

and more generally, we have the conformal commutation

$\left[\widetilde{Q}_{w},\,\Omega^{2}\,\widetilde{\mathcal{K}}_{(0)}\right]=0\,,$
(59)

valid in any conformal frame.

In particular, taking the Killing tensors ${\boldsymbol{k}}_{(j)}$ ($j>0)$ in
the Kerr–NUT–AdS metric ${\boldsymbol{g}}$, we find that they satisfy the
obstruction condition (54) with $f_{(j)}$ given by

$\displaystyle f_{(j)}$ $\displaystyle=$
$\displaystyle\frac{1}{4(1-D^{2})}\Bigl{[}2D\,k^{ab}_{(j)}R_{ab}+3\Box\Tr
k_{(j)}$ (60)
$\displaystyle+(D+1)(D-2)k_{(j-1)}^{ab}\Bigl{(}\alpha_{j}h_{a}{}^{c}\,(d\xi)_{cb}+\beta_{j}\xi_{a}\,\xi_{b}\Bigr{)}$
$\displaystyle-2R\,\Tr k_{(j)}\Bigr{]}\,.$

It can then easily be checked that333Of course, the expression (53) is only
defined in this coordinate invariant way in the $\Omega=1$ frame although its
coordinate expression will be unchanged no matter the frame. the corresponding
operators

${\cal K}_{w}^{(j)}\equiv Q_{w}(k_{(j)})\,,$ (61)

(IV), coincide with the operators $\mathcal{K}_{(j)}$, (22),

${\cal K}_{w}^{(j)}=\mathcal{K}_{(j)}\,.$ (62)

Since all of these operators commute with one another for $\Omega=1$, their
conformal versions $\widetilde{{\cal K}}_{w}^{(j)}$, (51) also mutually
commute in the spacetime $\widetilde{{\boldsymbol{g}}}$. Of course, these are
nothing else than the operators $\widetilde{\mathcal{O}}_{(j)}$, (IV), this
time, however, written in a conformally invariant way (IV).444Although the
formulae (43) and (IV) look rather different, they represent the same
operators, and in particular, the coordinate expressions for the operators
$\widetilde{\mathcal{O}}_{(j)}$ and $\widetilde{{\cal K}}_{w}^{(j)}$ will
coincide in any conformal frame. The apparent differences arise from how we
choose scale the Killing tensors. The remaining commutation relations are then
guaranteed by (59), since we define for $j=0$

$\widetilde{{\cal
K}}_{w}^{(0)}\equiv\widetilde{\mathcal{O}}_{(0)}=\Omega^{2}\widetilde{\mathcal{K}}_{(0)}\,,$
(63)

reflecting the fact that the metric transforms differently than the other
Killing tensors under the conformal transformation.

Similarly one can ‘lift’ the operators $\mathcal{L}_{(j)}$, (23), to the
conformal ones (as in (50) and c.f. (44) where the Killing vectors transform
differently)

$\mathcal{L}^{(j)}_{w}=-i\,l_{(j)}^{a}\nabla_{\\!a}+i\frac{w}{D}(\nabla_{\\!a}\,l_{(j)}^{a})\;,$
(64)

where the second term identically vanishes in the frame $\Omega=1$ where
$l_{(j)}^{a}$ are (full, _not_ conformal,) Killing vectors. Of course, these
will coincide with $\widetilde{\mathcal{P}}_{(j)}$, (44), in any coordinate
system.

To summarize, we have found a conformally invariant ‘generalization’
$\\{\mathcal{K}^{(j)}_{w},\mathcal{L}^{(j)}_{w}\\}$ of the symmetry operators
(22) and (23), with the two being equal in the Kerr-NUT–AdS conformal frame
${\boldsymbol{g}}$. Writing $\widetilde{\Phi}=\Omega^{w}\Phi$ in any conformal
frame $\widetilde{{\boldsymbol{g}}}$, these operators obey the following
eigenvalue problem:

$\displaystyle\widetilde{\mathcal{K}}^{(j)}_{w}\widetilde{\Phi}$
$\displaystyle=$ $\displaystyle C_{j}\widetilde{\Phi}\,,$ (65)
$\displaystyle\widetilde{\mathcal{L}}^{(j)}_{w}\widetilde{\Phi}$
$\displaystyle=$ $\displaystyle\Psi_{j}\widetilde{\Phi}\,,$ (66)

guaranteeing $R$-separability of $\widetilde{\Phi}$ in any of these frames.

##  V Associated Hamilton–Jacobi equation and its separability

We finally turn to study the natural extension of the Hamiltonian–Jacobi
equation that arises from the the geometric optics (WKB) approximation of the
conformal wave equation.

Consider the following ‘$\alpha$-modified conformal wave equation’:

$\bigl{(}\alpha^{2}\Box-\eta R\bigr{)}\Phi=0\,.$ (67)

Then, upon employing the geometric optics ansatz

$\Phi=\Phi_{0}\exp\Bigl(\frac{i}{\alpha}S\Bigr{missing})\,,$ (68)

while taking the WKB limit $\alpha\to 0$, we arrive at the corresponding
Hamilton–Jacobi equation:

$g^{ab}\partial_{a}S\partial_{b}S+\eta R=0\,.$ (69)

This equation is obviously not conformally invariant, however, it is
consistent with the particle Hamiltonian,

$H=g^{ab}p_{a}p_{b}+\eta R\,.$ (70)

See e.g. DeWitt:1952js ; Omote:1976fx  for how such a coupling to the Ricci
scalar can arise from quantum corrections. The equations of motion for this
Hamiltonian yield the following modified geodesic equation

$\frac{Dp_{a}}{d\lambda}=-\eta\partial_{a}R\,.$ (71)

Let us stress that the procedure of deriving (69) is similar to how one
arrives at the massive Hamilton–Jacobi equation starting from the massive
($\alpha$-modified) Klein–Gordon one, e.g. Sergyeyev:2007gf . There is,
however, a fundamental difference. Namely, the $\alpha$-modified equation (67)
is not conformally invariant, unless $\alpha=1$. This is the reason why the
WKB limit $\alpha\to 0$ does not produce a conformally invariant
Hamilton–Jacobi equation. If instead, one started with the conformal wave
equation, setting $\alpha=1$ in (67), the WKB approximation would then yield
the massless Hamilton–Jacobi equation, which of course is conformally
invariant.

In what follows we consider the Hamilton–Jacobi equation (69) of potential
physical interest and show its separability in the off-shell Kerr–NUT–AdS
spacetimes. Using the form of the inverse metric given by (19) for $j=0$, the
Hamilton–Jacobi equation (69) takes the following explicit form:

$\displaystyle\sum_{\mu=1}^{n}\\!\left[\\!\frac{X_{\mu}}{U_{\mu}}{S_{\mu}^{\prime
2}}+\frac{1}{U_{\mu}X_{\mu}}\\!\left(\sum_{k=0}^{n-1+\varepsilon}(-x_{\mu}^{2})^{n-1-k}\\!\Psi_{k}\right)^{\\!\\!2}\right]$
$\displaystyle+\varepsilon\,\frac{1}{cA^{\\!(n)}}\,\Psi_{n}^{2}+\eta\sum_{\mu=1}^{n}\frac{r_{\mu}}{U_{\mu}}=0\;,$
(72)

where we have used the additive separation ansatz:

$S=\sum_{\mu=1}^{n}S_{\mu}(x_{\mu})+\sum_{k}\Psi_{k}\psi_{k}\,.$ (73)

Using next the following identity:

$\frac{1}{A^{(n)}}=\sum_{\mu}\frac{1}{x_{\mu}^{2}U_{\mu}}\,,$ (74)

we can rewrite the previous equation as

$\sum_{\mu}\frac{G_{\mu}}{U_{\mu}}=0\,,$ (75)

where

$G_{\mu}=X_{\mu}S_{\mu}^{\prime
2}+\frac{1}{X_{\mu}}\\!\left(\sum_{k=0}^{n-1+\varepsilon}(-x_{\mu}^{2})^{n-1-k}\,\Psi_{k}\right)^{\\!\\!2}+\varepsilon\,\frac{\Psi_{n}^{2}}{cx_{\mu}^{2}}+\eta
r_{\mu}\,.$ (76)

To proceed, we use the separation Lemma: The most general solution of

$\sum_{\mu=1}^{n}\frac{f_{\mu}(x_{\mu})}{U_{\mu}}=0\,,$ (77)

where $U_{\mu}$ is defined in (II), is given by

$f_{\mu}=\sum_{k=1}^{n-1}C_{k}(-x_{\mu}^{2})^{n-1-k}\,,$ (78)

where $C_{k}$ are arbitrary (separation) constants. This yields the following
ordinary differential equations for the separated solution:

$\displaystyle X_{\mu}S_{\mu}^{\prime
2}+\frac{1}{X_{\mu}}\\!\left(\sum_{k=0}^{n-1+\varepsilon}(-x_{\mu}^{2})^{n-1-k}\\!\Psi_{k}\right)^{\\!\\!2}+\varepsilon\,\frac{\Psi_{n}^{2}}{cx_{\mu}^{2}}+\eta
r_{\mu}$ $\displaystyle=\sum_{k=1}^{n-1}C_{k}(-x_{\mu}^{2})^{n-1-k}\,.$ (79)

Inverting this expression and identifying the canonical momenta
${\boldsymbol{p}}={\boldsymbol{d}}S$ the corresponding constants of motion of
the modified geodesic equation (71) are given by

$C_{j}=k^{(j)}_{ab}p^{a}p^{b}+\eta R_{(j)}\,,$ (80)

where $R_{(j)}$ are given by (III). It would be interesting to understand what
these constants of motion represent physically, e.g. in a quantum system
DeWitt:1952js ; Omote:1976fx , as this would give a natural interpretation for
the functions $R_{(j)}$.

##  VI Discussion

In this paper we have built on the previous work Gray:2020rtr  to find
covariant forms of the symmetry operators (22) and (23) of the conformal wave
equation in the Kerr–NUT–AdS background (11). These operators are built out of
the principal Killing–Yano tensor, its symmetry descendants, and the curvature
tensor. Moreover their commutativity descends naturally from the commutation
properties of the Killing tensors and the special character of the Ricci
scalar functions $R_{(j)}$, (III). We then showed how to lift these to a full
set of conformally invariant mutually commuting symmetry operators $\\{{\cal
K}_{w}^{(j)},{\cal L}_{w}^{(j)}\\}$ that guarantee $R$-separability of the
conformal wave equation in any conformally related spacetime
$\widetilde{{\boldsymbol{g}}}$, providing thus a highly non-trivial example to
the beautiful theory developed in Michel:2013dfa .

The conformal wave equation (1) is characterized by a specific value of
$\eta$. In principle one can consider more general wave equations, where
$\eta$ takes any value. It is easy to see that all such equations still
separate in the Kerr–NUT–AdS backgrounds; the operators (22) and (23) commute
for any value of $\eta$. However, for general $\eta$ the corresponding wave
equations are not conformally invariant and will not separate in a generic
conformally related spacetime. In this case, one could use the conformal
properties outlined in appendix B to construct an equation which separates in
the conformal spacetime, however there is no clear physical interpretation for
such an equation.

We have also introduced a modified Hamilton–Jacobi equation for a single
particle with a Ricci scalar potential term. This equation naturally arises
from the WKB limit of the ‘$\alpha$-modified’ conformal wave equation. This
limit breaks the conformal invariance and the resulting equation no longer
enjoys conformal symmetry. We have shown that this equation also separates in
the Kerr–NUT-AdS spacetimes – the corresponding non-trivial constants of
motion are given by the Killing tensors and the scalar functions $R_{(j)}$,
giving a natural setting for the interpretation of the latter.

In future, we would like to study the physical implications of the newly
derived (non-minimal coupling) Hamilton–Jacobi equation. We also hope to
extend the present results to understand separability of conformal fields with
higher spin.

## Acknowledgements

We would like to thank T. Bäckdahl for pointing to us the extended
mathematical literature on the conformal wave symmetry operators and in
particular the ref. Michel:2013dfa . F.G. acknowledges support from NSERC via
a Vanier Canada Graduate Scholarship. Y.Y. is supported by Grant-in-Aid for
Scientific Research from the Ministry of Education, Culture, Sports, Science
and Technology, Japan No.19K03877. This work was also supported by the
Perimeter Institute for Theoretical Physics and by the Natural Sciences and
Engineering Research Council of Canada (NSERC). Research at Perimeter
Institute is supported in part by the Government of Canada through the
Department of Innovation, Science and Economic Development Canada and by the
Province of Ontario through the Ministry of Colleges and Universities.
Perimeter Institute and the University of Waterloo are situated on the
Haldimand Tract, land that was promised to the Haudenosaunee of the Six
Nations of the Grand River, and is within the territory of the Neutral,
Anishnawbe, and Haudenosaunee peoples.

##  Appendix A Covariant form of $R_{(j)}$

In this appendix we find the covariant form of $R_{(j)}$ in Kerr–NUT–AdS
spacetimes by starting from the explicit expressions in canonical coordinates
(11). To start with we need an expression for the Ricci tensor. It is rather
simple since it is diagonal in the orthonormanl basis of the metric

$\displaystyle{\boldsymbol{e}}^{\mu}$
$\displaystyle=\frac{{\boldsymbol{d}}x_{\mu}}{\sqrt{Q_{\mu}}}\,,\quad\hat{{\boldsymbol{e}}}^{\mu}=\sqrt{Q_{\mu}}\sum_{k}A^{(k)}_{\mu}{\boldsymbol{d}}\psi_{k}\,,$
$\displaystyle{\boldsymbol{e}}^{0}$
$\displaystyle=\sqrt{\frac{c}{A^{\\!(n)}}}\sum_{k}A^{(k)}{\boldsymbol{d}}\psi_{k}\,,$
(81)

where $Q_{\mu}=X_{\mu}/U_{\mu}$. In fact it is given by Hamamoto:2006zf ;
Frolov:2017kze

$\text{\bf
Ric}=-\sum_{\mu}\hat{r}_{\mu}\left({\boldsymbol{e}}^{\mu}{\boldsymbol{e}}^{\mu}+\hat{{\boldsymbol{e}}}^{\mu}\hat{{\boldsymbol{e}}}^{\mu}\right)-\varepsilon\hat{r}_{0}{\boldsymbol{e}}^{0}{\boldsymbol{e}}^{0}\,,$
(82)

where we have introduced

$\displaystyle\hat{r}_{\mu}$
$\displaystyle=\frac{\hat{X}_{\mu}^{\prime\prime}+\frac{\epsilon\hat{X}_{\mu}^{\prime}}{x_{\mu}}}{2U_{\mu}}+\sum_{\nu\neq\mu}\frac{x_{\nu}\hat{X}_{\nu}^{\prime}-x_{\mu}\hat{X}_{\mu}^{\prime}-(1-\epsilon)(\hat{X}_{\nu}-\hat{X}_{\mu})}{\left(x_{\nu}^{2}-x_{\mu}^{2}\right)U_{\nu}}\;\;,\;\;$
$\displaystyle\hat{r}^{0}$
$\displaystyle=\sum_{\nu}\frac{\hat{X}^{\prime}_{\nu}}{x_{\nu}U_{\nu}}\;,\quad\hat{X}_{\mu}=X_{\mu}+\varepsilon
c/x_{\mu}^{2}\,.$ (83)

Also in this basis the Killing tensors are diagonal too,

${\boldsymbol{k}}_{(j)}=\\!\sum_{\mu=1}^{n}\\!A^{\\!(j)}_{\mu}\left[{\boldsymbol{e}}_{\mu}{\boldsymbol{e}}_{\mu}+\hat{{\boldsymbol{e}}}_{\mu}\hat{{\boldsymbol{e}}}_{\mu}\right]+\varepsilon
A^{\\!(j)}{\boldsymbol{e}}_{0}{\boldsymbol{e}}_{0}\,.$ (84)

Hence, using the identity

$\sum_{\nu\neq\mu}\frac{A^{\\!(j)}_{\mu}-A^{\\!(j)}_{\nu}}{x^{2}_{\nu}-x^{2}_{\mu}}=(n-j)A^{\\!(j-1)}_{\mu}\;,$
(85)

we have

$\displaystyle R_{(j)}-k^{ab}_{(j)}R_{ab}$
$\displaystyle=\sum_{\mu}\left[\varepsilon\frac{A^{\\!(j-1)}_{\mu}x_{\mu}\hat{X}_{\mu}^{\prime}}{U_{\mu}}+2\sum_{\nu\neq\mu}\frac{A^{\\!(j)}_{\mu}}{x^{2}_{\nu}-x^{2}_{\mu}}\left(\frac{x_{\mu}\hat{X}^{\prime}_{\mu}-(1-\varepsilon)\hat{X}_{\mu}}{U_{\mu}}+\frac{x_{\nu}\hat{X}^{\prime}_{\nu}-(1-\varepsilon)\hat{X}_{\nu}}{U_{\nu}}\right)\right]$
$\displaystyle=\sum_{\mu}\left[\varepsilon\frac{A^{\\!(j-1)}_{\mu}x_{\mu}\hat{X}_{\mu}^{\prime}}{U_{\mu}}+2\frac{x_{\mu}\hat{X}^{\prime}_{\mu}-(1-\varepsilon)\hat{X}_{\mu}}{U_{\mu}}\sum_{\nu\neq\mu}\frac{A^{\\!(j)}_{\mu}-A^{\\!(j)}_{\nu}}{x^{2}_{\nu}-x^{2}_{\mu}}\right]$
$\displaystyle=2\sum_{\mu}\frac{A^{\\!(j-1)}_{\mu}}{U_{\mu}}\left([n-j+\varepsilon/2]x_{\mu}\hat{X}^{\prime}_{\mu}-(n-j)(1-\varepsilon)\hat{X}_{\mu}\right)\;.$
(86)

Furthermore, since the Killing vectors satisfy
$\nabla\tensor{[}_{(a}]{l}{{}^{(j)}_{b)}}=0$ the only information of their
derivatives is contained in their exterior derivative. In particular, since in
the orthonormal basis (A)

${\boldsymbol{l}}^{(j)}=\sum_{\mu}A^{\\!(j)}_{\mu}\sqrt{Q_{\mu}}\hat{{\boldsymbol{e}}}^{\mu}+\epsilon
A^{\\!(j)}\sqrt{\frac{c}{A^{\\!(n)}}}{\boldsymbol{e}}^{0}\,,$ (87)

we have that

$\displaystyle{\boldsymbol{d}}{\boldsymbol{l}}^{(j)}=$
$\displaystyle\sum_{\mu}\left[\left(A^{\\!(j)}_{\mu}Q_{\mu}\frac{X_{\mu}^{\prime}}{X_{\mu}}-\varepsilon\frac{2}{x_{\mu}}\frac{cA^{\\!(j)}}{A^{\\!(n)}}+2x_{\mu}\sum_{\nu\neq\mu}\frac{Q_{\mu}A^{\\!(j)}_{\mu}+Q_{\nu}A^{\\!(j)}_{\nu}}{x^{2}_{\nu}-x^{2}_{\mu}}\right]{\boldsymbol{e}}^{\mu}\wedge\hat{{\boldsymbol{e}}}^{\mu}\right.$
$\displaystyle\left.+\varepsilon
2x_{\mu}\sqrt{Q_{\mu}}\sqrt{\frac{c}{A^{\\!(n)}}}A^{\\!(j-1)}_{\mu}{\boldsymbol{e}}^{\mu}\wedge{\boldsymbol{e}}^{0}+\sum_{\nu\neq\mu}2x_{\nu}\sqrt{Q_{\mu}Q_{\nu}}\frac{A^{\\!(j)}_{\mu}-A^{\\!(j)}_{\nu}}{x^{2}_{\nu}-x^{2}_{\mu}}{\boldsymbol{e}}^{\nu}\wedge\hat{{\boldsymbol{e}}}^{\mu}\right)\;.$
(88)

Moreover introducing the Killing co-potential

$(D-2j-1){\boldsymbol{\omega}}^{(j)}_{ab}:=k_{(j)\,a}^{\,\,\,n}h_{nb}=\sum_{\mu}A^{\\!(j)}_{\mu}x_{\mu}{\boldsymbol{e}}^{\mu}\wedge\hat{{\boldsymbol{e}}}^{\mu}\;.$
(89)

which generates the Killing tensors Frolov:2017kze

$l^{a}_{(j)}=\nabla_{b}\,\omega^{ba}_{(j)}\,,$ (90)

we can calculate

$k_{(j)n}^{\,\,\,a}h^{nb}\,dl_{ab}^{(k)}=2\sum_{\mu}\frac{1}{U_{\mu}}\left(A^{\\!(j)}_{\mu}A^{\\!(k)}_{\mu}x_{\mu}\hat{X}^{\prime}_{\mu}+\sum_{\nu\neq\mu}\frac{2\hat{X}_{\mu}(A^{\\!(j)}_{\mu}A^{\\!(k)}_{\mu}x^{2}_{\mu}-A^{\\!(j)}_{\nu}A^{\\!(k)}_{\nu}x^{2}_{\nu})-\varepsilon
cA^{\\!(j)}_{\mu}\left(\frac{A^{\\!(k)}_{\nu}}{U_{\nu}}-\frac{A^{\\!(k)}_{\mu}}{U_{\mu}}\right)}{x^{2}_{\nu}-x^{2}_{\mu}}\right)\;.$
(91)

Notice the last term proportional to $\varepsilon$ vanishes when $k=0$.
Finally let us us calculate $\Box\Tr({\boldsymbol{k}}_{(j)})$. First, we have

$\Tr({\boldsymbol{k}}_{(j)})=\varepsilon
A^{\\!(j)}+\sum_{\mu}2A^{\\!(j)}_{\mu}=(2(n-j)+\varepsilon)A^{\\!(j)}\;.$ (92)

Since this expression only depends on $x_{\mu}$ we can use the form of the
wave operator (see (20) in Gray:2020rtr ) to write

$\displaystyle\nabla_{a}(k_{(j)}^{ab}\nabla_{b}\Tr[{\boldsymbol{k}}_{(j)}])$
$\displaystyle=\sum_{\mu}\frac{A^{\\!(j)}_{\mu}}{U_{\mu}}\left[X_{\mu}\partial^{2}_{\mu}\Tr({\boldsymbol{k}}_{(j)})+\partial_{\mu}\Tr({\boldsymbol{k}}_{(j)})\left(X_{\mu}^{\prime}+\frac{\varepsilon}{x_{\mu}}X_{\mu}\right)\right]$
$\displaystyle=4\sum_{\mu}\frac{A^{\\!(j)}_{\mu}A^{\\!(j-1)}_{\mu}}{U_{\mu}}[n-j+\frac{\varepsilon}{2}]\left(x_{\mu}X^{\prime}_{\mu}+(1+\varepsilon)X_{\mu}\right)\,.$
(93)

Putting this together we have

$\displaystyle\alpha_{j}k_{(j-1)n}^{\,\,\,a}h^{nb}\,dl_{ab}^{(0)}-\beta_{j}l_{(j-1)}^{a}\,l_{a}^{(0)}+\frac{D-4}{2(D-2)}\Box\Tr({\boldsymbol{k}}_{(j)})$
$\displaystyle=2\sum_{\mu}\frac{1}{U_{\mu}}\left(A^{\\!(j-1)}_{\mu}\left(\left[\alpha_{j}+\frac{(D-4)(n-j+\frac{\varepsilon}{2})}{D-2}\right]x_{\mu}\hat{X}^{\prime}_{\mu}-\left[\frac{\beta_{j}}{2}-\frac{(D-4)(n-j+\frac{\varepsilon}{2})(1+\varepsilon)}{D-2}\right]\hat{X}_{\mu}\right)-2\alpha_{j}\hat{X}_{\mu}\sum_{\nu\neq\mu}\frac{A^{\\!(j)}_{\mu}-A^{\\!(j)}_{\nu}}{x^{2}_{\nu}-x^{2}_{\mu}}\right)$
$\displaystyle=2\sum_{\mu}\frac{A^{\\!(j-1)}_{\mu}}{U_{\mu}}\left(\left[\alpha_{j}+\frac{(D-4)(n-j+\frac{\varepsilon}{2})}{D-2}\right]x_{\mu}\hat{X}^{\prime}_{\mu}-\left[\frac{\beta_{j}}{2}+2(n-j)\alpha_{j}-\frac{(D-4)(n-j+\frac{\varepsilon}{2})(1+\varepsilon)}{D-2}\right]\hat{X}_{\mu}\right)\,.$
(94)

Thus, using $\varepsilon=\\{0,1\\}$ we can choose the coefficients to be

$\displaystyle\alpha_{j}$
$\displaystyle=\frac{2(n-j+\frac{\varepsilon}{2})}{D-2}\;,$ (95)
$\displaystyle\beta_{j}$
$\displaystyle=\frac{4(n-j+\frac{\varepsilon}{2})}{D-2}(D-3-2(n-j+\frac{\varepsilon}{2}))\,.$
(96)

Thence we obtain our covariant expression for $R_{(j)}$

$\displaystyle R_{(j)}=$ $\displaystyle
k^{ab}_{(j)}R_{ab}+\frac{D-4}{2(D-2)}\Box\Tr({\boldsymbol{k}}_{(j)})$
$\displaystyle+\alpha_{j}k_{(j-1)n}^{\,\,\,a}h^{nb}\,dl_{ab}^{(0)}-\beta_{j}l_{(j-1)}^{a}\,l_{a}^{(0)}\;,$
(97)

which matches the form in the text upon noting
${\boldsymbol{l}}_{0}={\boldsymbol{\xi}}$ and $D=2n+\varepsilon$.

Moreover the derivative of $R_{(j)}$ is particularly nice. We can calculate

$\displaystyle\nabla_{a}R_{(j)}\overset{a=\nu}{=}\sum_{\mu=1}^{n}\frac{\partial_{\nu}r_{\mu}\,A^{(j)}_{\mu}}{U_{\mu}}+2x_{\nu}A^{\\!(j)}_{\nu}\sum_{\mu\neq\nu}\frac{\frac{r_{\nu}}{U_{\nu}}+\frac{r_{\mu}}{U_{\mu}}}{x^{2}_{\mu}-x^{2}_{\nu}}$
$\displaystyle=\frac{r^{\prime}_{\nu}\,A^{(j)}_{\nu}}{U_{\nu}}+2x_{\nu}A^{\\!(j)}_{\nu}\sum_{\mu\neq\nu}\frac{\frac{r_{\nu}}{U_{\nu}}+\frac{r_{\mu}}{U_{\mu}}}{x^{2}_{\mu}-x^{2}_{\nu}}\;.$
(98)

Notice that one can also construct

$\displaystyle k_{(j)\,a}^{\;\;\;\;\;\;\;b}\nabla_{b}R\,\overset{a=\nu}{=}$
$\displaystyle\,\sum_{\mu}A^{\\!(j)}_{\nu}\,\partial_{\nu}\frac{r_{\nu}}{U_{\nu}}$
$\displaystyle=\frac{r^{\prime}_{\nu}\,A^{(j)}_{\nu}}{U_{\nu}}+2x_{\nu}A^{\\!(j)}_{\nu}\sum_{\mu\neq\nu}\frac{\frac{r_{\nu}}{U_{\nu}}+\frac{r_{\mu}}{U_{\mu}}}{x^{2}_{\mu}-x^{2}_{\nu}}$
$\displaystyle=\nabla_{a}R_{(j)}\;.$ (99)

Thus we have found a covariant expression for our symmetry operators’
derivatives

$\kappa_{a}^{(j)}:=k_{(j)\,a}^{\;\;\;\;\;\;\;b}\nabla_{b}R=\nabla_{a}R_{(j)}\;.$
(100)

Clearly ${\boldsymbol{\kappa}}$ is closed and also locally exact in all
dimensions thus we can say that our $R_{(j)}$ are the potentials for
${\boldsymbol{\kappa}}^{(j)}$ i.e.

${\boldsymbol{\kappa}}^{(j)}={\boldsymbol{d}}R_{(j)}\;.$ (101)

##  Appendix B Conformal Transformations

Given the spacetime $(\mathcal{M},{\boldsymbol{g}})$ we now consider a
conformal transformation of the metric, Killing tensors, and scalar field
(${\boldsymbol{k}}_{(j)}\rightarrow\Omega^{-2}{\boldsymbol{k}}_{(j)}$,
$\Phi\rightarrow\Omega^{w}\Phi$ for $w=1-D/2$) to the conformal spacetime
$(\mathcal{M},{\boldsymbol{g}},\Omega)$. The goal of this section is to find a
conformally covariant form of our wave operators

$\bigl{(}\hat{\mathcal{K}}_{(j)}-\eta
R_{(j)}\bigr{)}\Phi\,,\quad\hat{\mathcal{K}}_{(j)}=\nabla_{a}k^{ab}_{(j)}\nabla_{b}\,,\quad\eta=\frac{1}{4}\frac{D-2}{D-1}\,.$
(102)

Using the conformal properties of the Ricci tensor and covariant derivatives,
we find the following transformations

$\displaystyle\Omega^{2}\,\hat{\mathcal{K}}_{(j)}\Phi\rightarrow$
$\displaystyle\Omega^{w}\left(\hat{\mathcal{K}}_{(j)}+w\nabla_{a}(k^{ab}_{(j)}\nabla_{b}\log\Omega)\right.$
$\displaystyle\left.+w(w-2+D)\nabla_{a}\log\Omega\,k^{ab}_{(j)}\nabla_{b}\log\Omega\right)\Phi$
(103)

and

$\displaystyle\Omega^{2}\,k^{ab}_{(j)}R_{ab}\rightarrow$ $\displaystyle
k^{ab}_{(j)}R_{ab}-\left[(D-2)k^{ab}_{(j)}+k^{c}_{(j)\,c}g^{ab}\right]\nabla_{a}\nabla_{b}\log\Omega$
$\displaystyle+(D-2)\left[k^{ab}_{(j)}-k^{c}_{(j)\,c}g^{ab}\right]\nabla_{a}\log\Omega\nabla_{b}\log\Omega\;.$
(104)

Thence we have

$\displaystyle\Omega^{2}\left(\hat{\mathcal{K}}_{(j)}\Phi-\eta
k^{ab}_{(j)}R_{ab}\Phi\right)/\Phi\rightarrow$
$\displaystyle(\hat{\mathcal{K}}_{(j)}\Phi-\eta
k^{ab}_{(j)}R_{ab}\Phi)/\Phi+w(\nabla_{a}k^{ab}_{(j)})\nabla_{b}\log\Omega+((w+\eta(D-2))k^{ab}_{(j)}{+}\eta
k^{c}_{(j)c}g^{ab})\left[\nabla_{a}\nabla_{b}\log\Omega\right]$
$\displaystyle+(w(w-2+D)-(D-2)\eta)k^{ab}_{(j)}+(D-2)\eta
k^{c}_{(j)c}g^{ab})\left[\nabla_{a}\log\Omega\nabla_{b}\log\Omega\right]$
$\displaystyle=\left(\hat{\mathcal{K}}_{(j)}\Phi-\eta
k^{ab}_{(j)}R_{ab}\Phi\right)/\Phi+w(\nabla_{a}k^{ab}_{(j)})\nabla_{b}\log\Omega-\eta
D\,\hat{k}^{ab}_{(j)}\left[(D-2)\nabla_{a}\log\Omega\nabla_{b}\log\Omega+\nabla_{a}\nabla_{b}\log\Omega\right]\,.$
(105)

Here we have introduced the traceless Killing tensor
$\hat{k}^{ab}_{(j)}=k^{ab}_{(j)}-k^{c}_{(j)\,c}g^{ab}/D$. Clearly this
vanishes when $j=0$ so the first operator is conformally invariant. Notice
that the last term contains two derivatives of the conformal factor, so
consider the term identically zero term (following from the Killing tensor
equation)

$\displaystyle\nabla_{a}\nabla_{b}\left(k_{(j)}^{ab}+\frac{1}{2}k^{c}_{(j)\,c}g^{ab}\right)\equiv
0\,.$ (106)

Under the transformation
${\boldsymbol{k}}_{(j)}\rightarrow\Omega^{2}{\boldsymbol{k}}_{(j)}$ this
becomes

$\displaystyle\Omega^{2}\,\nabla_{a}\nabla_{b}\left(k_{(j)}^{ab}+\frac{1}{2}k^{c}_{(j)\,c}g^{ab}\right)$
$\displaystyle\rightarrow\nabla_{a}\nabla_{b}\left(k_{(j)}^{ab}+\frac{1}{2}k^{c}_{(j)\,c}g^{ab}\right)+(D+2)(\nabla_{a}k^{ab}_{(j)})\nabla_{b}\log\Omega$
$\displaystyle+D\,\hat{k}^{ab}_{(j)}\left[(D-2)\nabla_{a}\log\Omega\nabla_{b}\log\Omega+\nabla_{a}\nabla_{b}\log\Omega\right]\,.$
(107)

So we have

$\displaystyle\Omega^{2}\left(\hat{\mathcal{K}}_{(j)}\Phi-\eta\left[k^{ab}_{(j)}R_{ab}-\left\\{\nabla_{a}\nabla_{b}\left(k_{(j)}^{ab}+\frac{1}{2}k^{c}_{(j)\,c}g^{ab}\right)\right\\}\right]\Phi\right)/\Phi\rightarrow$
$\displaystyle\left(\hat{\mathcal{K}}_{(j)}\Phi-\eta\left[k^{ab}_{(j)}R_{ab}-\left\\{\nabla_{a}\nabla_{b}\left(k_{(j)}^{ab}+\frac{1}{2}k^{c}_{(j)\,c}g^{ab}\right)\right\\}+(D-4)(\nabla_{a}k^{ab}_{(j)})\nabla_{b}\log\Omega\right]\Phi\right)/\Phi\,.$
(108)

Note that, as the covariant derivatives and Killing tensors in the second line
are in the $\Omega=1$ frame, we have
$(D-4)(\nabla_{a}k^{ab}_{(j)})\nabla_{b}\log\Omega=-(D-4)/2\,(\nabla_{a}k^{c}_{(j)\,c})\nabla_{b}\log\Omega$.
Thus this term will be canceled by the transformation of
$\Box\Tr({\boldsymbol{k}}_{(j)})$. That is,

$\frac{D-4}{2(D-2)}\Box\Tr({\boldsymbol{k}}_{(j)})\rightarrow\Omega^{-2}\left(\frac{D-4}{2(D-2)}\Box\Tr({\boldsymbol{k}}_{(j)})+\frac{D-4}{2}\nabla_{a}\left[\Tr({\boldsymbol{k}}_{(j)})\right]\nabla^{a}\log\Omega\right)\,.$
(109)

We now consider the conformal transformation of the final piece;

${\cal
R}_{(j)}:=\alpha_{j}k_{(j-1)n}^{\,\,\,a}h^{nb}\,dl_{ab}^{(0)}-\beta_{j}l_{(j-1)}^{a}\,l_{a}^{(0)}\,.$
(110)

Now, if ${\boldsymbol{k}}_{(j)}\rightarrow\Omega^{-2}{\boldsymbol{k}}_{(j)}$
consistency demands that
${\boldsymbol{h}}\rightarrow\Omega^{2}{\boldsymbol{h}}$ and that
$l^{a}_{(j)}\rightarrow\Omega^{-2}l^{a}_{(j)}$. That is, one can show on a $p$
form $\star\rightarrow\Omega^{d-2p}\star$. Assuming $h\rightarrow\Omega^{r}h$;
$h^{j}\rightarrow\Omega^{jr}h^{j}$ then $f^{(j)}=\star
h^{j}\rightarrow\Omega^{d-4j+jr}f^{(j)}$. So

$k^{(j)}_{ab}\propto f_{ac_{1}\dots c_{D-2j-1}}f_{b}^{c_{1}\dots
c_{D-2j-1}}\rightarrow\Omega^{2(d-4j+jr)+2(D-2j-1)}k^{(j)}_{ab}=\Omega^{2+2j(-2+r)}k^{(j)}_{ab}\;.$
(111)

Hence demanding for all $j$ that
$k^{(j)}_{ab}\rightarrow\Omega^{2}k^{(j)}_{ab}$ fixes $r=2$. Then, we are left
with ${\cal R}_{(j)}$ as a scalar density of weight $-2$:

${\cal R}_{(j)}\rightarrow\Omega^{-2}\,{\cal R}_{(j)}\;.$ (112)

Thus putting this all together

$\Omega^{2}\,\left[\left(\hat{\mathcal{K}}_{(j)}-\eta\left[R_{(j)}-\left\\{\nabla_{a}\nabla_{b}\left(k_{(j)}^{ab}+\frac{1}{2}k^{c}_{(j)\,c}g^{ab}\right)\right\\}\right]\right)\Phi\right]/\Phi\rightarrow\\\
\left[\left(\hat{\mathcal{K}}_{(j)}-\eta\left[R_{(j)}-\left\\{\nabla_{a}\nabla_{b}\left(k_{(j)}^{ab}+\frac{1}{2}k^{c}_{(j)\,c}g^{ab}\right)\right\\}\right]\right)\Phi\right]/\Phi\,,$
(113)

which gives us the form we use in the main text.

## References

  * (1) V. Frolov, P. Krtous and D. Kubiznak, _Black holes, hidden symmetries, and complete integrability_ , _Living Rev. Rel._ 20 (2017) 6 [1705.05482].
  * (2) T. Houri, T. Oota and Y. Yasui, _Closed conformal Killing-Yano tensor and Kerr-NUT-de Sitter spacetime uniqueness_ , _Phys. Lett. B_ 656 (2007) 214 [0708.1368].
  * (3) P. Krtous, V. P. Frolov and D. Kubiznak, _Hidden Symmetries of Higher Dimensional Black Holes and Uniqueness of the Kerr-NUT-(A)dS spacetime_ , _Phys. Rev._ D78 (2008) 064022 [0804.4705].
  * (4) T. Houri, T. Oota and Y. Yasui, _Generalized Kerr-NUT-de Sitter metrics in all dimensions_ , _Phys. Lett. B_ 666 (2008) 391 [0805.0838].
  * (5) V. P. Frolov, P. Krtouš and D. Kubizňák, _Separability of Hamilton-Jacobi and Klein-Gordon equations in general Kerr-NUT-AdS spacetimes_ , _JHEP_ 0702 (2007) 005 [hep-th/0611245].
  * (6) T. Oota and Y. Yasui, _Separability of Dirac equation in higher dimensional Kerr-NUT-de Sitter spacetime_ , _Phys. Lett._ B659 (2008) 688 [0711.0078].
  * (7) O. Lunin, _Maxwell’s equations in the Myers-Perry geometry_ , _JHEP_ 12 (2017) 138 [1708.06766].
  * (8) V. P. Frolov, P. Krtouš, D. Kubizňák and J. E. Santos, _Massive Vector Fields in Rotating Black-Hole Spacetimes: Separability and Quasinormal Modes_ , _Phys. Rev. Lett._ 120 (2018) 231103 [1804.00030].
  * (9) O. Lunin, _Excitations of the Myers-Perry Black Holes_ , _JHEP_ 10 (2019) 030 [1907.03820].
  * (10) F. Gray, I. Holst, D. Kubiznak, G. Odak, D. M. Pirvu and T. R. Perche, _Conformally Coupled Scalar in Rotating Black Hole Spacetimes_ , _Phys. Rev. D_ 101 (2020) 084031 [2002.05221].
  * (11) R. M. Wald, _General Relativity_. Chicago Univ. Pr., Chicago, USA, 1984.
  * (12) J. Hennig and R. Panosso Macedo, _Fully pseudospectral solution of the conformally invariant wave equation on a Kerr background_ , 2012.02240.
  * (13) J.-P. Michel, F. Radoux and J. Šilhan, _Second order symmetries of the conformal Laplacian_ , _SIGMA_ 10 (2014) 016 [1308.1046].
  * (14) B. S. DeWitt, _Point transformations in quantum mechanics_ , _Phys. Rev._ 85 (1952) 653.
  * (15) M. Omote, _Point Canonical Transformations and the Path Integral_ , _Nucl. Phys. B_ 120 (1977) 325.
  * (16) C. Destri, P. Maraner and E. Onofri, _On the definition of quantum free particle on curved manifolds_ , _Nuovo Cim. A_ 107 (1994) 237 [hep-th/9210027].
  * (17) D. Lian, L. Hu and Q. Liu, _Identification of Geometric Potential from Quantum Conditions for a Particle on a Curved Surface_ , 1701.08370.
  * (18) N. Hamamoto, T. Houri, T. Oota and Y. Yasui, _Kerr-NUT-de Sitter curvature in all dimensions_ , _J. Phys. A_ 40 (2007) F177 [hep-th/0611285].
  * (19) A. Sergyeyev and P. Krtous, _Complete Set of Commuting Symmetry Operators for Klein-Gordon Equation in Generalized Higher-Dimensional Kerr-NUT-(A)dS Spacetimes_ , _Phys. Rev. D_ 77 (2008) 044033 [0711.4623].
  * (20) S. Benenti, C. Chanu and G. Rastelli, _Remarks on the connection between the additive separation of the hamilton–jacobi equation and the multiplicative separation of the schrödinger equation. i. the completeness and robertson conditions_ , _Journal of Mathematical Physics_ 43 (2002) 5183.
  * (21) D. Kubiznak and P. Krtous, _On conformal Killing-Yano tensors for Plebanski-Demianski family of solutions_ , _Phys. Rev. D_ 76 (2007) 084036 [0707.0409].
  * (22) B. Carter, _Killing tensor quantum numbers and conserved currents in curved space_ , _Physical Review D_ 16 (1977) 3395.
  * (23) C. Boyer, E. Kalnins and W. Miller, _Symmetry and separation of variables for the helmholtz and laplace equations_ , _Nagoya Mathematical Journal_ 60 (1976) 35.
  * (24) E. G. Kalnins and W. Miller Jr, _Intrinsic characterisation of orthogonal r separation for laplace equations_ , _Journal of Physics A: Mathematical and General_ 15 (1982) 2699.
  * (25) N. Kamran and R. McLenaghan, _Separation of variables and symmetry operators for the conformally invariant klein-gordon equation on curved spacetime_ , _letters in mathematical physics_ 9 (1985) 65.
  * (26) C. Duval, P. Lecomte and V. Ovsienko, _Conformally equivariant quantization: existence and uniqueness_ , in _Annales de l’institut Fourier_ , vol. 49, 1999.
  * (27) M. Eastwood, _Higher symmetries of the laplacian_ , _Annals of mathematics_ (2005) 1645.
  * (28) M. Eastwood and T. Leistner, _Higher symmetries of the square of the laplacian_ , in _Symmetries and overdetermined systems of partial differential equations_ , pp. 319–338.  Springer, 2008.
  * (29) A. R. Gover and J. Šilhan, _Higher symmetries of the conformal powers of the laplacian on conformally flat manifolds_ , _Journal of mathematical physics_ 53 (2012) 032301.
  * (30) L. Andersson, T. Bäckdahl and P. Blue, _Second order symmetry operators_ , _Classical and Quantum Gravity_ 31 (2014) 135015\.


# Expert Selection in Distributed Gaussian Processes: A Multi-label
Classification Approach

Hamed Jalali
Department of Computer Science
University of Tuebingen
<EMAIL_ADDRESS>
&Gjergji Kasneci
Department of Computer Science
University of Tuebingen
<EMAIL_ADDRESS>

###### Abstract

By distributing the training process, local approximation reduces the cost of
the standard Gaussian Process. An ensemble technique combines local
predictions from Gaussian experts trained on different partitions of the data
by assuming a perfect diversity of local predictors. Although it keeps the
aggregation tractable, this assumption is often violated in practice. Taking
dependencies between experts enables ensemble methods to provide consistent
results. However, they have a high computational cost, which is cubic in the
number of experts involved. By implementing an expert selection strategy, the
final aggregation step uses fewer experts and is more efficient. Indeed, a
static selection approach that assigns a fixed set of experts to each new data
point cannot encode the specific properties of each unique data point. This
paper proposes a flexible expert selection approach based on the
characteristics of entry data points. To this end, we investigate the
selection task as a multi-label classification problem where the experts
define labels, and each entry point is assigned to some experts. The proposed
solution’s prediction quality, efficiency, and asymptotic properties are
discussed in detail. We demonstrate the efficacy of our method through
extensive numerical experiments using synthetic and real-world data sets.

##  1 Introduction and Related Works

Gaussian processes (GPs) [1] are interpretable and powerful Bayesian non-
parametric methods for non-linear regression. By applying Bayes’ theorem for
inference, the posterior predictive distribution of a GP is the best linear
unbiased estimator (BLUE) under the assumed model and provides proper
quantification of the prediction error uncertainty. GPs do not need
restrictive assumptions of the model and can estimate complex linear and non-
linear structures. While GPs are extensively used in practical cases [2, 3, 4,
5, 6, 7], their cubic training and quadratic prediction costs 111I.e., in the
size of the training set. limit their application to big data use cases [8].

To reduce computational expense, the prominent distributed Gaussian processes
(DGPs) divide the training data into numerous subsets. In this scenario, the
local inference is conducted for each partition separately, and then the local
estimations are combined through ensemble learning [9, 10, 11]. A local GP
with expertise in a particular section is called an expert. Experts share the
same hyper-parameters, thus accounting for implicit regularisation and
encountering overfitting [8, 12]. In a DGP, the _conditional independence_
(CI) assumption between partitions (i.e., between experts given the target)
allows for factorizing the global posterior distribution as a product of local
distributions. While this assumption reduces the computational cost, it
results in inconsistencies and suboptimal solutions [13]. Relaxing the
independence assumption raises the aggregation’s theoretical properties. If
the experts’ predictions are assumed to be random variables, their relative
correlations define dependencies between experts. The aggregated posterior
distribution, in this case, provides high-quality forecasts and is capable of
returning consistent results [14]. However, solutions that deal with the
consistency problem suffer from extra computational costs induced by the need
to find the inverse of the covariance matrix between experts for each test
point [12].

To boost the efficiency of dependency-based aggregation, authors in [15, 16]
discuss complexity reduction as an expert selection scenario that excludes a
subset of original experts and considers only the valuable experts in the
aggregation. For this purpose, they use the Gaussian graphical model (GGM),
where the nodes are original experts, and the edges are their interactions.
The nodes with fewer interactions are defined as unimportant experts and
excluded from the model. This approach can lower the complexity and provide a
good approximation for the original estimator. However, it is not flexible
concerning new entries, and the selected experts are fixed for all test
points.

The critical contribution of our work lies in selecting a dynamic and entry-
dependent subset of local experts for each new data point using multi-label
classification (MLC). We propose a flexible mechanism for each new observation
that designates related experts to provide local predictions. Multi-label
classification [17] is a generalization of multi-class classification, where
multiple labels may be assigned to each instance. To transform the distributed
learning case into a multi-label classification, the indices of the
partitions/experts are the labels/classes. The task is to assign some experts
to a new data point. Multi-class classification problem that allocates only
one expert per test point produces discontinuous separation boundaries between
sub-regions and therefore is not a proper solution for quantifying
uncertainties [18, 12].

In this work, we adopt two MLC algorithms to assign experts to data points:
k-nearest neighbors (KNN) and deep neural networks (DNN). For the first one,
we use the partition’s centroid as a substitute for the corresponding local
expert. By estimating the distance between a new entry point and the
centroids, we can find its $K$ nearest neighboring experts. For the second
approach, we train a neural network with a soft-max output layer and log-loss
(i.e., cross-entropy loss) using the train points and their related partition
indices that shows the partition they belong to. After training the DNN, we
send a new test point through the network, and the experts with higher
probability are assigned to this test point. Relative to consistent
aggregation methods that use dependency information, our approach keeps all
asymptotic properties of the original baseline and substantially provides
competitive prediction performance while leading to better computational costs
than other SOTA approaches, which use the dependency assumption.

The structure of the paper is as follows. Section 2 introduces the problem
formulation and related work. In Section 3, the proposed model and the
inference process are presented. Section 4 shows the experimental results, and
we conclude in Section 5.

##  2 Problem Set-up

Background: We consider the non-linear regression problem $y=f(x)+\epsilon$,
where $x\in R^{d}$ and $\epsilon\sim\mathcal{N}(0,\sigma^{2})$, and the
Gaussian likelihood is $p(y|f)=\mathcal{N}(f,\sigma^{2}I)$. The objective is
to learn the latent function f from a training set $\mathcal{D}=\\{X,y\\}$ of
size $n$. The GP then describes a prior distribution over the latent functions
as $f\sim GP\left(0,k(x,x^{\prime})\right)$, where $k(x,x^{\prime})$ is the
covariate function (kernel) with hyperparameters $\psi$, and $x,x^{\prime}\in
X$. The training step of GPs scales as $\mathcal{O}(n^{3})$ because it is
affected by the inversion and determinant of the $n\times n$ matrix
$\mathcal{C}$. It imposes limitations on the scalability of GPs for large data
sets. In distributed GP (DGP), one divides the full training data set into $M$
partitions (called experts)
$\mathcal{D}^{{}^{\prime}}=\\{\mathcal{D}_{1},\ldots,\mathcal{D}_{M}\\}$ and
trains standard GPs on these partitions. All GP experts are trained jointly
and share a single set of hyper-parameters $\theta=\\{\sigma^{2},\psi\\}$. For
a test set $X^{*}$ of size $n_{t}$, the predictive distribution of the $i$-th
expert $\mathcal{M}_{i}$ is
$p_{i}(y^{*}|\mathcal{D}_{i},X^{*})\sim\mathcal{N}(\mu_{i}^{*},\Sigma_{i}^{*})$.

Aggregation: Aggregating the experts in a DGP is typically based on the
assumption that they are conditionally independent (CI). Therefore, for a test
input $x^{*}$, the posterior distribution of DGP is given as:

$p(y^{*}|\mathcal{D},x^{*})\propto\prod_{i=1}^{M}p_{i}^{\beta_{i}}(y^{*}|\mathcal{D}_{i},x^{*}),$
(1)

where the weights $\beta=\\{\beta_{1},\ldots,\beta_{M}\\}$ describe the
importance and influence of the experts. The most popular aggregation methods
are the product of experts [19, 20] and Bayesian committee machine [21, 8,
22]. However, their predictions are not accurate enough, and CI-based
aggregation generally returns a sub-optimal solution [23].

To consider dependencies between experts in the aggregation step, assume that
the means of local predictive distributions $\mu^{*}_{i}$ are random variables
and $\mu^{*}(x^{*})=[\mu_{1}^{*}(x^{*}),\ldots,\mu_{M}^{*}(x^{*})]^{T}$ be an
$M\times 1$ vector that contains the centered predictions of $M$ experts for a
given test point $x^{*}\in X^{*}$. This allows us to consider correlations
between the experts’ predictions and latent variable $y^{*}$,
$k_{A}(x^{*})=Cov\left(\mu^{*}(x^{*}),y^{*}(x^{*})\right)$, and also leverage
internal correlations between experts,
$K_{A}(x^{*})=Cov\left(\mu^{*}(x^{*}),\mu^{*}(x^{*})\right)$. Thus, for the
test point, $x^{*}$ and sub-model predictions $\mu^{*}(x^{*})$, the aggregated
predictor is defined as

$y_{A}^{*}(x^{*})=k_{A}(x^{*})^{T}K_{A}(x^{*})^{-1}\mu^{*}(x^{*}).$ (2)

In [14], the authors showed that this linear estimator is the best linear
unbiased predictor (BLUP). This method is known as the nested pointwise
aggregation of experts (NPAE) and provides high-quality predictions. However,
it needs the inverse of the $M\times M$ auto-covariance matrix $K_{A}(x^{*})$,
and its time complexity is cubic in the number of experts, i.e.,
$\mathcal{O}(n_{t}M^{3})$, and it is not an efficient solution for large real-
world data sets.

##  3 Expert Selection in Local Approximation GPs

Despite the high aggregation accuracy in 2, its computational cost is a
challenge for using the method on large data sets. Expert selections can
improve the performance of dependency-based aggregation in two ways. First,
unrelated experts for a given data point can be excluded, and only informative
partitions can be considered to make the prediction. Second, mitigating the
number of experts reduces the prediction cost and enables the ensemble to be
used in large data sets.

###  3.1 Expert Selection Using Graphical Models

Gaussian graphical model (GGM) is the continuous form of pairwise Markov
random fields. It assumes the nodes of an undirected graph are random
variables, and the joint distribution of the random variables is multivariate
Gaussian distribution with zero mean and precision matrix $\Omega$,
$\mathcal{N}(0,\Omega^{-1})$. The elements of the precision matrix are the
unknown parameters and show interactions between experts (edges in the graph).
In the DGP framework, the experts’ predictions $\mu^{*}$ are the nodes in a
GGM. Authors in [15] proposed an expert selection method that uses only the
highly connected nodes in the graph as the most critical experts in final
aggregation defined in Equation 2. Although this method can mitigate the
computational cost, it selects a fixed and static set of experts for all test
points and, therefore, can not capture the specific properties of the new data
points.

###  3.2 Multi-label Classification for Flexible Expert Selection

Assigning experts to new entry points in a distributed learning model can be a
classification problem. Let’s assume that each expert is a class of
estimators. The selection problem for each test point $x^{*}$ is defined as a
multi-label classification task where each instance can be associated with
some classes. The main advantage of this method is its flexibility because the
selected experts depend on the given test point, and thus different experts
can be assigned to different test points. Assume $x^{*}$ is a new test point
and $\mathcal{M}=\\{\mathcal{M}_{1},\ldots,\mathcal{M}_{M}\\}$ is the Gaussian
experts set, and $\mathcal{L}=\\{1,\ldots,M\\}$ is the label set. The task is
to find
$\mathcal{M^{C}}(x^{*})=\\{\mathcal{M^{C}}_{1}(x^{*}),\ldots,\mathcal{M^{C}}_{K}(x^{*})\\}$
that represents $K$ selected experts to predict at $x^{*}$. We adopt two
prominent classification models to solve this multi-label task without
requiring problem transformations, K-nearest neighbors (KNN) and conventional
deep neural networks (DNN).

KNN: Let
$\mathcal{D}^{{}^{\prime}}=\\{\mathcal{D}_{1},\ldots,\mathcal{D}_{M}\\}$ be
the partitions based on a disjoint partitioning strategy, i.e. K-Means
clustering. Also, assume
$\mathcal{C}=\\{\mathcal{C}_{1},\ldots,\mathcal{C}_{M}\\}$ contains the
related centroids of the clusters in $\mathcal{D}$. For each test point
$x^{*}$, there is a $1\times M$ vector $dist(x^{*},\mathcal{C})$ , in which
the $i$’th element is the distance between $x^{*}$ and $\mathcal{D}_{i}$,
where $dist()$ is a distance metric. The adopted KNN selects $K$ experts with
the closest centroids to $x^{*}$.

DNN: For each $x_{i}\in X$, $i=1,\ldots,N$, the related partition label
$l_{i}\in\mathcal{L}$ is available as an output of the training step in DGP.
Therefore, instead of the original training set $(X,y)$, a new set of points
and labels $(X,\mathcal{L})$ is constructed for training the DNN. After that,
for each test point, $x^{*}\in X^{*}$, the network will provide a probability
vector $P^{\mathcal{L}}(x^{*})$ where $P^{\mathcal{L}}(x^{*})_{j}$ represents
the probability that $x^{*}$ belongs to the $j$’th expert. The $K$ partitions
with highest probabilities in $P^{\mathcal{L}}(x^{*})$ are assigned to
$x^{*}$.

##  4 Experiments

To evaluate the prediction quality of static and dynamic expert selection
methods, we use Concrete Compressive Strength
222https://archive.ics.uci.edu/ml/datasets/concrete+compressive+strength data
set. It is a small data set containing 1030 observations of 9 attributes (8
independent variables and one response variable). We use $\%90$ of the
observations for training and the rest for testing, where disjoint
partitioning is used to divide the data set into ten subsets. The quality of
predictions is evaluated in two ways, standardized mean squared error (SMSE)
and the mean standardized log loss (MSLL) [1]. Figure 1 depicts the prediction
quality of expert selection methods on the Concrete data set for different
values of $K$. The figure shows that multi-label-based expert selection models
provide higher quality predictions with lower deviation from the NPAE model.
The quality of the classification-based aggregations changes less as the
selection parameter $K$ increases. The case $K=10$ leads to the original NPAE
baseline (the dashed green lines), and both KNN and DNN return proper error
values in both plots.

(a)

(b)

Figure 1: Expert Selection prediction qualities of different experts selection
methods compared to original baseline, NPAE, from Concrete data set.

##  5 Conclusion

In this work, we have proposed a novel expert selection approach for
distributed learning with Gaussian agents, which leverages expert selection to
aggregate dependent local experts’ predictions. The available ensemble
baselines use all correlated experts in the aggregation step. It affects the
final predictions by the local predictions of weak experts or leads to
impractically high computational costs. Our proposed approach uses a multi-
label classification model, and the allocation of data points to experts is
defined by considering the experts as class labels. Unlike the available
expert selection method, the proposed model is more flexible to the model and
data changes and chooses a related group of experts at each entry point.
Excluding unrelated experts at each test point improves the prediction quality
and reduces computational costs. The classification methods in this work,
i.e., KNN and DNN, can be replaced with recent and more efficient solutions
proposed to solve the multi-label classification problem. The proposed
approach can be used for distributed and federated learning and does not
impose restricted assumptions.

## References

  * [1] C. E. Rasmussen and C. K. Williams.  Gaussian processes for machine learning.  MIT Press, 2006.
  * [2] Volker L. Deringer, Albert P. Bartók, Noam Bernstein, David M. Wilkins, Michele Ceriotti, and Gábor Csányi.  Gaussian process regression for materials and molecules.  Chemical Reviews, 121(16):10073–10141, 2021.
  * [3] Aaron Zeng, Hodde Ho, and Yao Yu.  Prediction of building electricity usage using gaussian process regression.  Journal of Building Engineering, 28(11):101054, 2020.
  * [4] A. Denzel and J. Kästner.  Gaussian process regression for geometry optimization.  he Journal of Chemical Physics, 148(9):094114, 2018.
  * [5] R.B. Gramacy.  urrogates: Gaussian Process Modeling, Design, and Optimization for the Applied Sciences.  Chapman and Hall/CRC, 2020.
  * [6] T. Le, K. Nguyen, V. Nguyen, T. D. Nguyen, and D. Phung.  Gogp: Fast online regression with gaussian processes.  IEEE International Conference on Data Mining, pages 257–266, 2017\.
  * [7] F. Tobar, T. D. Bui, and R. E. Turner.  Learning stationary time series using gaussian processes with nonparametric kernels.  In Advances in Neural Information Processing Systems, pages 3501–3509, 2015.
  * [8] M. P. Deisenroth and J. W. Ng.  Distributed gaussian processes.  International Conference on Machine Learning, pages 1481–1490, 2015\.
  * [9] D. Peteiro-Barral and B. Guijarro-Berdiñas.  A survey of methods for distributed machine learning.  Progress in Artificial Intelligence, 2(1):1–11, 2013.
  * [10] J. Zhang, M. Wang, Q. Li, S. Wang, X. Chang, and B. Wang.  A survey on distributed machine learning.  ACM Computing Surveys, 53(2):1–33, 2020.
  * [11] J. Liu, J. Huang, Y. Zhou, X. Li, S. Ji, H. Xiong, and D. Dou.  From distributed machine learning to federated learning: A survey.  arXiv preprint arXiv:2104.14362, 2021.
  * [12] X. Shen H. Liu, Y. S. Ong and J. Cai.  When gaussian process meets big data: A review of scalable gps.  IEEE Transactions on Neural Networks and Learning Systems, 31(11):4405–4423, 2020.
  * [13] B. Szabó and H. van Zanten.  An asymptotic analysis of distributed nonparametric methods.  Journal of Machine Learning Research, 20(87):1–30, 2019.
  * [14] D. Rullière, N. Durrande, F. Bachoc, and C. Chevalier.  Nested kriging predictions for datasets with a large number of observations.  Statistics and Computing, 28(4):849–867, 2018.
  * [15] H. Jalali, M. Pawelczyk, and G. Kasneci.  Model selection in local approximation gaussian processes: A markov random fields approach.  IEEE International Conference on Big Data (Big Data), pages 768–778, 2021.
  * [16] H. Jalali, M. Pawelczyk, and G. Kasneci.  Gaussian experts selection using graphical models.  arXiv preprint arXiv:2102.01496, 2021.
  * [17] Francisco Herrera, Francisco Charte, Antonio J. Rivera, and María J. del Jesus.  Multilabel Classification, pages 17–31.  Springer International Publishing, 2016.
  * [18] C. Park and D. Apley.  Patchwork kriging for large-scale gaussian process regression.  Journal of Machine Learning Research, 19(7):1–43, 2018.
  * [19] G. E. Hinton.  Training products of experts by minimizing contrastive divergence.  Neural Computation, 14(8):1771–1800, 2002.
  * [20] Y. Cao and D. J. Fleet.  Generalized product of experts for automatic and principled fusion of gaussian process predictions.  arXiv preprint arXiv:1410.7827, 2014.
  * [21] V. Tresp.  A bayesian committee machine.  Neural Computation, 12(11):2719–2741, 2000.
  * [22] H. Liu, J. Cai, Y. Ong, and Y. Wang.  Generalized robust bayesian committee machine for large-scale gaussian process regression.  International Conference on Machine Learning, pages 1–10, 2018\.
  * [23] A. Jaffe, E. Fetaya, B. Nadler, T. Jiang, and Y. Kluger.  Unsupervised ensemble learning with dependent classifiers.  In Artificial Intelligence and Statistics, pages 351–360, 2016\.