Superintegrable and scale invariant quantum mechanical systems with position
dependent mass

A. G. Nikitin111E-mail<EMAIL_ADDRESS>
Institute of Mathematics of National Academy of Sciences of Ukraine,
3 Tereshchenkivska Street, Kyiv 4, Ukraine, 01004

Abstract

Schrödinger equations with position dependent mass which are scale invariant
and admit second order integrals of motion are classified.

##  1 Introduction

It is well known that fundamental equations of mathematical physics admit
rather exten-ded symmetries. At the first place we can mention symmetries with
respect to continuous groups (Lie symmetries). In addition, there exists an
extended class of generalised symmet-ries (higher symmetries, supersymmetries,
hidden symmetries) connected with the menti-oned equations, see, e.g., book
[1] and paper [2].

An important class of equations of quantum mechanics (QM) is formed by the
superin-tegrable Schrödinger equations which admit more integrals of motion
than the number of degrees of freedom of the associated QM systems. As a rule
such equations also admit a wide Lie symmetry and different kinds of
generalised symmetries. A known example is the Schrödinger equations for the
Hydrogen atom which, in addition to its invariance with respect to the
rotation group, is supersymmetric and admits the hidden (Fock) symmetry with
respect to group O(4).

Searching for Lie symmetries of partial differential equations is a rather
popular business which generates a big number of publications. Rather
surprisingly, the correct group classification of Schrödinger equations for a
particle interacting with an external field appears only recently [3, 4, 5,
6], in spite of that the search for the related symmetries started long time
ago, see papers [7, 8, 9].

The systematic search for superintegrable QM systems was started in papers of
Yakov Smorodinsky and his collaborators [10], [11], who completely described
all inequivalent 2d Schrödinger equations which admit second order integrals
of motion. The modern trend is the studying of the related integrals of motion
of the third and even arbitrary orders [12]. See also [13] where the
determining equations for such symmetries were deduced.

The present paper is devoted to the classification of superintegrable
Schrödinger equa-tions with position dependent mass (PDM). Such equations have
very important applica-tions in the modern theoretical physics, being the
mathematical models of semiconductors [14], quantum liquids [15], quantum dots
[16], and many, many others quantum objects.

Higher symmetries of the 2d PDM Schrödinger equations are well known (see,
e.g., [17] and the references given therein). However, there are only some
particular results concerning the classification of such equations with three
independent variables [18, 19, 20, 21]. In the following we present the
complete classification of the special class of such equations, namely, of
equations which are scale invariant.

For the group classification of the stationary and time dependent PDM
Schrödinger equations and the related reaction-diffusion systems see papers
[22, 23] and [24, 25] corres-pondingly. In accordance with these papers such
equations cannot be Galilei invariant in contrast with the standard
Schrödinger equation and its generalisations for the case of particles with
arbitrary spin [26].

The classification of the higher symmetries of the 3d PDM Schrödinger
equations is a very complicated problem which is solved only for some special
classes of such equations. An important and rather extended class of them
includes equations which admit at least one parametric symmetry group. All
inequivalent equations from this class which admit first order integrals of
motion is fixed in [22].

In the present paper we classify the scale invariant 3d PDM Schrödinger
equations which admit at least one second order integral of motion and show
that they include many interesting and consistent equations. The found
equations form the important part of a more general class of superintegrable
equations which admit at least a minimal Lie symmetry.

##  2 PDM Schrödinger equations

We will study the stationary PDM Schrödinger equations of the following
generic form:

$\displaystyle\hat{H}\psi=E\psi,$ (1)

where

$\displaystyle\hat{H}=p_{a}f({\bf x})p_{a}+V({\bf x}).$ (2)

In equation (1) ${\bf x}=(x^{1},x^{2},x^{3}),$ $p_{a}=-i\partial_{a}$, а
$V({\bf x})$ та $f({\bf x})=\frac{1}{2m({\bf x})}$ are arbitrary functions of
${\bf x}$ which are associated with the potential and inverse mass
correspondingly, and summation is imposed over the repeated indices.

In paper [22] all inequivalent equations (1) admitting at least one first
order integral of motion (i.e., differential operator of first order commuting
with Hamiltonian $H$) are found. In [21] the complete classification of the
special class of superintegrable equations (1) was carried out, namely, the
equations invariant with respect to the rotation group.

In the present paper we classify superintegrable equations (1) which are scale
invariant. The corresponding arbitrary elements $V$ and $f$ take the form
[22]:

$\displaystyle f=r^{2}F(\varphi,\theta),\quad V=V(\varphi),$ (3)

where $F(.)$ and $V(.)$ are arbitrary functions, $\varphi$ і $\theta$ are the
Euler angles.

Hamiltonians (2) with special arbitrary elements (3) commute with the
generator of the scaling transformations which has the following form:

$\displaystyle D=x_{a}p_{a}-\frac{3{\mathrm{i}}}{2}.$

Our task is to find such of them which commute with at least one first order
differential operator.

##  3 Equivalence group

Changes of dependent and independent variables are called the equivalence
transformati-ons provided they keep the generic form of the differential
equation (in our case of equation (1)) up to the changes of the arbitrary
elements (in our case functions $f$ and $V$). The set of the equivalence
transformations has the gruppoid structure [27], and can include equivalence
groups and some discrete elements.

As it is shown in [22], the maximal continuous equivalence group of equation
(1) is C(3), i.e., the group of conformal transformations of 3d Euclidean
space. The generators of this group are the following first order differential
operators:

$\displaystyle\begin{split}&P^{a}=p^{a}=-i\frac{\partial}{\partial
x_{a}},\quad L^{a}=\varepsilon^{abc}x^{b}p^{c},\\\
&D=x_{n}p^{n}-\frac{3{\mathrm{i}}}{2},\quad
K^{a}=r^{2}p^{a}-2x^{a}D,\end{split}$ (4)

where $r^{2}=x_{1}^{2}+x_{2}^{2}+x_{3}^{2}$ and
$p_{a}=-i\frac{\partial}{\partial x_{a}}.$ The corresponding group
transformations (whose explicit form can be found in [22]) keep the generic
form of equations (1), (2) up to exact form of functions $f$ and $V$. The
important particular form of these transformations is the inversion:

$\displaystyle x_{a}\to\tilde{x}_{a}=\frac{x_{a}}{x^{2}},\quad\psi({\bf
x})\to\tilde{x}^{3}\psi(\tilde{\bf x})$ (5)

which acts on operators (4) in the following manner:

$\displaystyle P_{a}\to K_{a},\quad K_{a}\to P_{a},\quad L_{a}\to L_{a},\quad
D\to D.$ (6)

For the class of equations considered in the present paper the equivalence
group is reduced to the direct product of the rotations group whose generators
are componets of the orbital momenta $L_{a}$ and dilatation transformations
generated by $D$, since $L_{a}$ and $D$ commute with $H$ while the remaining
operators (4) do not have this property. Notice that the discrete equivalence
transformation (5) is kept also.

In the following we will use the rotations and the inverse transformation (5)
for optimisation of the requested calculation and selection of non-equivalent
versions of the studied equations.

##  4 Determining equations

The searched second order integrals of motion can be represented in the
following form:

$Q=\mu^{ab}\partial_{a}\partial_{b}+\xi^{a}\partial_{a}+\eta$ (7)

where $\mu^{ab}=\mu^{ba}$, $\xi^{a}$ and $\eta$ are unknown functions of $\bf
x$.

By definition operators $Q$ have to commute $\hat{H}$:

$[\hat{H},Q]\equiv\hat{H}Q-Q\hat{H}=0.$ (8)

Condition (8) represents the operator equation which can be satisfied when the
operators presented here act on an arbitrary twice differentiable function.
Evaluating the commuta-tor and equating the coefficients for the same
differential operators $\frac{\partial}{\partial
x_{a}},\frac{\partial^{2}}{\partial x_{a}\partial x_{b}}$ and
$\frac{\partial^{3}}{\partial x_{a}\partial x_{b}\partial x_{c}}$ we come to
the following system of determining equations:

$\displaystyle
5\left(\mu^{ab}_{c}+\mu^{ac}_{b}+\mu^{bc}_{a}\right)=\delta^{ab}\left(\mu^{nn}_{c}+2\mu^{cn}_{n}\right)+\delta^{bc}\left(\mu^{nn}_{a}+2\mu^{an}_{n}\right)+\delta^{ac}\left(\mu^{nn}_{b}+2\mu^{bn}_{n}\right),$
(9) $\displaystyle\left(\mu^{nn}_{a}+2\mu^{na}_{n}\right)f-5\mu^{an}f_{n}=0,$
(10) $\displaystyle
2f\eta_{a}+\xi^{a}_{n}f_{n}-\xi^{n}f_{an}+f\xi^{a}_{nn}+2\mu^{an}V_{n}-\mu^{mn}f_{mna}=0,$
(11)
$\displaystyle\left(\mu^{ab}_{nn}+\xi^{a}_{b}+\xi^{b}_{a}\right)f+\mu^{ab}_{n}f_{n}-\mu^{na}f_{nb}-\mu^{nb}f_{na}-\delta^{ab}\left(\mu^{mn}f_{mn}+\xi^{n}f_{n}\right)=0,$
(12) $\displaystyle
f\left(\mu^{mm}_{nn}+2\xi^{n}_{n}\right)+\left(\mu^{nn}_{m}-3\xi^{m}\right)f_{m}-5\mu^{mn}f_{mn}=0,$
(13) $\displaystyle(f\eta_{n})_{n}+\xi^{n}V_{n}+\mu^{mn}V_{mn}=0$ (14)

where $f_{n}=\frac{\partial f}{\partial x_{n}},\
\xi^{a}_{n}=\frac{\partial\xi^{a}}{\partial x_{n}}$, etc.

To find all operators (2) which admit integrals of motion (7) it is necessary
to find all inequivalent solutions of the very complicated system of equations
(9)–(14) for twelve unknowns $\mu^{ab},\xi^{a},\eta,f$ and $V$. Fortunately,
equations (12), (13) and (14) can be omitted since they are nothing but the
differential consequences of the remaining ones. In addition, variables ,
$\xi^{a}$ can be excluded. Indeed, diffirentiating (10) with respect to
$x_{a}$ and making the summation over the repeating index $a$, we come to
equation (13), while the same procedure with (11) lead to equation (14). On
the other hand, equation (13) can be deduced from (9) and (10), provided the
following relations are satisfied:

$\displaystyle\begin{split}&\hat{\xi}^{b}_{a}+\hat{\xi}^{a}_{b}=\frac{2}{3}\delta_{ab}\hat{\xi}^{n}_{n},\\\
&3\hat{\xi}^{n}f_{n}=2f\hat{\xi}^{n}_{n}\end{split}$ (15)

where $\hat{\xi}^{a}=\xi^{a}-\mu^{an}_{n}$ .

Relations (LABEL:de1) are nothing but the determining equations for the
coefficients of the first order symmetry operators of equations (1), found in
paper [22]. Since such operators and the corresponding equations (1) are known
[22], we suppose that functions $\hat{\xi}^{a}$ are trivial, and so
$\xi^{a}=\mu^{an}_{n}$. As a result equations (11) are reduced to the
following form:

$\displaystyle
2(f\eta^{a}+\mu^{ab}V_{b})+(\mu^{am}_{nm}f-\mu^{nm}f_{am})_{n}=0,$ (16)

and the corresponding integral of motion (7) can be rewritten as:

$\displaystyle Q=\partial_{b}\mu^{ab}\partial_{a}+\eta.$ (17)

The corresponding functions $\mu^{ab}$ have to satisfy the autonomous system
of equations (9), which specify the conformal Killing tensor. The latter one
is a linear combination of the following tensors (see, e.g., [28]):

$\displaystyle\begin{split}&\mu^{ab}_{1}=\lambda_{1}^{ab}+\delta^{ab}\frac{\lambda_{1}^{cd}x^{c}x^{d}}{x^{2}}\varphi_{1}(x),\\\
&\mu^{ab}_{2}=\lambda_{2}^{a}x^{b}+\lambda_{2}^{b}x^{a}+\delta^{ab}\lambda_{3}^{c}x^{c}\varphi_{2}(x),\\\
&\mu^{ab}_{3}=(\varepsilon^{acd}\lambda_{3}^{cb}+\varepsilon^{bcd}\lambda_{3}^{ca})x^{d},\\\
&\mu^{ab}_{4}=(x^{a}\varepsilon^{bcd}+x^{b}\varepsilon^{acd})x^{c}\lambda_{4}^{d},\\\
&\mu^{ab}_{5}=\delta^{ab}x^{2}\varphi_{5}(x)+k(x^{a}x^{b}-\delta^{ab}x^{2}),\\\
&\mu^{ab}_{6}=\lambda_{5}^{ab}x^{2}-(x^{2}\lambda_{5}^{bc}+x^{b}\lambda_{5}^{ac})x^{c}+\delta^{ab}\lambda_{6}^{cd}x^{c}x^{d}\varphi_{6}(x),\end{split}$
(18)
$\displaystyle\begin{split}&\mu^{ab}_{7}=(x^{a}\lambda_{7}^{b}+x^{b}\lambda^{a})x^{2}-4x^{a}x^{b}\lambda_{7}^{c}x^{c}+\delta^{ab}\lambda_{8}^{c}x^{c}x^{2}\varphi_{7}(x),\\\
&\mu^{ab}_{8}=2(x^{a}\varepsilon^{bcd}+x^{b}\varepsilon^{acd})\lambda_{8}^{dn}x^{c}x^{n}-(\varepsilon^{ack}\lambda_{8}^{bk}+\varepsilon^{bck}\lambda_{8}^{ak})x^{c}x^{2},\\\
&\mu^{ab}_{9}=\lambda_{9}^{ab}x^{4}-2(x^{a}\lambda_{9}^{bc}+x^{b}\lambda_{9}^{ac})x^{c}x^{2}+(4x^{a}x^{b}+k\delta^{ab}x^{2})\lambda_{9}^{cd}x^{c}x^{d}\\\
&+\delta^{ab}\lambda_{10}^{cd}x^{c}x^{d}x^{2}\varphi_{9}(x)\end{split}$ (19)

where $\lambda_{n}^{ab}=\lambda_{n}^{ba}$ and $\lambda_{n}^{a}$ are arbitrary
parameters, and $\varphi_{1},...,\varphi_{9}$ are arbitrary functions of
$x=r=\sqrt{x_{1}^{2}+x_{2}^{2}+x_{3}^{2}}$.

Thus our classification problem is reduced to finding inequivalent solutions
of equations (10) and (16) for unknowns $f$ and $V$, where $\mu^{ab}$ is a
linear combination of functions (18) and (19). The main difficulty is the
extended number of arbitrary parameters involved into the mentioned equations
which has to be reduced using the equivalence transformations.

##  5 Decoupling of the determining equations

Taking into account the scale invariance of equations (1) with arbitrary
elements (3) we can conclude that the Killing tensors $\mu^{ab}$ involved into
the determining equations should be homogeneous functions of $\bf x$. In other
words, equations , (10) and (16) are decoupled to five autonomous subsystems
which correspond to the Killing tensors including homogeneous polynomials of
fixed order $n$ with $n=0,1,2,3,4$ , and the corresponding functions
$\varphi_{1},\varphi_{2},...,\varphi_{9}$ should be constants.

Moreover since any equation (1), (2) which arbitrary element (3) is invariant
w.r.t. the inverse transformation (5), we can restrict ourselves to the
Killing tensors being polynomials of $x_{a}$ of order $n\leqslant 2$, which
are given by equations (18), since the symmetries with $n$=3 and $n$=4 are
equivalent to the symmetries with $n=1$ and $n=0$ correspondingly. It means
that it is sufficient to solve the determining equations (10) and (16) with
the following $\mu^{ab}$:

$\displaystyle\mu^{ab}=\lambda^{ab}+\kappa\delta^{ab}\frac{\tilde{\lambda}^{cd}x^{c}x^{d}}{x^{2}},$
(20)
$\displaystyle\begin{split}&\mu^{ab}=\lambda^{a}x^{b}+\lambda^{b}x^{a}-2\delta^{ab}\tilde{\lambda}^{c}x^{c}\\\
&+\mu^{a}x^{b}+\mu^{b}x^{a}+(\varepsilon^{acd}\lambda^{cb}+\varepsilon^{bcd}\lambda^{ca})x^{d},\end{split}$
(21) $\displaystyle\begin{split}&\mu^{ab}=\kappa
x^{a}x^{b}+(x^{a}\varepsilon^{bcd}+x^{b}\varepsilon^{acd})\lambda^{d}x^{c}+\delta^{ab}\tilde{\lambda}^{cd}x^{c}x^{d}\\\
&+\lambda^{ab}x^{2}-(x^{a}\lambda^{bc}+x^{b}\lambda^{ac})x_{c}.\end{split}$
(22)

Notice than any second order symmetry corresponding to $n=0$ and $n=1$ is
accompanied by the addition symmetry generated by the changes of variables
(6).

###  5.1 Symmetries independent on $x_{a}$

Let us start with the symmetries which correspond to Killing tensors
$\mu^{ab}$ fixed in (20). Substituting (20) into (13), we obtain the following
equation:

$\displaystyle\lambda^{ab}f_{b}+\kappa\frac{\tilde{\lambda}^{mn}x_{m}x_{n}}{r^{2}}f_{a}=0.$
(23)

Since $\mu^{ab}=\mu^{ba}$, up to rotation transformations there are three
inequivalent versions of nontrivial parameters $\lambda^{ab}$:

$\displaystyle\lambda^{11}=k_{1},\lambda^{22}=k_{2},\lambda^{33}=k_{3},$ (24)
$\displaystyle\lambda^{11}=k_{1},\lambda^{22}=k_{2},$ (25)
$\displaystyle\lambda^{33}=k_{3}.$ (26)

Substituting (24)-(26) into (23) and equating the coefficients for the same
independent variables $x_{m}$, we come to the conclusion that coefficient
$\kappa$ should be trivial, and

$\displaystyle f=0\ \ \text{ for version }(\ref{deq111}),$ (27) $\displaystyle
f=x_{3}^{2}\ \ \text{ for version }(\ref{deq121}),$ (28) $\displaystyle
f=(x_{1}^{2}+x_{2}^{2})F(\varphi)\ \ \text{ for version }(\ref{deq131})$ (29)

where $F(\varphi)$ is an arbitrary function of the Euler angle
$\varphi=\arctan\left(\frac{x_{2}}{x_{1}}\right).$

Substituting (56)-(29) into (16) we obtain the following solutions for the
latter equation:

$\displaystyle f=x_{3}^{2},\quad
V=\frac{x_{3}^{2}}{(x_{1}^{2}+x_{2}^{2})}F(\varphi),\quad\eta=-\frac{1}{x_{1}^{2}+x_{2}^{2}}F(\varphi),\quad\lambda^{11}=\lambda^{22}=1,$
(30) $\displaystyle f=x_{3}^{2},\quad
V=Const,\quad\eta=0,\quad\lambda^{11}=k_{1},\ \lambda^{22}=k_{2}\neq k_{1},$
(31) $\displaystyle f=(x_{1}^{2}+x_{2}^{2})F(\varphi),\quad
V=\frac{f}{x_{3}^{2}},\quad\eta=-\frac{1}{x_{3}^{2}},\quad\lambda^{33}=1.$
(32)

Thus up to the equivalence there are three versions of the scale invariant PDM
Schrödin-ger equations which admit second order integrals of motion
independent on $\bf x$. They include arbitrary elements $f$ and $V$ fixed in
formulae (30) and (31).

###  5.2 Symmetries linear in independent variables

Consider now the symmetries generated by the Killing tensors (22) which are
linear in $x_{a}$. The related determining equations (13) are reduced to the
following form:

$\displaystyle 2\lambda^{a}f=\mu^{ab}f_{b}.$ (33)

To simplify calculations we will use the identity
$2f=x_{1}f_{1}+x_{2}f_{2}+x_{3}f_{3}$ which makes it possible to reduce (33)
to the following homogeneous system of linear algebraic equations for
deraviatives $f_{a}$:

$\displaystyle M^{ab}f_{b}=0,$ (34)

where

$\displaystyle M^{ab}=\mu^{ab}-\lambda^{a}x^{b}-\mu^{a}x_{b}.$ (35)

Notice that operator (7) with $\mu^{ab}$ given in (31) is a bilinear
combination of generators of conformal group , $C(3)$, i.e.,

$\displaystyle\begin{split}&Q=a_{+}P_{2}L_{3}+a_{-}L_{2}P_{3}+b_{+}P_{1}L_{3}+b_{-}P_{3}L_{1}+c_{+}P_{1}L_{2}\\\
&+c_{-}P_{2}L_{1}+\lambda^{a}p_{a}D+\tilde{d}_{1}P_{1}L_{1}+\tilde{d}_{2}P_{2}L_{2}+\tilde{d}_{3}L_{3}P_{3},\end{split}$
(36)

where
$a_{\pm}=\lambda^{23}\pm\lambda^{1},b_{\pm}=\lambda^{31}\pm\lambda^{2},c_{\pm}=\lambda^{12}\pm\lambda^{3},\tilde{d}_{a}=\lambda^{bb}-\lambda^{cc}$,
(a, b, c) is the cycle (1, 2, 3).

Using rotation transformations coefficients $a_{1},b_{1}$ and $c_{2}$ can be
reduced to zero. As a result the components of tensor $M^{ab}$ (35) take the
following form:

$\displaystyle\begin{split}&M^{11}=-2cx_{3}+\lambda_{1}x_{1},\
M^{12}=\lambda_{2}x_{1}+d_{3}x_{3},M^{21}=\lambda_{1}x_{2}+d_{3}x_{3},\\\
&M^{22}=\lambda_{2}x_{2},\ M^{23}=d_{1}x_{1}+\lambda_{3}x_{2}+bx_{3},\\\
&M^{13}=ax_{3}+d_{2}x_{2}+\lambda_{3}x_{1},\
M^{31}=cx_{1}+d_{2}x_{2}+(a+\lambda_{1})x_{3},\\\
&M^{32}=(b+\lambda_{2})x_{3}+d_{1}x_{1},\
M^{33}=\lambda_{3}x_{3}-2ax_{1}-2bx_{2}\end{split}$ (37)

where we redenote $a_{-}=a,\ b_{-}=b,\ c_{+}=c,\
d_{1}=\tilde{d}_{2}-\tilde{d}_{3},\ d_{2}=\tilde{d}_{3}-\tilde{d}_{1},\
d_{3}=\tilde{d}_{1}-\tilde{d}_{2}$.

Equation (34) admits nontrivial solution iff the determinant of the matrix
whose entries $M^{ab}$ are given by equations (37) is equal to zero It is the
case provided one of the following conditions is satisfied:

$\displaystyle a=b=c=d_{a}=0,$ (38) $\displaystyle d_{1}=-d_{2},c=0,\
\lambda_{1}=\lambda_{2}=0,$ (39) $\displaystyle ad_{2}=-bc,\
b\lambda_{1}=-2\lambda_{3},$ (40) $\displaystyle b\lambda_{1}=a\lambda_{2},\
c=0,$ (41) $\displaystyle a=b=d_{a}=\lambda_{2}=0,$ (42)
$\displaystyle\lambda_{1}d_{1}=-c\lambda_{2},\
d_{1}\lambda_{3}=-a\lambda_{2},\ b=0.$ (43)

Substituting consequently the versions presented in (38)– (43) into equations
(34) and (16) with matrices (37) we find the corresponding functions $f$, $V$
and $\eta$.

Let us start with version (38). In this case the symmetry operator includes a
linear combination of generators $P_{1}D,P_{2}D$ and $P_{3}D$ with
coefficients $\lambda_{1},\lambda_{2}$ and $\lambda_{3}$ correspon-dingly.
Rotating the coordinate system, two of these coefficients (say, $\lambda_{1}$
and $\lambda_{2}$) can be nullified which reduce the tensor components given
in (37) to the following form:

$\displaystyle M^{3a}=x_{a},\ a=1,2,3,$

while the remaining components are reduced to zero. The related equations (34)
and (16) take the following form:

$\displaystyle x_{a}f_{3}=0,$ (44) $\displaystyle f\eta_{a}+x_{a}V_{3}=0.$
(45)

The generic solution of equation (44) for function $f$ given in (3) is:

$\displaystyle f=\tilde{r}^{2}F(\varphi),$ (46)

where $F(\varphi)$ is an arbitrary function of the Euler angle. Integrating
the related equation (45) for $V$ given in (3), we find the generic form of
functions $\eta$ and $V$:

$\displaystyle
V=G(\varphi)+c_{1}\frac{x_{3}}{\tilde{r}}F(\varphi)+c_{2}\frac{x_{3}}{r}F(\varphi)$
$\displaystyle\eta=\frac{c_{1}}{\tilde{r}}+\frac{c_{2}}{r}$

where $c_{1}$ and $c_{2}$ are integration constants, and $G(\varphi)$ is an
additional arbitrary function.

The found solutions are represented in the firs item of Table 1, where the
additional integral of motion obtained by the inversion transformation is
fixed also.

Considering version (39) we come to the following nontrivial components of
tensor $M^{ab}$ (37):

$\displaystyle M^{13}=ax_{3}+dx_{2}+\lambda_{3}x_{1},\
M^{23}=\lambda_{3}x_{2}-dx_{1}-bx_{3},$ $\displaystyle M^{31}=ax_{3}+dx_{2},\
M^{32}=-dx_{1}-bx_{3},\ M^{33}=\lambda_{3}x_{3}+2bx_{2}-2ax_{1}.$

Up to rotation transformations we can set $b=0$, and the related equations
(34) are reduced to the following ones:

$\displaystyle f_{3}=0,\quad af_{1}=0,\quad d(x_{1}f_{2}-x_{2}f_{1})=0.$ (47)

In accordance with (47) one out of two coefficients $a$ and $d$ should be
trivial provided the mass is not a constant. The corresponding solutions of
equation (47) look as:

$\displaystyle f=x_{2}^{2},\ \text{}d=0,\ a\neq 0,$ (48) $\displaystyle
f=\tilde{r}^{2},\ \text{}a=0,\ d\neq 0.$ (49)

In the cases (48) and (49) equations (16) take the following forms:

$\displaystyle\begin{split}&x_{2}^{2}\eta_{1}+(\lambda_{3}x_{1}+ax_{3})V_{3}=0,\\\
&x_{2}^{2}\eta_{2}+\lambda_{3}x_{2}V_{3}=0,\\\
&x_{2}^{2}\eta_{3}+ax_{3}V_{1}+(\lambda_{3}x_{3}-2ax_{1})V_{3}=0\end{split}$
(50)

and

$\displaystyle\begin{split}&\tilde{r}^{2}\eta_{1}=-(dx_{2}+\lambda_{3}x_{1})V_{3},\\\
&\tilde{r}^{2}\eta_{2}=(dx_{1}-\lambda_{3}x_{2})V_{3},\\\
&\tilde{r}^{2}\eta_{3}=d(x_{1}V_{2}-x_{2}V_{1})-\lambda_{3}x_{3}V_{3}.\end{split}$
(51)

correspondingly.

If all arbitrary parameters $a,\lambda_{3}$ and $d$ are nontrivial, equations
(50) and (51) have only constant solutions for $V$ and $\eta$. Equations (50)
for $a=0$ and equations (51) for $d=0$ are reduced to the particular cases of
equations (45) and (46). Their solutions are represented in Items 2 of Table 1
where $F(\varphi)=0$ for the system (51) and $F(\varphi)=\sin(\varphi)$ for
the system (51).

Obtained in analogous manner other versions of Hamiltonians (2) which admit
second order integrals of motion linear in the independent variables are
collected in Tables 1 and 2 where the systems which admit more than one
integral of motion are indicated also.

##  6 Symmetries bilinear in independent variables

Le last task is to find PDM Schrödinger equations admitting integrals of
motion which are generated by tensors (22). These integrals of motion can be
represented as:

$\displaystyle
Q=\nu^{ab}(\\{K_{a},P_{b}\\}+\\{P_{b},K_{a}\\})+\tilde{\lambda}^{ab}Q^{ab},$
(52)

where $\nu^{ab}=\lambda^{ab}+\varepsilon_{abc}\lambda_{c},$
$\tilde{\lambda}^{ab}$ are arbitrary koefficients,
$Q^{ab}=P_{c}x_{a}x_{b}P_{c}$, and the symbols $\\{.,.\\}$ denote
anticommutators.

Operators $Q^{ab}$ also can be expressed via operators (4), since the
following identities are true:

$\displaystyle\begin{split}&\\{L_{a},L_{b}\\}+\\{P_{a},K_{b}\\}=2Q^{ab},\quad
a\neq b,\\\
&\\{P_{1},K_{1}\\}+\\{P_{2},K_{2}\\}+L_{3}^{2}=2Q^{33}.\end{split}$ (53)

The corresponding determining equations (13) take the following form:

$\displaystyle 2f(\lambda^{ab}-\tilde{\lambda}^{ab})x_{b})=\mu^{ab}f_{b}.$
(54)

In analogy with (33)-(35) it is convenient to rewrite equations (54) in the
form (34), where

$\displaystyle M^{ab}=\mu^{ab}-\lambda^{ac}x_{c}x_{b}.$ (55)

Like in the case of equation (37), using the rotation transformations many of
arbitary coefficients in (55) can be nullified and the coefficient matrix can
be reduced to one of the canonical forms.The related matrices (55) are reduced
to linear combinations of the following matrices $N^{ab}$ and
$\tilde{N}^{ab}$:

$\displaystyle M=\nu^{ab}N^{ab}+\tilde{\nu}^{ab}\tilde{N}^{ab},$ (56)

where $\tilde{N}^{ab}=x_{a}x_{b}I$, $I$ is the unit matrix, and the nontrivial
entries $N^{ab}_{cd}$ of matrices $N^{ab}$ are given in the following formulae
(refer to (22):

$\displaystyle\begin{split}&N_{11}^{12}=-x_{1}x_{2},\
N_{11}^{13}=-x_{1}x_{3},\ N_{11}^{22}=N_{11}^{33}=x_{1}^{2},\\\
&N_{22}^{11}=N_{22}^{33}=x_{2}^{2},\ N_{22}^{21}=-x_{2}x_{1},\
N_{22}^{23}=-x_{2}x_{3},\\\ &N_{33}^{31}=-x_{3}x_{1},\
N_{33}^{32}=-x_{3}x_{2},\ N_{33}^{33}=\tilde{r}^{2},\\\
&N^{12}_{11}=-2x_{1}x_{2},\ N^{12}_{12}=r^{2},\ N^{12}_{21}=r^{2}-2x_{2}^{2},\
N^{12}_{31}=-2x_{2}x_{3},\\\ &N^{21}_{12}=r^{2}-2x_{1}^{2},\
N^{21}_{21}=r^{2},\ N^{21}_{22}=-2x_{1}x_{2},\
N^{21}_{32}=-2x_{1}x_{2}.\end{split}$ (57)

Matrices (57) are degenerated and so the corresponding equations (34) and (16)
have nontrivial solutions. Inequivalent versions of these matrices are
presented in the following formulae:

$\displaystyle N=\nu_{1}N^{11}+\nu_{2}N^{22}+\nu_{3}N^{33},$ (58)
$\displaystyle N=\nu_{4}(N^{11}+N^{22})+\nu_{5}N^{33}+\nu_{6}N^{12},$ (59)
$\displaystyle N=\nu_{7}(N^{11}+N^{22}+N^{33})+\nu_{8}N^{12}+\nu_{9}N^{23},$
(60) $\displaystyle
N=\nu_{10}(N^{11}+N^{22})+\nu_{11}N^{33}+\nu_{12}(N^{12}-N^{21}),$ (61)
$\displaystyle
N=\nu_{13}(N^{11}+N^{22}+N^{33})+\nu_{14}(N^{12}-N^{21})+\nu_{15}(N^{23}-N^{32})$
(62)

where $\nu_{1},\nu_{2},...,\nu_{15}$ are arbitrary linear coefficients.

Thus the classification of the integrals of motion bilinear in $x_{a}$ is
reduced to search of generic solutions of equations (34) and (16), where $M$
are matrices given by formulae (56) and (58)-(62). We will not present the
details of the requested cumbersome calculations but restrict ourselves to one
particular example of them.

Let matrix $M$ is reduced to (61) with nontrivial $\nu_{10}$ and $\nu_{12}$.
Then equations (34) have nontrivial solutions $f$ if $\nu_{11}=0$, and
equations (16) have nontrivial solutions for $V$ if $\nu_{10}=\nu_{12}$ or
$\nu_{10}\nu_{12}=0.$

If $\nu_{10}=\nu_{11}=0$ then equations (34) and (16) take the following form:

$\displaystyle x_{1}f_{2}-x_{2}f_{1}=0$

and

$\displaystyle f\eta_{a}=-2x_{a}(x_{1}V_{2}-x_{2}V_{1}),\ a=1,2,3$

correspondingly. Their generic solutions look as:

$\displaystyle f=r^{2}F(\theta),\ V=cF(\theta)\varphi+G(\theta),\
\eta=-c\ln(r^{2}).$

If $\nu_{10}\neq 0,\nu_{12}=\nu_{11}=0$ we have the following versions of
equations (34) and (16):

$\displaystyle\begin{split}&f_{1}=0,f_{2}=0,\\\
&f\eta_{1}+(x_{2}^{2}+x_{3}^{2})V_{1}-x_{1}x_{2})V_{2}=0,\\\
&f\eta_{2}-x_{1}x_{2}V_{1}+(x_{1}^{2}+x_{3}^{2})V_{2}=0,\\\
&f\eta_{3}-x_{1}x_{3}V_{1}-x_{2}x_{3}V_{2}=0,\end{split}$

and the related solutions are:

$\displaystyle f=x_{3}^{2},\
V=\frac{x_{3}^{2}}{\tilde{r}^{2}}F(\varphi)+G(\theta),\
\eta=-\frac{r^{2}}{\tilde{r}^{2}}F(\varphi)-G(\theta).$

Table 1. Inverse masses, potentials and the related integrals of motion
defined up to arbitrary functions.

№ | $f$ | $V$ | Integrals of motion
---|---|---|---
1 | $\tilde{r}^{2}F(\varphi)$ | $F(\varphi)G(\theta)$ | $\\{P_{3},K_{3}\\}+4G(\theta)$
2 | $\tilde{r}^{2}F(\varphi)$ | $cF(\varphi)\frac{x_{3}}{r}+G(\varphi)$ | $\begin{array}[]{c}\\{P_{3},D\\}-\frac{2c}{r},\ \\{K_{3},D\\}-{2c}{r}\end{array}$
3 | $\tilde{r}^{2}F(\varphi)$ | $cF(\varphi)\frac{x_{3}}{r}$ | $\begin{array}[]{c}\\{P_{3},K_{3}\\}+\frac{4c}{r},\ \\{P_{3},D\\}-\frac{2c}{r},\\\ \\{K_{3},D\\}-{2c}{r}\end{array}$
4 | $\tilde{r}^{2}F(\varphi)$ | $F(\varphi)$ | $\begin{array}[]{c}P_{1}^{2}+P_{2}^{2}+\frac{1}{x_{3}^{2}}F(\varphi),\\\ K_{1}^{2}+K_{2}^{2}+\frac{r^{4}}{x_{3}^{2}},\\\ \\{P_{3},K_{3}\\},\ \\{P_{3},D\\},\ \\{K_{3},D\\}\end{array}$
5 | $\tilde{r}^{2}F(\varphi)$ | $\frac{\tilde{r}^{2}}{x_{3}^{2}}F(\varphi)$ | $\begin{array}[]{c}P_{3}^{2}+\frac{1}{x_{3}^{2}},\ K_{3}^{2}+\frac{r^{4}}{x_{3}^{2}},\\\ \\{P_{3},K_{3}\\}+\frac{4\tilde{r}^{2}}{x_{3}^{2}}\end{array}$
6 | $\tilde{r}^{2}F(\theta)$ | $G(\varphi)F(\theta)+R(\theta)$ | $\begin{array}[]{c}L_{3}^{2}+2G({\varphi})\end{array}$
7 | $r^{2}F(\theta)$ | $cF(\theta)\varphi+G(\theta)$ | $\begin{array}[]{c}\\{L_{3},D\\}+2c\ln(r)\end{array}$
8 | $x_{3}^{2}$ | $\frac{x_{3}^{2}}{\tilde{r}^{2}}F(\varphi)+G(\theta)$ | $\begin{array}[]{c}\\{P_{1},K_{1}\\}+\\{P_{2},K_{2}\\}\\\ +\frac{4r^{2}}{\tilde{r}^{2}}F(\varphi)+4G(\theta),\\\ L_{3}^{2}+2F({\varphi})\end{array}$
9 | $x_{3}^{2}$ | $c\frac{x_{3}^{2}\varphi}{r^{2}}+G(\theta)$ | $\begin{array}[]{c}\\{P_{1},K_{1}\\}+\\{P_{2},K_{2}\\}+4\frac{r^{2}}{\tilde{r}^{2}}\varphi+4G(\theta),\\\ L_{3}^{2}+2{\varphi},\ \\{L_{3},D\\}+c\ln(r)\end{array}$
10 | $x_{2}^{2}$ | $F(\varphi)-\frac{x_{1}^{2}}{\tilde{r}^{2}}G(\theta)$ | $\begin{array}[]{c}\\{K_{1},P_{1}\\}-2J_{2}^{2}+\frac{4x_{1}^{2}}{\tilde{r}^{2}}G(\theta)-4F(\varphi),\\\ \\{K_{3},P_{3}\\}+4G(\theta)\end{array}$
11 | $x_{3}^{2}$ | $\tilde{V}$ | $\begin{array}[]{c}\\{P_{1},K_{1}\\}+(b+1-a)\\{P_{2},K_{2}\\}\\\ +2(a-1)J_{3}^{2}+4\tilde{\eta}\end{array}$

Finally, if $\nu_{10}=\nu_{12}$, $\nu_{11}=0$, we have the following equations
(34) and (16):

$\displaystyle f_{1}=0,\ f_{2}=0,$ $\displaystyle
f\eta_{1}+(x_{2}^{2}+x_{3}^{2}-2x_{1}x_{2})V_{1}+(2x_{1}^{2}-x_{1}x_{2})V_{2}=0,$
$\displaystyle
f\eta_{2}-(x_{1}x_{2}+2x_{2}^{2})V_{1}+(x_{1}^{2}+2x_{1}x2+x_{3}^{2})V_{2}=0,$
$\displaystyle
f\eta_{3}-(x_{1}x_{3}+2x_{2}x3)V_{1}+(2x_{1}x_{3}-x_{2}x_{3})V_{2}=0,$

and the related solutions are $f=x_{3}^{2},\ V=F(\theta),\ \eta=-F(\theta)$.

Consider also the most complicated case when the symmetry operator (52) is
reduced to the following bilinear form of generators of group C(3):

$\displaystyle Q=\\{K_{1},P_{1}\\}+\mu\\{K_{2},P_{2}\\}+\nu L_{3}^{2}.$ (63)

Table 2. Inverse masses, potentials and the related integrals of motion
defined up to arbitrary coefficients.

№ | $f$ | $V$ | Integrals of motion
---|---|---|---
1 | $r^{2}$ | $c\frac{r^{2}}{x_{3}^{2}}$ | $\begin{array}[]{c}\ L_{3},\\{L_{1},L_{2}\\}+4c\frac{x_{1}x_{2}}{x_{3}^{2}}\end{array}$
2 | $x_{3}^{2}$ | $\frac{cx_{3}^{2}}{\nu^{2}r^{2}+\nu((\mu+1)x_{3}^{2}+\mu x_{1}^{2}+x_{2}^{2})+\mu x_{3}^{2}}$ | $\\{P_{1},K_{1}\\}+\mu\\{P_{2},K_{2}\\}+4\nu\frac{r^{2}}{x_{3}^{2}}V$
3 | $x_{3}^{2}$ | $\begin{array}[]{c}c_{1}\frac{x_{3}^{2}}{x_{1}^{2}}+c_{2}\frac{x_{3}^{2}}{x_{2}^{2}}+c_{3}\frac{bx_{1}^{2}+ax_{2}^{2}}{x_{3}^{2}}\\\ +c_{4}\frac{x_{3}^{2}}{2(ax_{2}^{2}-bx_{1}^{2})+(a-b)x_{3}^{2}}\end{array}$ | $\begin{array}[]{c}\\{P_{1},K_{1}\\}+(b+1-a)\\{P_{2},K_{2}\\}\\\ +2(a-1)J_{3}^{2}+4c_{1}\frac{ax_{2}^{2}+x_{3}^{2}}{x_{1}^{2}}\\\ +4c_{2}\frac{(bx_{1}^{2}+(b-a+1)x_{3}^{2}}{x_{2}^{2}}+4c_{3}\frac{ax_{1}^{2}+b(b-a+1)x_{2}^{2}}{x_{3}^{2}}\\\ +4c_{4}\frac{(2ax_{2}^{2}+x_{3}^{2})}{a(x_{3}^{2}+2x_{2}^{2})+b(x_{3}^{2}+2x_{1}^{2})}\end{array}$
4 | $x_{2}^{2}$ | $c_{1}\frac{x_{2}^{2}}{x_{3}^{2}}+c_{2}\frac{x_{1}}{\tilde{r}}$ | $\begin{array}[]{c}\\{P_{1},D\\}-\\{P_{3},L_{2}\\}-2c_{2}\frac{1}{\tilde{r}}+4c_{1}\frac{x_{1}}{x_{3}^{2}},\\\ \\{K_{1},D\\}-\\{K_{3},L_{2}\\}-2c_{2}\frac{r^{2}}{\tilde{r}}+4c_{1}\frac{r^{2}x_{1}}{x_{3}^{2}}\end{array}$
5 | $x_{2}^{2}$ | $c\frac{x_{2}^{2}}{x_{3}^{2}}$ | $\begin{array}[]{c}\\{P_{1},D\\}-\\{P_{3},L_{2}\\}+4c\frac{x_{1}}{x_{3}^{2}},\ K_{1},\\\ \\{K_{1},D\\}-\\{K_{3},L_{2}\\}+4c\frac{r^{2}x_{1}}{x_{3}^{2}},\ P_{1}\end{array}$
6 | $x_{2}^{2}$ | $c\frac{x_{1}}{\tilde{r}}$ | $\begin{array}[]{c}\\{P_{1},D\\}-\\{P_{3},L_{2}\\}-2c\frac{1}{\tilde{r}},\ K_{3},\\\ \\{K_{1},D\\}-\\{K_{3},L_{2}\\}-2c\frac{r^{2}}{\tilde{r}},\ P_{3}\end{array}$
7 | $x_{3}^{2}$ | $c_{1}\frac{x_{3}^{2}x_{2}}{\tilde{r}x_{1}^{2}}+c_{2}\frac{x_{3}^{2}}{x_{1}^{2}}$ | $\begin{array}[]{l}\\{L_{3},P_{1}\\}+4c_{2}\frac{x_{2}}{x_{1}^{2}}+2c_{1}\frac{2x_{2}^{2}+x_{1}^{2}}{\tilde{r}x_{1}^{2}},\\\ \\{L_{3},K_{1}\\}+4c_{2}\frac{x_{2}r^{2}}{x_{1}^{2}}+2c_{1}\frac{(2x_{2}^{2}+x_{1}^{2})r^{2}}{\tilde{r}x_{1}^{2}}\end{array}$
8 | $x_{3}^{2}$ | $c\frac{x_{3}^{2}}{x_{1}^{2}}$ | $\begin{array}[]{l}\\{L_{3},P_{1}\\}+4c\frac{x_{2}}{x_{1}^{2}},\ \\{L_{3},K_{1}\\}+4c\frac{x_{2}r^{2}}{x_{1}^{2}},\\\ \\{P_{1},K_{1}\\}+\frac{4cr^{2}}{x_{1}^{2}},\ P_{2},\ K_{2}\end{array}$
9 | $x_{3}^{2}$ | $c\frac{\tilde{r}^{2}}{r^{2}}$ | $\begin{array}[]{c}\\{K_{1},P_{2}\\}+4c\frac{x_{1}x_{2}}{r^{2}},\ \\{K_{2},P_{1}\\}+4c\frac{x_{1}x_{2}}{r^{2}},\\\ \\{P_{1},K_{1}\\}+4c\frac{x_{1}^{2}}{r^{2}},\ L_{3},\ \\{P_{2},K_{2}\\}+4c\frac{x_{2}^{2}}{r^{2}}\end{array}$
10 | $x_{3}^{2}$ | $c\frac{x_{3}^{2}}{r^{2}}$ | $\begin{array}[]{c}\\{P_{2},K_{1}\\}-4c\frac{x_{1}x_{2}}{r^{2}},\ \\{P_{1},K_{2}\\}-4c\frac{x_{1}x_{2}}{r^{2}},L_{3},\\\ \\{P_{1},K_{1}\\}-\frac{4cx_{1}^{2}}{r^{2}},\ \\{P_{2},K_{2}\\}-\frac{4cx_{2}^{2}}{r^{2}}\end{array}$
11 | $\tilde{r}^{2}$ | $c_{1}e^{-2\varphi}\frac{r^{2}+x_{3}^{2}}{\tilde{r}^{2}}+c_{2}e^{-\varphi}\frac{x_{3}}{\tilde{r}}$ | $\begin{array}[]{l}\\{P_{3},(L_{3}+D)\\}-4c_{1}e^{-2\varphi}\frac{x_{3}}{\tilde{r}^{2}}-2c_{2}e^{-\varphi}\frac{1}{\tilde{r}},\\\ \\{K_{3},(L_{3}+D)\\}-4c_{1}e^{-2\varphi}\frac{r^{2}x_{3}}{\tilde{r}^{2}}-2c_{2}e^{-\varphi}\frac{r^{2}}{\tilde{r}}\end{array}$
12 | $\tilde{r}^{2}$ | $c_{1}e^{2\varphi}\frac{r^{2}+x_{3}^{2}}{\tilde{r}^{2}}+c_{2}e^{\varphi}\frac{x_{3}}{\tilde{r}}$ | $\begin{array}[]{l}\\{P_{3},(L_{3}-D)\\}+4c_{1}e^{2\varphi}\frac{x_{3}}{\tilde{r}^{2}}+2c_{2}e^{\varphi}\frac{1}{\tilde{r}},\\\ \\{K_{3},(L_{3}-D)\\}+4c_{1}e^{2\varphi}\frac{r^{2}x_{3}}{\tilde{r}^{2}}+2c_{2}e^{\varphi}\frac{r^{2}}{\tilde{r}}\end{array}$
13 | $\tilde{r}^{2}$ | $c\frac{x_{3}}{r}$ | $\begin{array}[]{c}\\{P_{3},K_{3}\\}+4c\frac{x_{3}}{r},\ \\{P_{3},D\\}-\frac{2c_{2}}{r},\\\ \\{K_{3},D\\}-\frac{2c_{1}r^{2}}{\tilde{r}}-{2c_{2}}{r},\ L_{3}\end{array}$
14 | $\tilde{r}^{2}$ | $\frac{\tilde{r}^{2}}{x_{3}^{2}}$ | $\begin{array}[]{c}P_{3}^{2}+\frac{1}{x_{3}^{2}},\ K_{3}^{2}+\frac{r^{4}}{x_{3}^{2}},\ L_{3}\end{array}$
15 | $\tilde{r}^{2}$ | $c$ | $P_{3},\ L_{3},\ D,\ K_{3}$
16 | $x_{3}^{2}$ | $c$ | $P_{1},\ P_{2},\ K_{1},\ K_{2},D,L_{3}$
17 | $r^{2}$ | $c$ | $L_{1},\ L_{2},\ L_{3},D$

The corresponding matrix $M$ is degenerated, and its nonzero entries (37) take
the following form:

$\displaystyle M^{11}=2(x_{3}^{2}+(1+\nu)x_{2}^{2}),\
M^{12}=-2(\mu+\nu)x_{1}x_{2},\ M^{21}=-2(1+\nu)x_{1}x_{2},$ $\displaystyle
M^{22}=2(\mu+\nu)x_{1}^{2}+\mu x_{3}^{2},\ M^{13}=-2x_{1}x_{3},\ M^{23}=-2\mu
x_{2}x_{3}.$

The related equations (34) are solved by $f=x_{3}^{2}$, and the corresponding
equations (13): take the following form:

$\displaystyle\begin{split}&(ax_{2}^{2}+x_{3}^{2})V_{1}-bx_{1}x_{2}V_{2}=x_{3}^{2}\eta_{1},\\\
&(bx_{1}^{2}+(b+1-a)x_{3}^{2})V_{2}-bx_{1}x_{2}V_{1}=x_{3}^{2}\eta_{2}\end{split}$
(64)

where $a=1+\nu$ and $b=\mu+\nu$. Notice that the third component of equations
(13) in our case is a consequence of the system (64) since $\eta$ should
satisfy the condition $x_{a}K_{a}=0$.

By definition potential $V$ should be scale invariant and so can be treated as
a function of two scale invariant variables $y_{1}=\frac{x_{1}}{x_{3}}$ and
$y_{2}=\frac{x_{2}}{x_{3}}$. The system (64) is compatible provided the
following second order equation for $V$ is satisfied:

$\displaystyle(ay_{2}^{2}-by_{1}^{2}+a-b)V_{y_{1}y_{2}}+y_{1}y_{2}(aV_{y_{1}y_{1}}-bV_{y_{2}y_{2}})+3(ay_{2}V_{y_{1}}-by_{1}V_{y_{2}})=0.$
(65)

The generic solution of the system (64) can be found in a a closed form only
for special combinations of parameters $a$ and $b$, namely, $a=b$ and $a=0$
(or $b=0$ which is the same). They are represented in Items 8-10 of Table 1.
For $a$ and $b$ arbitrary we were able to find only particular solutions which
are presented in Items 2 and 3 of Table 2. The arbitrary solutions are denoted
as $\tilde{V}$ and $\tilde{\eta}$ and represented formally in Item 11 of Table
1.

Analogously it is possible to solve the equations corresponding to the
remaining inequi-valent matrices $M$. Doing this it is necessary to extend
matrices $N$ by adding special matrices $\tilde{N}$ (refer to equation (24))
which keeps $M$ degenerated. The obtained in this way results are collected in
two tables one of which includes systems defined up to arbitrary functions,
the other one includes arbitrary coefficients. The presen-ted list of PDM
systems admitting second order integrals of motion is complete up to rotation
transformations.

In the tables $F(.),G(.)$ and $R(.)$ are arbitrary functions of the arguments
specified in brackets, $c,c_{1},c_{2},\mu$ and $\nu$ are arbitrary real
parameters $\varphi$ and $\theta$ are Euler angles,
$r^{2}=x_{1}^{2}+x_{2}^{2}+x_{3}^{2},\ \tilde{r}^{2}=x_{1}^{2}+x_{2}^{2},$
$P_{a},K_{a},D$ and $L_{3}$ are operators collected in (4), and the summation
is imposed over the repeating indices $a$ by values 1, 2 and 3. The symbol
$\\{A,B\\}$ denotes the anticommutator of operators $A$ and $B$, i.e.,
$\\{A,B\\}=AB+BA.$ For the systems fixed in the three latest items of Table 2
the second order integrals of motion are reduced to bilinear combinations of
the first order symmetries indicated there.

##  7 Concluding remarks

We find all inequivalent PDM Schrödinger equations which are scale invariant
and admit second order integrals of motion. In the classification tables
thirteen versions of such equations are collected, part of which include
arbitrary parameters which can be treated as coupling constants (see Table 2).
The remaining equations presented in Table 1 are defined up to arbitrary
functions of the reduced numbers of independent variables.

Notice that the systems presented in Items 7, 12 of Table 1 and Items 1, 3, 7,
8, 12, 13 are invariant w.r.t. the rotations around the third coordinate axis.

Since all found hamiltonians by definition commute with the dilatation
generator $D$, the related Schrödinger equations are integrable if they admit
two second order integrals of motion, and superintegrable if the number of
such integrals of motion is more extended. In accordance with the tables there
exist a rather extended number of integrable and superintegrable PDM systems
in the considered class.

Thus we present a number of new integrable PDM Schrödinger equations. The main
value of this result is its completeness, since we find all inequivalent
equations in the considered class. It can be considered as the first step in
the description of all PDM quantum systems which admit second order integrals
of motion and are invariant w.r.t. at least one parametrical Lie group.

A natural next step is construction of exact solutions of the obtained
integrable and superintagrable systems and the description of analogous
systems invariant w.r.t. all inequivalent one parametric lie groups. The
latter task is not too cumbersome since the number of such groups is not too
extended and is equal to four [22]. Let us remind that the exact solutions for
the PDM Schrödinger equations admitting at least five parametric symmetry
groups were found in [29].

## References

  * [1] V. I. Fushchich, and A. G. Nikitin. Symmetries of equations of quantum mechanics. New York: Allerton Press (1994).
  * [2] W. I Fushchych and A. G. Nikitin, On the new invariance groups of the Dirac and Kemmer-Duffin-Petiau equations, Lettere al Nuovo Cimento 19, 347-352 (1977).
  * [3] A. G. Nikitin, The maximal ”kinematical” invariance group for an arbitrary potential revised, Journal of Mathematical Physics, Analysis, Geometry 14, 519-531 (2018).
  * [4] A. G. Nikitin, Symmetries of Schrödinger equation with scalar and vector potentials, J. Phys. A: 53, 455202 (2020).
  * [5] A. G. Nikitin, Symmetries of the Schrödinger-Pauli equation for neutral particles, J. Math. Phys. 62, 083509 (2021).
  * [6] A. G, Nikitin Symmetries of the Schrödinger-Pauli equations for charged particles and quasirelativistic Schrödinger equations, J. Phys. A 55, 115202 (2022).
  * [7] U. Niederer, The maximal kinematical invariance group of the free Schrödinger equations, Helv. Phys. Acta, 45, 802–810 (1972).
  * [8] R. L. Anderson, S. Kumei, C. E. Wulfman, Invariants of the equations of wave mechanics. I., Rev. Mex. Fis. 21, 1–33 (1972).
  * [9] C. P. Boyer, The maximal kinematical invariance group for an arbitrary potential, Helv. Phys. Acta 47, 450–605 (1974).
  * [10] P. Winternitz, J. Smorodinsky, M. Uhlir̆ and I. Fris̆, Symmetry groups in classical and quantum mechanics, Sov. J. Nucl. Phys. 4, 444–450 (1967).
  * [11] A. Makarov, J. Smorodinsky, Kh. Valiev and P. Winternitz A systematic search for non-relativistic systems with dynamical symmetries, Nuovo Cim. A 52, 1061–1084 (1967)
  * [12] Ian Marquette and Pavel Winternitz. Higher order quantum superintegrability: a new Painleve conjecture. Integrability, Supersymmetry and Coherent States. Springer, Cham, pp. 103-131. (2019).
  * [13] A. G. Nikitin, Higher-order symmetry operators for Schrödinger equation. In CRM Proceedings and Lecture Notes (AMS), 37 , pp. 137–144 (2004).
  * [14] Oldwig von Roos, Position-dependent effective masses in semiconductor theory, Phys. Rev. B 27, 7547 (1983).
  * [15] A. de Saavedra, F. Boronat, A. Polls and A. Fabrocini, Effective mass of one He 4 atom in liquid He 3, Phys. Rev. B 50, 4248 (1994).
  * [16] P. Harrison, Quantum Wells, Wires and Dots (Wiley, New York, 2000).
  * [17] R. Heinonen, E. G. Kalnins, W. Miller Jr and E. Subag, Structure relations and Darboux contractions for 2D 2nd order superintegrable systems, SIGMA 11, 043 (2015).
  * [18] B. K. Berntson, E. G. Kalnins, and W. Miller Jr. Toward classification of 2nd order superintegrable systems in 3-dimensional conformally flat spaces with functionally linearly dependent symmetry operators, SIGMA: Symmetry, Integrability and Geometry: Methods and Applications, 16, 135. (2020).
  * [19] A. Ballesteros, A. Enciso, F. J. Herranz, O. Ragnisco and D. Riglioni, Superintegrable oscillator and Kepler systems on spaces of nonconstant curvature via the Stäckel transform, SIGMA 7 048 (2011).
  * [20] O. Ragnisco and D. Riglioni, A Family of Exactly Solvable Radial Quantum Systems on Space of Non-Constant Curvature with Accidental Degeneracy in the Spectrum, SIGMA 6 097 ( 2010).
  * [21] A. G. Nikitin, Superintegrable and shape invariant systems with position dependent mass, J. Phys. A: Math. Theor. 48 335201 (2015).
  * [22] A. G. Nikitin and T. M. Zasadko, Superintegrable systems with position dependent mass, J. Math. Phys. 56 042101 (2015).
  * [23] A. G. Nikitin, Kinematical invariance groups of the 3d Schrödinger equations with position dependent masses, J. Math. Phys. 58, 083508 (2017).
  * [24] A. G. Nikitin, Group classification of systems of non-linear reaction-diffusion equations with general diffusion matrix. I. Generalized Ginzburg-Landau equations, J. Math. Analysis and Appl. 324, 615-628 (2007).
  * [25] A. G. Nikitin, Group classification of systems of nonlinear reaction-diffusion equations with triangular diffusion matrix, Ukr. Math. J. 59, 439-458 (2007).
  * [26] A. G. Nikitin and V. I. Fushchich, Equations of motion for particles of arbitrary spin invariant under the Galileo group, Theor. and Math. Phys. 44, 584-592 (1980)
  * [27] O. O. Vaneeva, R. O. Popovych and C. Sophocleous, Equivalence transformations in the study of integrability, Physica Scripta, 89, 038003 (2014).
  * [28] A. G. Nikitin, Generalized Killing tensors of arbitrary valence and order, Ukr. Math. J. 43 734–743 (1991).
  * [29] A.G. Nikitin, Exact solvability of PDM systems with extended Lie symmetries, Proceedings of the Institute Mathematics of the Nat. Acad. Sci. of Ukraine 16, No 1, 1–18 (2019).


Citation

Y. Logan, K. Kokilepersaud, G. Kwon and G. AlRegib, C. Wykoff, H. Yu, ”Multi-
Modal Learning Using Physicians Diagnostics for Optical Coherence Tomography
Classification,” IEEE International Symposium on Biomedical Imaging (ISBI),
accepted on Jan. 7 2022.

Review

Date of accept: 7 Jan 2022

Bib

@ARTICLE{Logan2022_ISBI,
author={Y. Logan, K. Kokilepersaud, G. Kwon and G. AlRegib, C. Wykoff and H.
Yu},
journal={IEEE International Symposium on Biomedical Imaging},
title={Multi-Modal Learning Using Physicians Diagnostics for Optical Coherence
Tomography Classification},
year={2022}

Copyright

©2022 IEEE. Personal use of this material is permitted. Permission from IEEE
must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works.

Contact

<EMAIL_ADDRESS>OR<EMAIL_ADDRESS>
http://ghassanalregib.com/

# Multi-Modal Learning Using Physicians Diagnostics for Optical Coherence
Tomography Classification

###### Abstract

In this paper, we propose a framework that incorporates experts diagnostics
and insights into the analysis of Optical Coherence Tomography (OCT) using
multi-modal learning. To demonstrate the effectiveness of this approach, we
create a medical diagnostic attribute dataset to improve disease
classification using OCT. Although there have been successful attempts to
deploy machine learning for disease classification in OCT, such methodologies
lack the experts insights. We argue that injecting ophthalmological
assessments as another supervision in a learning framework is of great
importance for the machine learning process to perform accurate and
interpretable classification. We demonstrate the proposed framework through
comprehensive experiments that compare the effectiveness of combining
diagnostic attribute features with latent visual representations and show that
they surpass the state-of-the-art approach. Finally, we analyze the proposed
dual-stream architecture and provide an insight that determine the components
that contribute most to classification performance.

Index Terms—  Multi-modal Learning, Diagnostic Attributes, Latent
Representation, Autoencoder, OCT

##  1 Introduction

Optical coherence tomography (OCT) is a non-invasive imaging technique that
provides high-resolution cross-sectional images of the retina. Eye care
professionals use it to evaluate the health of the posterior segment of the
eye and to detect signs of various vitreoretinal pathologies. The physician’s
assessment when interpreting an ophthalmic image, provides valuable
information that can be a source of knowledge to a model while training.
However, current methods are not using this information. Existing approaches
for OCT classification typically entail transfer learning with pretrained
models [1, 2], as shown in Fig. 1. Though transfer learning often yields
better performance than medical experts, very rarely does this transition from
research to real clinical settings. This is because transfer learning models
make decisions based on unrelated, non-medical, internal weights that are not
directly optimised for OCT. This lack of medical intuition behind the
decision-making process inspires a lack of trust from the medical community.

Fig. 1: High-level overview of proposed multi-modal learning framework using
Physicians diagnostics as an attribute.

In this paper, we aim to mitigate these concerns by considering a physician’s
assessment as another modality to incorporate with OCT for richer feature
representations and improve retinal diagnosis performance. A modality, in this
context, is an information channel such as OCT. A multi-modal machine learning
problem uses information sourced from multiple communication channels to train
models [3]. Existing multi-modal research with OCT combine it with other
medical imaging technologies [4]. However, combining OCT with diagnostic
information is largely unexplored. Diagnostic attributes are text-based and
they characterize diseased OCT on patterns a physician would recognize.
Combining OCT with diagnostic attributes adds medical context to a model from
a multi-modal machine learning perspective. To the best of our knowledge, this
is the first attempt to combine ophthalmic imagery with diagnostic attributes
for OCT classification.

We highlight the difference between existing and proposed methods for OCT
classification in Fig. 1. Our approach involves fusion of heterogeneous
modalities trained from the outset with random weights. The challenge is to
transform the data into a meaningful format for the model to join these
modalities. Combining the uni-modal information into the same representation
space with a dual-stream autoencoder accomplishes this. The dissimilar natures
of data make it challenging to optimally exploit any complementarity and
redundancies present. Additionally, there is no closed-form relationship
between each data type so the ideal mapping of the modalities may not exist
[3]. However, multi-modal learning with diagnostic attributes adds robustness
to the feature space as they are not susceptible to variations in image
quality like OCT. Also, in practice, networks trained from random
initialization improve generalization as they avert rote memorization and
emphasize learning over time.

The diagnostic attributes dataset was acquired by mining a plethora of medical
documentation describing attributes indicative of retinal pathologies. During
data collection, each attribute noted was chronicled in multiple medical
publications. These attributes allow integration of medical insight in the
framework to improve OCT classification.

Fig. 2: A Healthy cross-sectional scan compared to scans with Diabetic Macular
Edema (DME), Choroidal Neovascularization (CNV) and Drusen retinal pathologies
with arrows highlighting some attributes of each disease.

In this paper, we create a diagnostic attribute dataset consisting of features
highly representative of three retinal pathologies. We also demonstrate how
these attributes are created. We then propose to characterize these diseases
among healthy ophthalmic images through use of this dataset along with OCT. We
create a dual-stream architecture and compare it to other architectures
trained without diagnostic attributes to show the generalizability of
diagnostic attributes for pathology classification. More specifically, we fuse
diagnostic attributes with spectral-domain OCT images in order to improve
disease classification performance. The contributions of this paper are three
fold:

  1. i

We create an ophthalmic imagery diagnostic attribute dataset with features
indicative of three retinal pathologies.

  2. ii

We propose a dual-stream autoencoder to learn joint representations for
diagnostic attributes and OCT.

  3. iii

We show that combining diagnostic attribute features with latent visual
representations surpasses the state-of-the-art approach.

##  2 Related Works

Fig. 3: Architectures for OCT classification: (A) Single-stream autoencoder
(B) Dual-stream autoencoder

Using multiple modalities for different vision tasks within the medical domain
is not new. In the ophthalmological space, this traditionally entails
combining modalities from different imaging techniques. For example, the
authors of [5] utilized information jointly between fundus images and OCT
scans in order to perform segmentation of the optic disc. More recently, work
in [6] takes a retinal image along with a structure map and show that two
modalities are better than one. Outside OCT, multiple modalities have also
been used to improve classification tasks where representations of MR and PET
scans were fused for Alzheimer’s diagnosis [7]. Neither of these utilize a
doctor’s assessment as an additional modality to integrate more supervision in
the learning framework.

Diagnostic text-based attributes allow an ophthalmologist’s assessment to be
injected into the model’s learning. However, this research task is largely
unexplored. A seldom example in [8] uses clinical textual notes as a method to
improve diagnosis performance. However, notes contain irrelevant and possibly
inaccurate information that require pre-processing. Text-based attributes from
medical reports are most commonly used with images for multi-modal document
retrieval tasks [9, 10]. The work in [11] presents an argument for merging
medical modalities due to the incomplete picture that a single modality can
represent about a patient’s condition. However, the methods they present fail
to capture direct diagnostic attributes which may offer a complete
representation.

Stacked autoencoders are unsupervised neural models that have also been used
to create feature representations for each modality and share them to learn
joint a representation [12]. We use a dual-stream autoencoder to learn joint
representations for diagnostic attributes and OCT. In the remaining sections
of this paper, we explore how diagnostic attributes impact the the learning
framework for retinal disease classification and how our autoencoder
architecture improves classification.

##  3 Multi-modal Representations

Autoencoders are the core elements of our proposed multi-modal OCT
classification. More specifically, we employ a dual-stream autoencoder,
inspired by [12], that utilizes an encoder and decoder pair for both the
visual and attribute input vectors as shown in Fig. 3. When OCT, $x$, is given
to the visual encoder, $f_{v}$, it generates compressed latent features
$z_{img}=f_{v}(x)$ that the visual decoder, $g_{v}$, uses to reconstruct its
output $\hat{x}=g_{v}(z_{img})$. When diagnostic attributes, $a$, are fed to
the attribute encoder, $f_{a}$, it too creates a compact latent variable
$z_{att}=f_{a}(a)$ that the attribute decoder, $g_{a}$, uses to reconstruct
its output $\hat{a}=g_{a}(z_{att})$.

The visual $\\{{x_{1}},...,{x_{M}}\\}$ and attribute
$\\{{a_{1}},...,{a_{M}}\\}$ inputs have both been projected onto the same
latent space and both latent representations are fed back into their own
decoder in order to reconstruct their original representation. Minimizing the
reconstruction error $\mathcal{L}_{recon}$ allows the reconstructed image and
attributes to be less differing from their input. It is defined as:

$\displaystyle\mathcal{L}_{recon}=||x-\hat{x}||_{2}+||a-\hat{a}||_{2}$ (1)

Second, these representations are then fed into the decoder of the other
latent vector. In other words, the visual latent representation undergoes
reconstruction in the attribute decoder and vice versa. The cross-
reconstruction error, $\mathcal{L}_{cross}$, enforces the model to learn the
association between OCT and diagnostic attributes by cross-reconstructing each
modality accurately. It is calculated as follows:

$\displaystyle\mathcal{L}_{cross}=||x-\hat{x}_{a}||_{2}+||a-\hat{a}_{x}||_{2}$
(2)

We also impose the latent loss, $\mathcal{L}_{dist}$, which is the cross
entropy loss of $l_{2}$ distance between latent representations of OCT and its
attribute. This loss enforces alignment between the OCT latent representation
and its attribute. $\mathcal{L}_{dist}$ is calculated as follows with
$z_{att}^{k}$ being the latent variable for the $k^{th}$ class attribute:

$\displaystyle\mathcal{L}_{dist}=-\log\left(\frac{\exp(-||z_{img}-z_{att}||_{2})}{\sum_{k=1}^{4}\exp(-||z_{img}-z_{att}^{k}||_{2})}\right)$
(3)

Training of the autoencoders is optimised by minimizing the following loss
function that has scale factors $\lambda_{cross}=1$ and
$\lambda_{dist}=0.005$:

$\displaystyle\mathcal{L}=\mathcal{L}_{recon}+\lambda_{cross}\mathcal{L}_{cross}-\lambda_{dist}\mathcal{L}_{dist}$
(4)

Disease characterization is performed with a classifier that uses latent
visual features $z_{img}$ to classify a total of four classes: healthy and
three OCT patterns indicative of specific disease states. Its loss is weighted
categorical cross entropy, $\mathcal{L}_{class}$. It is calculated as follows
with $\hat{y}[class]$ being the classifier’s predicted output for the target
$class$, $\hat{y}[k]$ being its predicted output for the $k^{th}$ class and
$w[class]$ is the inverse frequency weight for the target $class$:

$\displaystyle\mathcal{L}_{class}=w[class]\left[-\log\left(\frac{\exp([\hat{y}[class])}{\sum_{k=1}^{4}\exp(\hat{y}[k])}\right)\right]$
(5) $\displaystyle w[class]=\frac{Number\ of\ images\ in\ largest\
class}{Number\ of\ images\ in\ class}$ (6)

##  4 Attribute Dataset Collection

The diagnostic attributes dataset is intended to improve the task of OCT
classification. Therefore, all attributes considered were documented in
multiple medical publications [13, 14, 15, 16] as highly expressive of the
three pathologies. For example, by OCT, drusen are visualized as uniform,
typically dome-shaped hyperreflectvie collections beneath the retinal pigment
epithelium (RPE). We identified 26 visual indicators used by physicians when
making diagnoses based on review of OCT images. A subset of these are shown in
Table 1.

The dataset consists of diagnostic attributes for four classes: healthy,
drusen, choroidal neovascularization (CNV) and DME. Each class is associated
with an attribute vector with a length equivalent to the number of attributes.
For each position within the vector, a 1 represents the presence of that
attribute and a 0 represented its absence. The attribute for healthy OCT had 0
for all disease attributes and 1 indicating healthy. These vectors were then
associated with each image based on their class and used within the
architecture described in Section 3.

Pathology |  Manifestation in OCT
---|---
DME |  • Reduced intraretinal reflectivity • Hard exudates
CNV |  • RPE detachment • Vitreous hemmorhage
Drusen |  • Deposits underneath retina • Bumpy surface
Table 1: A subset of the diagnostic attributes from each retinal disease. Metric | Classifier | Single-Stream | Dual-Stream
---|---|---|---
Average Precision | 0.513 | 0.664 | 0.703
Average Recall | 0.499 | 0.672 | 0.690
Average F1 Score | 0.445 | 0.643 | 0.669
Table 2: Metrics comparing model performances.  Metric |  $\Delta$ Performance
---|---
Avg Precision | +2.10%
Avg Recall | +0.85%
Avg F1 Score | +1.75%
Table 3: Impact of diagnostic attributes combined with existing pretrained Inception V3 state-of-the-art approach in [17].  Ablation | Recall
---|---
No regularization | 0.665
Dropout | 0.675
Batch normalization | 0.641
Batch normalization and dropout | 0.690
Table 4: Ablation study on regularization terms for the dual-stream model

Fig. 4: Confusion matrices showing performance of the baseline (left) and
dual-stream (right) classification models.

##  5 Experiments

To evaluate the effectiveness of diagnostic attributes for disease
classification, we compare the proposed dual-stream autoencoder with two
different models. First, we train a classifier directly on OCT as a baseline.
The classifier contains two linear layers with relu activation, batch
normalization and dropout between. For the second comparison model, we train
the single-stream autoencoder as shown in Fig. 3A and later train the same
classifier with its $z_{img}$ features. Comparison between the single-stream
and the dual-stream autoencoders will highlight the contribution of diagnostic
attributes when the latent representations are utilized for disease
classification. All encoders and decoders used are symmetric with two linear
layers having a relu activation between. To investigate whether diagnostic
attribute features combined with visual features can surpass existing methods,
first, we re-implement the Inception V3 architecture, $\phi^{SOTA}$,
pretrained on ImageNet that achieved SOTA OCT classification [17] in Pytorch.
Then, we train a multilayer perceptron ,$\phi^{MLP}$, on attributes and
minimize the following objective function inspired by [18] to combine
diagnostic attributes:

$\displaystyle\mathcal{L}_{combo}=\mathcal{L}_{class}^{SOTA}+\mathcal{L}_{class}^{MLP}+\mathcal{L}_{PC}$
(7)
$\displaystyle\mathcal{L}_{PC}=1[\hat{y}^{MLP}=y]\frac{1}{2}||\phi^{SOTA}(x_{i})-\phi^{MLP}(x_{i})||_{2}^{2}$
(8)

Finally, an ablation study was conducted to show which regularization terms
contribute to learning better joint representations with diagnostic
attributes. Fig 2 shows examples for each class in the dataset [17]. Arrows
indicate sample locations for pathology in each disease. The training set
consisted of 50280 healthy, 10488 DME, 36345 CNV and 7756 Drusen samples.
Validation and test sets consisted of 861 and 250 samples from each class,
respectively. During classification, the loss function is weighted by
$w[class]$ to counteract the imbalanced dataset. All OCT were resized to
$512\times 128$, normalized and flattened before being given as input.

##  6 Results and discussion

We report average precision, recall and f1 scores. The positive impact of
combining diagnostic attributes with the SOTA method [17] is illustrated in
Table 3. The results of all other models are in Table 2. The classifier
trained on visual features from the dual-stream autoencoder outperforms both
baseline and single-stream models in all metrics. This means there is value in
incorporating meaningful knowledge from a doctor’s assessment into multi-modal
learning framework. Secondly, the baseline experiment having the poorest
performance highlights that latent representations are better for the task of
OCT disease classification.

Relative to the baseline, the dual-stream model was the best at distinguishing
healthy, Drusen, DME and CNV in OCT as shown in Fig. 4. Drusen is a
challenging class to discriminate as it contains the fewest samples in the
dataset. Furthermore, the ”bumpy” structural nature of Drusen in OCT is
sometimes subtle and can be misclassified. Nevertheless, the dual-stream model
was still the best at recognizing it.

To investigate which regularization terms of the dual-stream model contributed
most to its success, an ablation study was conducted by keeping one
regularization component and removing all others. Regularization involved
standardizing the mini-batch inputs with batch normalization and randomly
discarding some network units via dropout [19, 20]. Table 4 shows both
regularizations together had the highest impact on model recall. Following
this, dropout was the most influential regularization.

##  7 Conclusion

In this paper, we propose a framework that integrates supervision from medical
insights through multi-modal learning. We created a diagnostic attribute
dataset to intermix the diagnostic knowledge from ophthalmologists into the
machine learning framework. Our experiments validate the effectiveness of
diagnostic attributes for OCT classification. In addition, we perform an
ablation study to appreciate regularization terms that support learning
descriptive joint representations for OCT and diagnostic attributes. Finally,
we demonstrate that attributes combined with existing methods surpass SOTA.
Encouraged by these results, we will further investigate learning methods to
fully utilize diagnostic attributes for the characterization of pathologies.

##  8 Acknowledgments

This material is based upon work supported by the National Science Foundation
Graduate Research Fellowship under Grant No. DGE-1650044.

##  9 Compliance with Ethical Standards

This is a numerical simulation study for which no ethical approval was
required.

## References

  * [1] Hokuto Hirano, Akinori Minagi, and Kazuhiro Takemoto,  “Universal adversarial attacks on deep neural networks for medical image classification,”  BMC medical imaging, vol. 21, no. 1, pp. 1–13, 2021.
  * [2] Tae Keun Yoo, Joon Yul Choi, and Hong Kyu Kim,  “Feasibility study to improve deep learning in oct diagnosis of rare retinal diseases with few-shot classification,”  Medical & Biological Engineering & Computing, vol. 59, no. 2, pp. 401–415, 2021.
  * [3] Tadas Baltrušaitis, Chaitanya Ahuja, and Louis-Philippe Morency,  “Multimodal machine learning: A survey and taxonomy,”  IEEE transactions on pattern analysis and machine intelligence, vol. 41, no. 2, pp. 423–443, 2018.
  * [4] Rainer A Leitgeb and Bernhard Baumann,  “Multimodal optical medical imaging concepts based on optical coherence tomography,”  Frontiers in Physics, vol. 6, pp. 114, 2018.
  * [5] Mohammad Saleh Miri, Michael D Abràmoff, Kyungmoo Lee, Meindert Niemeijer, Jui-Kai Wang, Young H Kwon, and Mona K Garvin,  “Multimodal segmentation of optic disc and cup from sd-oct and color fundus photographs using a machine-learning graph-based approach,”  IEEE transactions on medical imaging, vol. 34, no. 9, pp. 1854–1866, 2015.
  * [6] Weixin Luo, Zaiwang Gu, Jiang Liu, and Shenghua Gao,  “Encoding structure-texture relation with p-net for anomaly detection in retinal images,”  2020\.
  * [7] Siqi Liu, Sidong Liu, Weidong Cai, Hangyu Che, Sonia Pujol, Ron Kikinis, Dagan Feng, Michael J Fulham, et al.,  “Multimodal neuroimaging feature learning for multiclass diagnosis of alzheimer’s disease,”  IEEE Transactions on Biomedical Engineering, vol. 62, no. 4, pp. 1132–1140, 2014.
  * [8] Zhi Qiao, Xian Wu, Shen Ge, and Wei Fan,  “Mnn: multimodal attentional neural networks for diagnosis prediction,”  Extraction, vol. 1, pp. A1, 2019.
  * [9] Yikuan Li, Hanyin Wang, and Yuan Luo,  “A comparison of pre-trained vision-and-language models for multimodal representation learning across medical images and reports,”  arXiv preprint arXiv:2009.01523, 2020.
  * [10] Tzu-Ming Harry Hsu, Wei-Hung Weng, Willie Boag, Matthew McDermott, and Peter Szolovits,  “Unsupervised multimodal representation learning across medical images and reports,”  arXiv preprint arXiv:1811.08615, 2018.
  * [11] Shih-Cheng Huang, Anuj Pareek, Saeed Seyyedi, Imon Banerjee, and Matthew P Lungren,  “Fusion of medical imaging and electronic health records using deep learning: a systematic review and implementation guidelines,”  NPJ digital medicine, vol. 3, no. 1, pp. 1–9, 2020.
  * [12] Jiquan Ngiam, Aditya Khosla, Mingyu Kim, Juhan Nam, Honglak Lee, and Andrew Y Ng,  “Multimodal deep learning,”  in ICML, 2011.
  * [13] Betül Tuğcu and Hakan Özdemir,  “Imaging methods in the diagnosis of optic disc drusen,”  Turkish Journal of Ophthalmology, vol. 46, no. 5, pp. 232, 2016\.
  * [14] Bartosz L Sikorski, Grazyna Malukiewicz, Joanna Stafiej, Hanna Lesiewska-Junk, and Dorota Raczynska,  “The diagnostic function of oct in diabetic maculopathy,”  Mediators of inflammation, vol. 2013, 2013.
  * [15] Sandeep Grover, Ravi K Murthy, Vikram S Brar, and Kakarla V Chalam,  “Comparison of retinal thickness in normal eyes using stratus and spectralis optical coherence tomography,”  Investigative ophthalmology & visual science, vol. 51, no. 5, pp. 2644–2647, 2010.
  * [16] JA Montero and JM Ruiz-Moreno,  “Macular thickness in patients with choroidal neovascularization determined by rta and oct3: comparative results,”  Eye, vol. 19, no. 1, pp. 72–76, 2005.
  * [17] Daniel Kermany, Kang Zhang, and Michael Goldbaum,  “Large dataset of labeled optical coherence tomography (oct) and chest x-ray images,”  Mendeley Data, v3 http://dx. doi. org/10.17632/rscbjbr9sj, vol. 3, 2018.
  * [18] Sijie Yan, Yuanjun Xiong, Kaustav Kundu, Shuo Yang, Siqi Deng, Meng Wang, Wei Xia, and Stefano Soatto,  “Positive-congruent training: Towards regression-free model updates,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2021, pp. 14299–14308.
  * [19] Sergey Ioffe and Christian Szegedy,  “Batch normalization: Accelerating deep network training by reducing internal covariate shift,”  in International conference on machine learning. PMLR, 2015, pp. 448–456.
  * [20] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov,  “Dropout: a simple way to prevent neural networks from overfitting,”  The journal of machine learning research, vol. 15, no. 1, pp. 1929–1958, 2014.


# A New Higgs Boson with Electron-Muon Flavor-Violating Couplings

R. Primulando<EMAIL_ADDRESS>Center for Theoretical Physics,
Department of Physics, Parahyangan Catholic University, Jalan Ciumbuleuit 94,
Bandung 40141, Indonesia  J. Julio<EMAIL_ADDRESS>National Research and
Innovation Agency, KST B. J. Habibie, South Tangerang 15314, Indonesia  N.
Srimanobhas<EMAIL_ADDRESS>High Energy Physics Research Unit,
Department of Physics, Faculty of Science, Chulalongkorn University,
Pathumwan, Bangkok 10330, Thailand  P. Uttayarat<EMAIL_ADDRESS>Department of Physics, Srinakharinwirot University, 114 Sukhumvit 23 Rd.,
Wattana, Bangkok 10110, Thailand

###### Abstract

A recent CMS search for a new resonance decaying to $e\mu$ in the mass range
110 GeV to 160 GeV finds an excess of events at 146 GeV. We interpret the
search results in the context of the type-III two-Higgs-doublet-model. We find
that the excess is moderately constrained by low-energy lepton-flavor-
violation constraints, in particular the $\mu\to e\gamma$ decay. We also find
the bounds from CMS search can be superior to the low-energy constraints for
the scalar mass between 110 GeV and 150 GeV, suggesting the importance of this
mass range for future searches.

##  I Introduction

As the most recently discovered particle of the Standard Model (SM) [1, 2],
the 125-GeV Higgs boson, $h$, is the least studied fundamental particle. After
the discovery, immense amount of works on precision measurements of the $h$
properties has been done at the Large Hadron Collider (LHC). It is found that
the $h$ properties agree quite well with the SM expectations [3, 4]. In
particular, since its discovery in the $\gamma\gamma$, $W^{+}W^{-}$ and $ZZ$
channels, the fermionic decay channels $b\bar{b}$ [5, 6], $\tau^{+}\tau^{-}$
[7], and $\mu^{+}\mu^{-}$ [8] have been established. The consistency of these
measurements with the SM predictions is one of the biggest triumphs of the SM.
Searches for the remaining predicted decay channels, e.g., $e^{+}e^{-}$ [9,
10] and $c\bar{c}$ [11, 12], are underway. Despite all these, one should not
discount other possible decay modes that are not predicted by the SM.
Discovering any of these decays will be a clear signal of new physics.

Lepton-flavor-violating (LFV) couplings are prototypical examples of new
physics. The LFV couplings of the $h$ lead to LFV decays
$h\to\ell\ell^{\prime}$, which are correlated with low-energy LFV decays of
charged lepton $\ell\to\ell^{\prime}\gamma$ and $\ell\to 3\ell^{\prime}$.
Additionally, the $\mu\to e$ conversion in atomic nuclei constrains any
processes involving the LFV $e$-$\mu$ coupling. In the case of LFV tau sector,
it has been demonstrated that the LHC searches for $h\to\tau\ell$ provide more
stringent constraints than the low-energy processes [13, 14]. This indicates
that collider searches for $h\to\tau\ell$, together with their heavy scalar
counterparts, are important in constraining new physics parameter space [15,
16, 17, 18, 19, 20, 21, 22].

The collider search for $h\to e\mu$, on the other hand, has often been
overlooked as tools for constraining new physics parameter space. This is due
to a conventional wisdom that such a channel is correlated with $\mu\to
e\gamma$ and $\mu\to e$ conversion, which provide more stringent constraints
[13, 14]. In this work, we show that such expectation does not generally apply
to a new resonance decaying into $e^{\pm}\mu^{\mp}$. In particular, when the
new particle mixes with the $h$ and its mass is less than 160 GeV, its
contribution naturally cancels that of the $h$ in the $\mu\to e\gamma$ and
$\mu\to e$ conversion, weakening the bounds from these processes.

The strongest bound on $h\to e\mu$ derived from the full LHC Run II data is
set by the CMS with $\text{BR}(h\to e^{\pm}\mu^{\mp})\leq 4.4\times 10^{-5}$
[23]. Beside searching for $h\to e^{\pm}\mu^{\mp}$, CMS also looks for the LFV
decay of a new resonance, $H$, for 110 GeV $<m_{H}<$ 160 GeV. The CMS finds an
excess at $m_{H}\sim 146$ GeV with a $3.8\sigma$ $(2.8\sigma)$ local (global)
significance; the best fit cross-section is found to be $\sigma\left(pp\to
H\to e^{\pm}\mu^{\mp}\right)=3.82_{-1.09}^{+1.16}$ fb. However, the
corresponding ATLAS search [10] does not find any significant excess at 146
GeV. Hence, the nature of the CMS excess remains inconclusive.

The CMS excess will inevitably induce $\mu\to e$ conversion in nuclei by a
tree-level exchange of $H$. One can estimate the $\mu\to e$ conversion rate in
a simplified model where $H$ is assumed to couple only to $e$-$\mu$ and to
gluon via an effective coupling. We find that in this scenario the induced
rate for $\mu\to e$ conversion in the gold nucleus can be as low as
$10^{-16}$, which is safely below current experimental bound. However, in a
more complete model, there will be more states that contribute to the process.
Hence, it is compelling to investigate this in a more realistic model.

In this work, we compare the CMS search results against the corresponding low-
energy bounds in the context of the type-III two-Higgs-doublet model (2HDM).
While not advocating the presence of the excess, we try to examine whether it
holds up against the current and the future $\mu\to e\gamma$ and $\mu\to e$
conversion bounds. We also highlight the importance of collider searches for
LFV decays at a relatively low mass, i.e., between 100 GeV and 150 GeV.

##  II Type-III 2HDM

The type-III 2HDM is conveniently described in the Higgs basis [24], where the
two Higgs doublets, $H_{1,2}$, are given by

$H_{1}=\begin{pmatrix}G^{+}\\\\[1.99997pt]
\dfrac{v+h_{1}+iG}{\sqrt{2}}\end{pmatrix},\quad
H_{2}=\begin{pmatrix}H^{+}\\\\[1.99997pt]
\dfrac{h_{2}+iA}{\sqrt{2}}\end{pmatrix}.$ (1)

Here $v$ is the vacuum expectation value, and $G^{+}$ and $G$ are the would-be
Goldstone bosons. For simplicity, we assume a $CP$ symmetry in the scalar
sector so that the two $CP$-even states, $h_{1}$ and $h_{2}$, do not mix with
the $CP$-odd state $A$. Both $h_{1}$ and $h_{2}$ can mix through

$\begin{pmatrix}h_{1}\\\
h_{2}\end{pmatrix}=\begin{pmatrix}c_{\alpha}&s_{\alpha}\\\
-s_{\alpha}&c_{\alpha}\end{pmatrix}\begin{pmatrix}h\\\ H\end{pmatrix},$ (2)

where $c_{\alpha}$ ($s_{\alpha}$) stands for $\cos\alpha$ ($\sin\alpha$). Here
we identify $h$ with the 125-GeV Higgs boson. Since the properties of the $h$
agree with the SM predictions [25, 26, 27], the mixing angle $s_{\alpha}$ is
expected to be small. The masses of the extra Higgs bosons $H$, $A$, and
$H^{+}$ are, in principle, arbitrary.

The Yukawa couplings of $H_{1}$ generate fermion masses while the Yukawa
couplings of $H_{2}$ give rise to potential flavor violations. The most
general form of the Yukawa couplings of $H_{1}$ and $H_{2}$ to the leptons is
given by

$\begin{split}\mathcal{L}_{yuk}&\supset-\frac{\sqrt{2}m_{i}}{v}\bar{L}_{i}\ell_{Ri}H_{1}-\sqrt{2}Y_{ij}\bar{L}_{i}\ell_{Rj}H_{2}+\text{h.c.},\end{split}$
(3)

where $L$ denotes the lepton doublet, $m_{i}$ is the $i$th generation charged
lepton mass, and $i,j=e,\mu,\tau$ indicate lepton generations.

In this work, we will focus on LFV in the $e$-$\mu$ sector, so we take only
$Y_{e\mu}$ and $Y_{\mu e}$ in Eq. (3) to be nonzero. The couplings $Y_{e\mu}$
and $Y_{\mu e}$, in principle, can be complex. However, the imaginary part of
$Y_{e\mu}Y_{\mu e}$ is strongly constrained by the electron electric dipole
moment measurement [28]. Therefore, for simplicity, we assume that both
$Y_{e\mu}$ and $Y_{\mu e}$ are real in this work.

In our minimal scenario, only $h$ and $H$ can be singly produced via gluon
fusion and vector boson fusion processes. This makes $H$ the only relevant
resonance for CMS and ATLAS LFV searches. The $H$ production cross-sections
and non-LFV partial decay widths can be obtained from the would-be SM Higgs
boson values by scaling them with a factor of $s_{\alpha}^{2}$. On the other
hand, the $h$ production cross-sections and non-LFV partial decay widths are
reduced by a factor of $c_{\alpha}^{2}$ from their SM values. In our analysis,
we use the SM-like Higgs cross-sections and decay widths provided by the LHC
Higgs Cross Section Working Group [29].

##  III Low-energy LFV constraints

Naturally, the LFV couplings $Y_{e\mu}$ and $Y_{\mu e}$, together with the
scalar mixing $s_{\alpha}$, will induce $\mu\to e\gamma$, which proceeds via
loops. The partial decay width is given by

$\Gamma(\mu\to
e\gamma)=\frac{\alpha_{em}m_{\mu}^{5}}{64\pi^{4}}\left(|c_{L}|^{2}+|c_{R}|^{2}\right),$
(4)

where $\alpha_{em}$ is the electromagnetic fine-structure constant. The Wilson
coefficients $c_{L,R}$ arise at one- and two-loop level. The one-loop
contributions are given by

$\displaystyle c_{L}^{(1)}=-\frac{s_{2\alpha}}{24}\frac{m_{\mu}Y_{\mu e}}{v}$
$\displaystyle\left[\frac{1}{m_{h}^{2}}\left(4+3\ln\frac{m_{\mu}^{2}}{m_{h}^{2}}\right)\right.$
$\displaystyle\quad\left.-\frac{1}{m_{H}^{2}}\left(4+3\ln\frac{m_{\mu}^{2}}{m_{H}^{2}}\right)\right],$
(5)

with $s_{2\alpha}\equiv\sin 2\alpha$. The main contributions to the Wilson
coefficients $c_{L,R}$ are dominated by two-loop photon-exchange diagrams
involving the $W$ boson and top quark. They are given by

$\displaystyle c_{L}^{(2W)}$
$\displaystyle\simeq-\frac{\alpha_{em}s_{2\alpha}Y_{\mu e}}{8\pi
vm_{\mu}}\left[\frac{f(z_{Wh})-g(z_{Wh})}{2z_{Wh}}+3f(z_{Wh})\right.$
$\displaystyle\quad\left.+\frac{23}{4}g(z_{Wh})+\frac{3}{4}h(z_{Wh})-(h\to
H)\right]$ (6) $\displaystyle c_{L}^{(2t)}$
$\displaystyle\simeq\frac{\alpha_{em}s_{2\alpha}Y_{\mu e}}{3\pi
vm_{\mu}}\left[f(z_{th})-f(z_{tH})\right],$ (7)

with $z_{ab}=m_{a}^{2}/m_{b}^{2}$. The loop functions $f$, $g$, and $h$ can be
found in Ref. [30]. Note the cancellation between the $h$ and the $H$
contributions due to the mixing of Eq. (2). Such a cancellation weakens the
$\mu\to e\gamma$ constraint as $m_{H}$ approaches $m_{h}$. It should be noted
that the loop functions in Eqs. (5)–(7) are typically $\mathcal{O}(1)$. Hence,
the one-loop coefficient is parametrically suppressed by
$\,m_{\mu}^{2}/(\alpha_{em}m_{h(H)}^{2})$ compared to the two-loop ones. The
Wilson coefficient $c_{R}$ is obtained from $c_{L}$ by replacing $Y_{\mu e}\to
Y_{e\mu}$.

In addition to such $W$ and top contributions, we have included contributions
from the so-called “set C” diagrams [30]. Numerically, they give a correction
of order 10% to the Wilson coefficients. Other contributions, i.e., two-loop
$Z$-exchange diagrams are found to be more suppressed compared to the others.

The most stringent constraint on the $\mu\to e\gamma$ decay is provided by the
MEG experiment with $\text{BR}(\mu\to e\gamma)\leq 4.2\times 10^{-13}$ [31].
The upgraded MEGII experiment, which is currently taking data, is expected to
push the bound down to $6\times 10^{-14}$ [32].

The Wilson coefficients $c_{L}$ and $c_{R}$ also lead to $\mu\to e$ conversion
in atomic nuclei. In addition, the $\mu\to e$ conversion also gets tree-level
contributions mediated by the $CP$-even Higgs bosons. In our scenario, the
conversion rate is given by [33] (see also [14])

$\Gamma(\mu\to
e)=m_{\mu}^{5}\left|\frac{e\,c_{R}}{16\pi^{2}}D+\sum_{N=p,n}g_{L}^{N}S^{N}\right|^{2}+(L\leftrightarrow
R),$ (8)

where $D$ and $S^{N}$ are overlap integrals, whose values are given in [33].
The effective coupling $g_{L}^{N}$ is given by

$g_{L}^{N}=s_{2\alpha}Y_{e\mu}\left(\frac{1}{m_{H}^{2}}-\frac{1}{m_{h}^{2}}\right)\frac{m_{N}}{v}\sum_{q}f^{(q,N)},$
(9)

where $f^{(q,N)}$ and $m_{N}$ denote the nucleon form factor and the nucleon
mass, respectively. The form factors are given in Refs. [34, 35], and in above
equation, they are summed over all quark flavors. The effective coupling
$g_{R}^{N}$ can be obtained from $g_{L}^{N}$ by replacing $Y_{e\mu}\to Y_{\mu
e}$. As in the $\mu\to e\gamma$ case, the $\mu\to e$ conversion constraint
also suffers from the same blind spot when $m_{H}$ approaches $m_{h}$.
Currently, the strongest constraint comes from conversion in the gold nucleus,
which sets $\Gamma(\mu\to e$)/$\Gamma$(captured) $<7\times 10^{-13}$ [36].

##  IV Collider Searches for LFV decays

The LFV couplings $Y_{e\mu}$ and $Y_{\mu e}$, together with the mixing angle
$s_{\alpha}$, lead to LFV decays of the $h$ and $H$. The partial decay width
$h\to e^{\pm}\mu^{\mp}$ is given by

$\Gamma(h\to
e^{\pm}\mu^{\mp})=\frac{s_{\alpha}^{2}m_{h}}{8\pi}\left(\left|Y_{e\mu}\right|^{2}+\left|Y_{\mu
e}\right|^{2}\right).$ (10)

The partial decay width for $H\to e^{\pm}\mu^{\mp}$ can be obtained from Eq.
(10) by changing $s_{\alpha}\to c_{\alpha}$ and $m_{h}\to m_{H}$.

Both ATLAS and CMS Collaborations have searched for such LFV decay of $h$
without any positive signal [37, 10]. Assuming the SM production cross-
sections for the $h$, the ATLAS has set an upper bound $\text{BR}({h\to
e^{\pm}\mu^{\mp}})\leq 6.2\times 10^{-5}$, while the CMS has set a slightly
stronger constraint $\text{BR}({h\to e^{\pm}\mu^{\mp}})\leq 4.4\times
10^{-5}$. As mentioned earlier, in our scenario, the $h$ production cross-
sections get corrected by $c_{\alpha}^{2}$, so does the bound on the $h\to
e^{\pm}\mu^{\mp}$ branching ratio. Hence, the CMS upper bound now reads
$c_{\alpha}^{2}\text{BR}({h\to e^{\pm}\mu^{\mp}})\leq 4.4\times 10^{-5}$.

The collider search for $H\to e^{\pm}\mu^{\mp}$ was done by the CMS
Collaboration [37]. The result was presented in terms of the $H$ production
cross-section times the branching ratio into $e^{\pm}\mu^{\mp}$. In our model,
it is formulated as

$\sigma\times\text{BR}({H\to
e^{\pm}\mu^{\mp}})=\frac{s_{\alpha}^{2}\sigma_{SM}\Gamma(H\to
e^{\pm}\mu^{\mp})}{s_{\alpha}^{2}\Gamma_{SM}+\Gamma(H\to e^{\pm}\mu^{\mp})},$
(11)

where $\sigma_{SM}$ and $\Gamma_{SM}$ are the cross-section and the total
decay width of the SM-like Higgs with mass $m_{H}$.

###  IV.1 The CMS 146-GeV excess

The CMS search for a resonance decaying to $e^{\pm}\mu^{\mp}$ has observed an
excess of events around 146 GeV with a local significance of 3.8$\sigma$ [37].
However, the corresponding ATLAS search does not find any excess at 146 GeV.
By performing a simple profiled likelihood ratio analysis, we find that the
CMS signal is disfavored by the ATLAS measurement at the 2.5$\sigma$ level.

Taking both CMS and ATLAS results into consideration, we perform a naive
combination of the two searches using the counting experiment. The combined
tool with asymptotic approximation [38] is used to calculate the best fit
cross-section and its significance. The data and background Monte Carlo (MC)
from both experiments, as well as the CMS signal MC for $m_{H}=146$ GeV, are
extracted from the Data-MC plots in Refs. [37, 10]. In the case of ATLAS MC,
we estimate the signal model using Madgraph5 [39] to generate parton-level
productions and $H$ decays, followed by showering and hadronization with
Pythia8 [40]. Finally, the detector responses are simulated using Delphes3
[41]. The overall uncertainty is estimated to be 1%, as mentioned in Ref. [10]
for the $e\mu$ channel. We find that the combination reduces the local
significance of the excess to 3.3$\sigma$, with the best fit cross section
$\sigma\times\text{BR}({H\to e^{\pm}\mu^{\mp}})=2.92^{+0.91}_{-0.89}$ fb.

Figure 1: The region in the
$\sin\alpha$–$\sqrt{\left|Y_{e\mu}\right|^{2}+\left|Y_{\mu e}\right|^{2}}$
plane compatible with the CMS excess (red) and the CMS-ATLAS combination
(cyan) at the 1$\sigma$ level. The region surrounded by the red (cyan) dashed
lines is excluded by CMS multilepton search unless $m_{A}\gtrsim 800$ GeV. The
yellow, magenta, and blue lines are the upper bounds from the $h\rightarrow
e^{\pm}\mu^{\mp}$, $\mu\to e$ conversion, and $\mu\rightarrow e\gamma$
searches, respectively.

The region in $s_{\alpha}$–$\sqrt{Y_{e\mu}^{2}+Y_{\mu e}^{2}}$ parameter space
consistent with the CMS excess, together with the combined CMS and ATLAS data,
at 1$\sigma$ level is shown in Fig. 1. The preferred regions of parameter
space are compared against constraints from the $h\to e\mu$, $\mu\to e\gamma$,
and $\mu\to e$ conversion searches. Note that the small 1$\sigma$ region
preferred by the CMS excess consistent with the $\mu\to e\gamma$ bound is
centered around $s_{\alpha}=0.014$ and $\sqrt{Y_{e\mu}^{2}+Y_{\mu
e}^{2}}=6.6\times 10^{-4}$. In that figure, we also show the projected
sensitivity of the MEGII experiment. Should the origin of the CMS excess be a
new particle, the MEGII experiment will observe the $\mu\to e\gamma$ decay.

In addition to the aforementioned constraints, the 146-GeV excess is also
constrained by the CMS search for a new resonance decaying into a pair of
$ZZ^{*}$ [42] and a pair of tau leptons [43]. In both cases, the bounds are
$\mathcal{O}(100)$ fb, which correspond to $s_{\alpha}\lesssim 0.2$.

The 146-GeV excess is also constrained by multilepton searches via $HA$ pair
production. The $A$ will decay predominantly into $e\mu$, $hZ$, or $HZ$ if the
latter is kinematically open. A significant portion of the $H$, on the other
hand, will decay into $e\mu$ for sufficiently small $s_{\alpha}$. If this is
the case, the pair produced $HA$ will result in four-lepton signatures, which
are strongly constrained by the CMS multilepton searches [44]. From our
analysis, we find that the CMS multilepton constraints imply either
$m_{A}\gtrsim 800$ GeV or $s_{\alpha}\gtrsim 0.025$ provided that $A\to HZ$ is
kinematically open.

###  IV.2 Bounds for other masses

While the excess at 146 GeV seen by CMS could be a result of a new particle,
we should not draw a definite conclusion on its nature without more data.
Thus, we also study the possibility that the CMS search finds no excess
anywhere in the search region 110 GeV $\leq m_{H}\leq$ 160 GeV. In this
scenario, we take the 95% confidence level upper limits on
$\sigma\times\text{BR}({H\to e\mu})$ reported by CMS and compare them against
constraints from $\mu\to e\gamma$ and $\mu\to e$ conversion.

Figure 2: The comparison between CMS observed limit (solid red), CMS expected
limit (dashed red), ATLAS observed limit (black), current $\mu\to e\gamma$
constraint (solid blue), current $\mu\to e$ conversion (magenta), and the
projected $\mu\to e\gamma$ limit from MEGII (dotted blue) on the CMS search
region. The bounds are projected onto the
$m_{H}-s_{2\alpha}\sqrt{\left|Y_{e\mu}\right|^{2}+\left|Y_{\mu e}\right|^{2}}$
plane. For the CMS bounds, only the lower end of the projection is shown.

The bounds from collider search, $\mu\to e\gamma$, and $\mu\to e$ conversion
experiments are functions of the mixing angle $s_{\alpha}$ and the LFV
couplings $\sqrt{\left|Y_{e\mu}\right|^{2}+\left|Y_{\mu e}\right|^{2}}$. While
the $\mu\to e\gamma$ and $\mu\to e$ conversion constraints depend on the
product $s_{2\alpha}\sqrt{\left|Y_{e\mu}\right|^{2}+\left|Y_{\mu
e}\right|^{2}}$, the collider constraint is a more complicated function. For
illustrative purposes, the collider bound, for certain value of $H$ mass,
would trace out a curve in the
$s_{\alpha}$–$\sqrt{\left|Y_{e\mu}\right|^{2}+\left|Y_{\mu e}\right|^{2}}$
plane similar to the preferred region for the 146-GeV excess in Fig. 1.

In order to compare these two types of constraints on the entire CMS search
region, we project the constraints onto the
$s_{2\alpha}\sqrt{\left|Y_{e\mu}\right|^{2}+\left|Y_{\mu e}\right|^{2}}$ line.
Such a projection maps the $\mu\to e\gamma$ and $\mu\to e$ conversion
constraints to a point. On the other hand, the bound from CMS search gets
mapped into a bounded-from-below interval. If the minimum of such an interval
is smaller than the $\mu\to e\gamma$ ($\mu\to e$ conversion) projection, it is
said that the collider bound is more constraining than the $\mu\to e\gamma$
($\mu\to e$ conversion) bound for some part of the
$s_{\alpha}$–$\sqrt{\left|Y_{e\mu}\right|^{2}+\left|Y_{\mu e}\right|^{2}}$
parameter space. In Fig. 2 we show the comparison between the collider bounds
and the $\mu\to e\gamma$ and $\mu\to e$ conversion constraints for 110 GeV
$\leq m_{H}\leq$160 GeV. Note that, in Fig. 2, only the minimum of the
collider constraint projection is shown for each $m_{H}$ value. As can be seen
from the plot, for $m_{H}\lesssim$ 140 GeV, the collider search can be more
constraining than even the projected MEGII bound. Hence, we advocate both the
ATLAS and CMS Collaborations to continue searching for $H\rightarrow
e^{\pm}\mu^{\mp}$ in that low-$m_{H}$ region.

##  V Conclusion and Discussion

In this work, we have analyzed the LHC searches for a new resonance decaying
to $e^{\pm}\mu^{\mp}$ in the context of the type-III 2HDM. This analysis is
motivated by a recent CMS search [37] that explores the LFV decay of a new
scalar boson between 110 GeV $<m_{H}<160$ GeV. Within the search region, the
CMS finds a possible excess at $m_{H}=146$ GeV with a local (global)
significance of $3.8\sigma$ $(2.8\sigma)$. A simplistic combination of the CMS
and ATLAS searches reduces the local significance to 3.3$\sigma$, with
$\sigma\times\text{BR}({H\to e^{\pm}\mu^{\mp}})=2.92^{+0.91}_{-0.89}$ fb. In
the type-III 2HDM context, the 146-GeV excess is only moderately constrained
by the current $\mu\to e\gamma$ data, while the future MEGII search will probe
the whole parameter region preferred by the excess.

It is interesting to note that the 146-GeV excess could be related to another
excess from the so-called “multilepton anomalies” [45]. We leave this
investigation for possible future work.

In the event that the excess is due to an upward fluctuation in the data, we
analyze the bounds on $\sigma\times\text{BR}({H\to e^{\pm}\mu^{\mp}})$
provided by CMS over the whole search region. The comparison between the CMS
bounds and the low-energy counterparts is shown in Fig. 2. From the plot, we
see that for 110 GeV $\leq m_{H}\lesssim$ 140 GeV, the current CMS bound is
better than the current and even the projected MEGII bounds for $\mu\to
e\gamma$. Thus, we encourage both CMS and ATLAS Collaborations to continue
searching for LFV decays of a new resonance in this low-mass region.

Lastly, we note that even though we take $Y_{ee}$ and $Y_{\mu\mu}$ to be zero
in this work, it does not mean we need a large hierarchy between the flavor-
violating and the flavor-conserving couplings. In fact, $Y_{ee/\mu\mu}$ can be
comparable to $Y_{e\mu/\mu e}$ and still be consistent with low-energy
constraints. We set them to zero in the spirit of minimal scenario.

###### Acknowledgements.

R. P. was supported by the Parahyangan Catholic University under grant no.
III/LPPM/2023-02/32-P. J. J. was supported in part by the Indonesia Toray
Science Foundation. P. U. was supported in part by the Mid-Career Research
Grant from the National Research Council of Thailand under contract no.
N42A650378. N. S. was supported by Thailand NSRF via PMU-B under grant number
B37G660013. R. P. thanks The Abdus Salam International Centre for Theoretical
Physics for kind hospitality while this work was being completed.

## References

  * Aad  _et al._ [2012] G. Aad _et al._ (ATLAS), Phys. Lett. B 716, 1 (2012), arXiv:1207.7214 [hep-ex] .
  * Chatrchyan  _et al._ [2012] S. Chatrchyan _et al._ (CMS), Phys. Lett. B 716, 30 (2012), arXiv:1207.7235 [hep-ex] .
  * Tumasyan  _et al._ [2022a] A. Tumasyan _et al._ (CMS), Nature 607, 60 (2022a), arXiv:2207.00043 [hep-ex] .
  * ATLAS Collaboration [2022] ATLAS Collaboration, Nature 607, 52 (2022), [Erratum: Nature 612, E24 (2022)], arXiv:2207.00092 [hep-ex] .
  * Aaboud  _et al._ [2018] M. Aaboud _et al._ (ATLAS), Phys. Lett. B 786, 59 (2018), arXiv:1808.08238 [hep-ex] .
  * Sirunyan  _et al._ [2018a] A. M. Sirunyan _et al._ (CMS), Phys. Rev. Lett. 121, 121801 (2018a), arXiv:1808.08242 [hep-ex] .
  * Aad  _et al._ [2016a] G. Aad _et al._ (ATLAS, CMS), JHEP 08, 045 (2016a), arXiv:1606.02266 [hep-ex] .
  * Sirunyan  _et al._ [2021] A. M. Sirunyan _et al._ (CMS), JHEP 01, 148 (2021), arXiv:2009.04363 [hep-ex] .
  * CMS Collaboration [2022a] CMS Collaboration, (2022a), arXiv:2208.00265 [hep-ex] .
  * Aad  _et al._ [2020a] G. Aad _et al._ (ATLAS), Phys. Lett. B 801, 135148 (2020a), arXiv:1909.10235 [hep-ex] .
  * Aad  _et al._ [2022] G. Aad _et al._ (ATLAS), Eur. Phys. J. C 82, 717 (2022), arXiv:2201.11428 [hep-ex] .
  * CMS Collaboration [2022b] CMS Collaboration, (2022b), arXiv:2205.05550 [hep-ex] .
  * Blankenburg  _et al._ [2012] G. Blankenburg, J. Ellis, and G. Isidori, Phys. Lett. B 712, 386 (2012), arXiv:1202.5704 [hep-ph] .
  * Harnik  _et al._ [2013] R. Harnik, J. Kopp, and J. Zupan, JHEP 03, 026 (2013), arXiv:1209.1397 [hep-ph] .
  * Davidson and Grenier [2010] S. Davidson and G. J. Grenier, Phys. Rev. D 81, 095016 (2010), arXiv:1001.0434 [hep-ph] .
  * Crivellin  _et al._ [2013] A. Crivellin, A. Kokulu, and C. Greub, Phys. Rev. D 87, 094031 (2013), arXiv:1303.5877 [hep-ph] .
  * Kopp and Nardecchia [2014] J. Kopp and M. Nardecchia, JHEP 10, 156 (2014), arXiv:1406.5303 [hep-ph] .
  * Buschmann  _et al._ [2016] M. Buschmann, J. Kopp, J. Liu, and X.-P. Wang, JHEP 06, 149 (2016), arXiv:1601.02616 [hep-ph] .
  * Primulando and Uttayarat [2017] R. Primulando and P. Uttayarat, JHEP 05, 055 (2017), arXiv:1612.01644 [hep-ph] .
  * Altmannshofer  _et al._ [2016] W. Altmannshofer, J. Eby, S. Gori, M. Lotito, M. Martone, and D. Tuckler, Phys. Rev. D 94, 115032 (2016), arXiv:1610.02398 [hep-ph] .
  * Primulando  _et al._ [2020] R. Primulando, J. Julio, and P. Uttayarat, Phys. Rev. D 101, 055021 (2020), arXiv:1912.08533 [hep-ph] .
  * Barman  _et al._ [2023] R. K. Barman, P. S. B. Dev, and A. Thapa, Phys. Rev. D 107, 075018 (2023), arXiv:2210.16287 [hep-ph] .
  * Hayrapetyan  _et al._ [2023] A. Hayrapetyan _et al._ (CMS), (2023), arXiv:2305.18106 [hep-ex] .
  * Georgi and Nanopoulos [1979] H. Georgi and D. V. Nanopoulos, Phys. Lett. B 82, 95 (1979).
  * Aad  _et al._ [2016b] G. Aad _et al._ (ATLAS, CMS), JHEP 08, 045 (2016b), arXiv:1606.02266 [hep-ex] .
  * Sirunyan  _et al._ [2019] A. M. Sirunyan _et al._ (CMS), Eur. Phys. J. C 79, 421 (2019), arXiv:1809.10733 [hep-ex] .
  * Aad  _et al._ [2020b] G. Aad _et al._ (ATLAS), Phys. Rev. D 101, 012002 (2020b), arXiv:1909.02845 [hep-ex] .
  * Andreev  _et al._ [2018] V. Andreev _et al._ (ACME), Nature 562, 355 (2018).
  * de Florian  _et al._ [2016] D. de Florian _et al._ (LHC Higgs Cross Section Working Group),  2/2017 (2016), 10.23731/CYRM-2017-002, arXiv:1610.07922 [hep-ph] .
  * Chang  _et al._ [1993] D. Chang, W. S. Hou, and W.-Y. Keung, Phys. Rev. D 48, 217 (1993), arXiv:hep-ph/9302267 .
  * Baldini  _et al._ [2016] A. M. Baldini _et al._ (MEG), Eur. Phys. J. C 76, 434 (2016), arXiv:1605.05081 [hep-ex] .
  * Baldini  _et al._ [2018] A. M. Baldini _et al._ (MEG II), Eur. Phys. J. C 78, 380 (2018), arXiv:1801.04688 [physics.ins-det] .
  * Kitano  _et al._ [2002] R. Kitano, M. Koike, and Y. Okada, Phys. Rev. D 66, 096002 (2002), [Erratum: Phys.Rev.D 76, 059902 (2007)], arXiv:hep-ph/0203110 .
  * Crivellin  _et al._ [2014] A. Crivellin, M. Hoferichter, and M. Procura, Phys. Rev. D 89, 093024 (2014), arXiv:1404.7134 [hep-ph] .
  * Bishara  _et al._ [2016] F. Bishara, J. Brod, P. Uttayarat, and J. Zupan, JHEP 01, 010 (2016), arXiv:1504.04022 [hep-ph] .
  * Bertl  _et al._ [2006] W. H. Bertl _et al._ (SINDRUM II), Eur. Phys. J. C 47, 337 (2006).
  * CMS Collaboration [2023] CMS Collaboration, (2023), CMS-PAS-HIG-22-022 .
  * Cowan  _et al._ [2011] G. Cowan, K. Cranmer, E. Gross, and O. Vitells, The European Physical Journal C 71 (2011), 10.1140/epjc/s10052-011-1554-0.
  * Alwall  _et al._ [2014] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, JHEP 07, 079 (2014), arXiv:1405.0301 [hep-ph] .
  * Bierlich  _et al._ [2022] C. Bierlich _et al._ ,  (2022), 10.21468/SciPostPhysCodeb.8, arXiv:2203.11601 [hep-ph] .
  * de Favereau  _et al._ [2014] J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaître, A. Mertens, and M. Selvaggi (DELPHES 3), JHEP 02, 057 (2014), arXiv:1307.6346 [hep-ex] .
  * Sirunyan  _et al._ [2018b] A. M. Sirunyan _et al._ (CMS), JHEP 06, 127 (2018b), [Erratum: JHEP 03, 128 (2019)], arXiv:1804.01939 [hep-ex] .
  * CMS Collaboration [2022c] CMS Collaboration, (2022c), arXiv:2208.02717 [hep-ex] .
  * Tumasyan  _et al._ [2022b] A. Tumasyan _et al._ (CMS), JHEP 04, 147 (2022b), arXiv:2106.14246 [hep-ex] .
  * Crivellin  _et al._ [2021] A. Crivellin, Y. Fang, O. Fischer, A. Kumar, M. Kumar, E. Malwa, B. Mellado, N. Rapheeha, X. Ruan, and Q. Sha, (2021), arXiv:2109.02650 [hep-ph] .