Riemann hypothesis: is Bender-Brody-Müller Hamiltonian a new line of attack? There is a beautiful paper in Physical Review Letters [PRL 118, 130201 (2017), DOI:10.1103/PhysRevLett.118.130201] by Carl Bender, Dorje Brody, and Markus Müller (BBM) on a Hamiltonian approach to the Riemann Hypothesis. The paper  is surprisingly easy to follow
for a physicist.
BBM define a Hamiltonian
$$\hat{H} = (1- e^{-i \hat{p}} )^{-1} ( \hat{x} \hat{p} + \hat{p} \hat{x} ) (1- e^{-i \hat{p}})$$
where $p=- i \partial_x$ is the momentum operator in $\hbar=1$ units.
The authors show that eigenfunctions of $\hat{H}$ vanishing at infinity must be be in the form of Hurwitz theta function, $\psi_z= - \zeta(z,x+1)$, so that $$\hat{H} \psi_z = i (2 z -1)\psi_z $$
Imposing a boundary condition $\psi_z(0)=0$, by the virtue of $\zeta(z,1)=\zeta(z)$, they show that all non-trivial zeros of the Riemann $\zeta$ function must be eigenvalues of $\hat{H}$ with the imposed boundary conditions.
BBM call this result a "a complex extended
version of the Berry-Keating conjecture" and go on to provide heuristic arguments that all eigenvalues of $\hat{H}$ are real.
How promising is this new development in the context of solving the Riemann Hypothesis?
Steven Strogatz seems optimistic.
Update [19.10.2020]: The authors of [1] have published additional comments [3
, 4], including a response [4] to @Jean Bellissard's comment [5] that grew out of his answer below.
 A: The paper in question is written at the physical level of rigor, is published in a physics journal, and doesn't seem to make claims that it establishes mathematically rigorous results. However, I'm not convinced by the criticism so far that it's impossible to turn the main idea of the paper into a mathematically well-defined $L^2$ spectral problem whose solution would imply RH.
Perhaps something along the following lines may work.
1) The eigenfunctions $\psi_z(x)$ of $\hat{H}$ are indeed not in $L^2(0, \infty)$ but it seems that one can map the ones that correspond to the non-trivial zeros of $\zeta$ as follows. Let $w(x)$ be a positive smooth function on $(0,\infty)$ that coincides with $1/x$ near $0$ and with $1/x^{\frac{3}{2}}$ near infinity and let $\hat{W}$ be the operator given by multiplication with $w$. Then $\hat{W}\psi_z \in L^2(0, \infty)$ if and only if $z$ is a non-trivial zero of the Riemann zeta function. Indeed, $\hat{W}\psi_z$ is integrable at $0$ if and only if $\psi_z(0)=0$, i.e. if $z$ is a zero of $\zeta$. Further, known asymptotic properties of the Hurwitz zeta function (already mentioned in the BBM paper, namely sublinear growth for $z$ in the critical strip and much faster growth when $z$ is a trivial zero) imply that $\hat{W}\psi_z$ is square-integrable at infinity only when the zero $z$ is non-trivial.
(The behavior of $\psi_z(x)$ at infinity is also  discussed in Jean Bellissard's  Comment paper  arxiv.org/abs/1704.02644, using the integral representation of the Hurwitz zeta function.)
2) Now consider the operator $\hat{H}_W=\hat{W}\hat{H}\hat{W}^{-1}$, so that, at least formally, $\hat{W}\psi_z$ are eigenfunctions of $\hat{H}_W$ for every non-trivial zero $z$. To make this rigorous, one has to show that $\hat{H}_W$ is a densely defined closable operator on $L^2(0,\infty)$ and that $\hat{W}\psi_z$ are in the domain of its closure. (That of course would include settling quite a few technicalities regarding the proper definition and domain of $\hat{\Delta}^{-1}$, the inverse of the difference operator considered in the BBM paper, but I don't see why this would be apriori impossible.)
3) Assuming that was described in 2) can be done, one gets an embedding of the nontrivial zeros into the discrete spectrum of $\hat{H}_W$. Now It's known that the Berry-Keating Hamiltonian $\hat{x} \hat{p} + \hat{p} \hat{x}$ is selfadjoint and hence has real spectrum, so one may try to use the fact that $\hat{H}_W$ is similar to $\hat{x} \hat{p} + \hat{p} \hat{x}$ to obtain estimates on the resolvent operator that imply that the spectrum of $\hat{H}_W$ is contained in the spectrum of $\hat{x} \hat{p} + \hat{p} \hat{x}$.
Finally, of course it may be that  $\hat{H}_W$ has lots of other spectrum including complex one, while its discrete spectrum is real, so that RH may still be true but not accessible via this line of attack.
A: I tried to put my own argument on this blog recently but I failed. Let me try again.
1)- The original definition of $\zeta(z,x) $ is 
$$\zeta(z,x) = \sum_{n=1}^\infty \frac{1}{(n+x)^z}
$$
Then $\zeta(z,x)-\zeta(z,x-1)=-1/x^z$. Then, using the integral representation, this relation extends by analiticity to the maximal domain in which both sides are defined and holomorphic. For our purpose, it is enough to consider $x>0$ and $\Re{z}>0$. Let $\zeta_z$ denote the map $x\to \zeta(z,x)$. Then, using the generator $A$ of the dilation operator, it is easy to check that the eigenvalues have the form $1/2+\imath E$. For $A=(pq+qp)/2=-\imath xd/dx-\imath 1/2$. This gives $A(1/x^z)=\imath (z-1/2)(1/x^z)$. Going from $A$ to $H$ is done using the operator $\Delta$, leading to the formal equation $H\zeta_{1/2+\imath E}=E\zeta_z$. The Dirichlet b.c. at the origin forces E to satisfy $\zeta(1/2+\imath E)=0$.
2)- The first problem I see is that $\zeta_z$ does not seem to belong to the Hilbert space ${\mathcal H}=L^2(0,\infty)$. Using a classical method by Hadamard 
$$\int_0^\infty |\zeta(z,x)|^2 dx=\frac{1}{1-2\Re{z}}\zeta(2\Re{z}-1)$$
which converges for $\Re{z}>1$. Using the integral representation, a similar firmula can be obtained, but I failed to show that square integrability hold for $\Re{z}=1/2$. Hence I see a problem here.
3)- The other problem is the definition of $\Delta$. Let $S$ denote the translation operator, defined non rigorously by $Sf(x)=f(x-1)$. If restricted to $\mathcal H$, it is not unitary, because it is defined only for $x>1$. We can defined it by imposing $Sf(x)=0$ for $x\in [0,1]$. Then it is only a partial isometry. For $S^\ast f(x)=f(x+1)$ leading to $S^\ast S=I$, but $S S^\ast =P$ is the projection onto $L^2(1,\infty)$. Using these notations $\Delta= I-S=(S^\ast-I)S$. But we get a problem here: the function $\zeta_z$ is defined for $x>-1$, and the extension to the interval $(-1,0]$ is explicitly used in 1)-. So we cannot use 
$S$. But then what is the operator $e^{-\imath \hat{p}}$ used by the authors? 
4)- If $\hat{p}$ is the usual operator
$$\hat{p}=-\imath \frac{d}{dx}$$
then there is a problem with its domain of definition. On $L^2(\mathbb R)$, it is selfadjoint as can be seen by using Fourier transform. But on $\mathcal H$, it is not. This is a classical exercise found in the very old book of Courant-Hilbert. Namely one can always define it in the set of $L^2$ functions with $L^2$ derivative, vanishing at $x=0$. Then it is symmetric. If so its adjoint is defined on the same space but without the vanishing at $x=0$. Not only the adjoint is not symmetric but its set of eigenvalues is the open lower half plane. This is because if $f_z(x)=e^{zx}$, then $\hat{p}^\ast f_z=-\imath z f_z$, while $f_z\in \mathcal{H}$ for $\Re{z}<0$. The same argument shows that $+\imath$ cannot be an eigenvalue. This means that the "defect indices", namely the dimension of the eigenspaces with eigenvalues $\pm\imath$, are not equal. Then, the von Neumann theorem show that the operator $\hat{p}$ has no selfajoint extension. If not selfadjoint, the definition of its exponential becomes a problem, because the functional calculus is nit defined in general. 
5)- The previous argument can be rephrased in terms of the operator $S$. Its adjoint admits a lot of eigenvalues, namely the points inside the unit disk. 
In conclusion, the sloppyness of the definitions used but the authors leads to a complete mess. Nothing is correct in this paper. 
As long as physicists use algebra, or algorithmic arguments, they can find outstanding results. But when it comes to analysis, they may loose their judgment, and grave mistakes show up at the corner. Analysis is not easily amenable to algoritmic descriptions. And this is precisely where the power of Mathematics lies: by manipulating infinities, Mathematics goes way beyond the Church-Turing definition of computability. And what is Analysis if not manipulating infinities, through limits, convergence and the likes?
A: I didn't get everything written in this paper.
But keep in mind this line of attack has some good chances to be way too simple : all this works the same for $F(s) = \alpha L(s,\chi_5)+\overline{\alpha}L(s,\overline{\chi_5})= 2\sum\limits_{n=1}^\infty \Re(\alpha \;\chi_5(n)) n^{-s}$ where $\chi_5$ is the non-real character modulo $5$ and $\alpha \in \mathbb{C}$. 
$F(s)$ has the same kind of  integral representation and functional equation as $\zeta(s)$, so we can write for it the same kind of differential operator. But the RH obviously fails for $F(s)$ (it doesn't have an Euler product)
A: In my previous answer I left out the issue regarding how one can make sense rigorously of the operators $\hat{\Delta}=1-e^{-i \hat{p}}$ and $\hat{\Delta}^{-1}=(1- e^{-i \hat{p}} )^{-1}$ used in the BBM paper. I'd like to elaborate on this a bit. 
1) As it's well-known (and explained in Jean Bellissard's answer) the momentum operator $\hat{p}$ doesn't admit selfadjoint extensions on $L^2(0,\infty)$, hence one cannot use (at least not directly) functional calculus to make sense of $\hat{\Delta}$. In the BBM paper $\hat{\Delta}$ is interpreted as a difference operator and $\hat{\Delta}^{-1}$ as its inverse defined via infinite series (apparently, following Euler rather closely), but this approach is not satisfactory for many reasons.
2) It would be much more fitting to the general quantization idea to think of $\hat{\Delta}$ and $\hat{\Delta}^{-1}$ as pseudodifferential operators. Since this cannot be done directly, I suggest to proceed as follows, using appropriate cut-off functions to obtain smooth compactly supported symbols. Let $\kappa_{\varepsilon}(p)$ be a family of functions in $C^{\infty}_{c}(\mathbb{R})$ which is 1 on  $[-1/ \varepsilon, 1/ \varepsilon]$ and $0$ away from this interval, so that $\kappa_{\varepsilon} \rightarrow 1 $ as $\varepsilon \rightarrow 0$. Let $\hat{\Delta}_{\varepsilon}$ be the family of pseudodifferential (in fact smoothing) operators with symbols $(1-e^{-i p})\kappa_{\varepsilon}(p)$. Further, let $\chi_{\varepsilon}(p)$ be another family in $C^{\infty}_{c}(\mathbb{R})$  which is $0$ on $[-\varepsilon, \varepsilon]$ as well as outside $[-1/ \varepsilon, 1/ \varepsilon]$ and 1 in between, so that still one has $\chi_{\varepsilon} \rightarrow 1 $ as $\varepsilon \rightarrow 0$. Let $\hat{\Delta}_{\varepsilon}^{-1}$  be the family of pseudodifferential operators with symbols $(1-e^{-i p})^{-1}\chi_{\varepsilon}(p)$. (There is a plenty of research on pseudodifferential operators on the half-line and I suspect some of it may be very relevant to this particular context.)
3) Now an appropriate limit (say, weak graph limit) as $\varepsilon \rightarrow 0$ of these families of operators may or may not exist but even if doesn't not exist, it looks plausible to me that one can show that the formal relations among $\hat{\Delta}$, $\hat{\Delta}^{-1}$  and the Hurwitz zeta function given in the BBM paper hold ''asymptotically'' as $\varepsilon \rightarrow 0$. Thus the same should be true for the formal eigenvalue/eigenfunction equation, so that if one defines $\hat{H}_{\varepsilon}=\hat{\Delta}_{\varepsilon}^{-1}(\hat{x} \hat{p} + \hat{p} \hat{x})\hat{\Delta}_{\varepsilon}$, one would have
$$\lim_{\varepsilon \rightarrow 0}\hat{H}_{\varepsilon}\psi_z=i(2z-1)\psi_z  $$
as, say, weak limit of distributions.
4) Finally, one can look at the operators $\hat{H}_{\varepsilon,W}=\hat{W}\hat{H}_{\varepsilon}\hat{W}^{-1}$, where $\hat{W}$ was defined in my previous answer. The desired reality of the eigenvalues would then follow from a uniform (in $\varepsilon$) bound on the norms of resolvent operators of $\hat{H}_{\varepsilon,W}$ for every fixed non-real $z$. There are well-known sufficient conditions for the $L^2$-boundedness of a pseudodifferential operator in terms of its symbol.
A: It's an interesting article.   The main problem I can see is that their boundary condition invokes the zeta function itself,  so this approach may be incorporating the RH itself from the very beginning.   The approach doesn't seem to bring in anything deep about the primes nor number theory.  
Having worked on this problem over several years,  I remain skeptical about the Hilbert-Polya idea.   Berry-Keating did some very interesting and important work.   However the random matrix theory connection would seem to cast doubts on Hilbert-Polya rather than support it:  if the statistics of the zeros can be explained by  random matrix theory,  it is hard to imagine a simple,  non-random  hamiltonian that would reproduce this spectrum.   But this is not a rigorous counter-argument of course.  
A: I am a number theorist who has thought about the Riemann Hypothesis for many years.  This paper has elementary mistakes and there are simple counterexamples to the approach taken.  It would not have been published in a mathematics journal, and I am disappointed (but not surprised) it has been published in PRL.  Bellissard has listed some of the mistakes in the analysis (see arXiv:1704.02644) and the number theory follies would be clear to any graduate student in the field. 
It is interesting sociologically that so many well known physicists are blogging or tweeting about this when they could just consult experts in their universities.  For those experts, this is a joke.  
A dark day for science.
A: The current situation is as follows. The zetagrid project no longer exists. The Z function has a pole in Re (z) = 1, but is analytic for Re (z) > 1, 
and divergent for Re (z) < 1. Riemann analytically extended the function so that it could be analytic in the region 0 < Re (z) < 1, which is what is known 
as the critical strip or research interval. Where we have that for Re (z) = 1/2 there are infinite zeros (Hardy, 1914) and it is assumed (RH) that all the zeros 
of the Z are there. However, in November of 2017, this server showed that there are infinite zeros of the Z in the critical interval distributed in 
different families of lines, with which the Riemann Hypothesis is false! In addition, the existence of the zeroes by reflection predicted by the 
Riemann functional equation was confirmed and according to which if in the interval 0 < Re (z) < 1/2 there is some zero, this zero must be reflected 
in the interval 1/2 < Re (z) < 1
enter image description here
0.44 + 1977.19-> 0.56 + 1977.19
0.45 + 1329.09-> 0.55 + 1329.09
0.45 + 1415.59-> 0.55 + 1415.59
0.46 + 25.0-> 0.54 + 25.0
0.46 + 37.6-> 0.54 + 37.6
0.47 + 14.09-> 0.53 + 14.09
0.47 + 21.0-> 0.53 + 21.0
0.48 + 14.09-> 0.52 + 14.09
0.48 + 21.0-> 0.52 + 21.0
0.49 + 14.09-> 0.51 + 14.09
0.49 + 21.0-> 0.51 + 21.0
all of them are zeros of the Z (by reflection), which also are not on the line Re (z) = 1/2, as it is easy to see and check using my expressions 
for the Z function. I leave that task to you.
A: Some of the main discrepancies regarding the Bender-Brody-Müller (BBM) conjecture are with regard to the domain of the operator, and the convergence of the eigenspectrum.  By considering the problem as a Schrödinger operator, one may obtain a convergent eigenspectrum.  Such solutions are discussed here, and here.
A: Just a thought I'd like to elaborate here on, inspired by @AndréLeClair's answer (its too long for a comment). 
The basic question is: What is the connection between random matrix theory (where some basic statistics of the zeros of the zeta function are recovered) and the (barely) $\mathcal{PT}$ symmetric (rather then self-adjoint) Bender-Brody-Müller (BBM) operator?
It seems, ideally random matrix theory would be exactly the Heisenberg matrix formulation of the BBM "Schrödinger type" formulation. 
Can that possibly be true considering the BBM operator looks so $non$-$random$?
Maybe: The Heisenberg uncertainty principle tells us that no quantum object can have a defined momentum and position in the same state. Measuring momentum and position one after the other leads inevitably to a result for which only probability distributions rather than defined values can be given. 
Would it make sense to interpret the Berry-Keating (or therefrom derived) operators containing products $\hat{p}\hat{x}$ such that one would force the system into eigenstates with localised (defined) product of momentum and position? One then would have to read these terms like $\hat{p}(\hat{x}(\psi))$ Leading to kind of a paradoxical state, which reflects itself already in that these operators are not self-adjoint. 
(I am not sure if this argument has already been considered by BBM).
A: We can prove Riemann hypothesis by using the operator $\hat{G} + \hat{G}^\dagger$, defined as
$$\hat{G} = \hat{p}\hat{\Delta}_-\hat{x}\hat{\Delta}_+$$
$$\hat{G}^\dagger = \hat{\Delta}_-\hat{x}\hat{\Delta}_+ \hat{p}$$
$\hat{\Delta}_+$ is the forward difference operator and $\hat{\Delta}_-$ is the backward difference operator.
If $r$ and $s$ are zeros of the Riemann zeta with real parts greater than $\frac{1}{2}$, we can show that $\hat{p}$ is self adjoint and that
$$\left< \varphi_r \left| \hat{G} \right| \varphi_s \right> = \left< \varphi_r \left| \hat{G}^\dagger \right| \varphi_s \right> $$
$\hat{G} + \hat{G}^\dagger$ is related to the BBM Hamiltonian $\hat{H} = \hat{\Delta}_+^{-1} \left( \hat{x} \hat{p} + \hat{p} \hat{x} \right) \hat{\Delta}_+$ by
$$ \hat{G} + \hat{G}^\dagger = \hat{\Delta}_- \hat{\Delta}_+ \hat{H} $$
Because $\varphi_s$ are eigenfunctions of $\hat{H}$ with eigenvalues $-i\hbar (1 - 2s)$, the expected value of $\hat{G} + \hat{G}^\dagger$ is
$$\left< \varphi_s \left| \hat{G} + \hat{G}^\dagger \right| \varphi_s \right>   = -i\hbar (1 - 2s) \left< \varphi_s \left| \hat{\Delta}_- \hat{\Delta}_+ \right| \varphi_s \right>$$
And $\left< \varphi_s \left| \hat{\Delta}_- \hat{\Delta}_+ \right| \varphi_s \right>$ converges for $Re(s) > \frac{1}{2}$.
This equation is valid if $s$ is a zero of Riemann zeta with real part greater than $\frac{1}{2}$. So, if we assume that a zero exists with real part greater than $\frac{1}{2}$, then the eigenvalue $-i\hbar (1 - 2s)$ has to be real, so the real part of $s$ has to be $\frac{1}{2}$. This contradictions means that there can not be zeros with real part greater than $\frac{1}{2}$. The Riemann hypothesis is true!
Edit: The last equation would still be valid for $s$ not being a zero of Riemann zeta. But then $\hat{G}^\dagger = \hat{\Delta}_-\hat{x}\hat{\Delta}_+ \hat{p}$ would not be the Hermitian of $\hat{G} = \hat{p}\hat{\Delta}_-\hat{x}\hat{\Delta}_+$ and then $\left< \varphi_s \left| \hat{G} + \hat{G}^\dagger \right| \varphi_s \right>$ would probably not be real. You can read in my paper why the hermiticy is there for the boundary condition $\varphi_s(1) = \zeta(s) = 0$ with $Re(s) > \frac{1}{2}$.
