What's wrong with this Hamiltonian matrix? Suppose we have a Hamiltonian matrix:
$$H=
\begin{pmatrix}
0&\tanh x-\partial_x\\
\tanh x+\partial_x&0\\
\end{pmatrix}
$$
Obviously, $H^\dagger_{ij}=H_{ji}$.
Two of the eigen-states of this matrix are:
$$\Psi^{(\pm)}=
\begin{pmatrix}
\tanh x\\
\pm 1\\
\end{pmatrix}
$$
from which one can easily find $H\Psi^{(\pm)}=\pm\Psi^{(\pm)}$
But one can also find 
$$
\sum_{i=1}^2\int
\Psi^{(+)\dagger}_i\Psi^{(-)}_i
\mathrm{d}x
\ne 0
$$ 
which indicates that the eigen-states are not othogonal to each other.
One may ascribe this perculiarity to the Hermicity of Hamiltonian, but it looks fine, even when we replace this $\tanh x$ by $1$, everything goes well. 
Can any one tell me where the problem is? thx.
 A: Let's recall why eigenvectors of Hermitian operators with different eigenvalues are orthogonal. If $A\vert\lambda_i\rangle=\lambda_i\vert \lambda_i\rangle$, then we have:
$$\langle \lambda_1\vert \lambda_2\rangle=\frac{1}{\lambda_1}\langle\lambda_1\vert A\vert\lambda_2\rangle=\frac{1}{\lambda_2}\langle\lambda_1\vert A\vert \lambda_2\rangle,$$
so if $\lambda_1\not =\lambda_2$ and both are nonzero then we have $\langle\lambda_1\vert\lambda_2\rangle=0$. Let's see what happens when we try this for our eigenvectors.
$$
\begin{split}
\langle\lambda_1\vert\lambda_2\rangle&=\begin{pmatrix}\tanh x&1\end{pmatrix}\begin{pmatrix}\tanh x\\-1\end{pmatrix}=\int -\text{sech}^2x\,dx=-2\\
\frac{1}{\lambda_2}\langle\lambda_1\vert A\vert\lambda_2\rangle&=\frac{1}{-1}\begin{pmatrix}\tanh x&1\end{pmatrix}\begin{pmatrix}0&\tanh x-\partial_x\\\tanh x+\partial_x &0\end{pmatrix}\begin{pmatrix}\tanh x\\-1\end{pmatrix}\\&=\int -\text{sech}^2x\,dx=-2\\
\frac{1}{\lambda_1}\langle\lambda_1\vert A\vert \lambda_2\rangle&=\left(\begin{pmatrix}0&\tanh x-\partial_x\\\tanh x+\partial_x &0\end{pmatrix}\begin{pmatrix}\tanh x\\1\end{pmatrix}\right)^\dagger\begin{pmatrix}\tanh x\\-1\end{pmatrix}\\
&=\begin{pmatrix}\tanh x\\1\end{pmatrix}^\dagger \begin{pmatrix}\tanh x\\-1\end{pmatrix}=-2
\end{split}
$$
At first glance perhaps everything looks OK, since they all are all equal. But something is wrong, since $\lambda_1\not =\lambda_2$, yet we have $\lambda_1 C=\lambda_2 C\not =0$!
We are led to the conclusion that $I_1\not =I_2$:
$$\begin{split}I_1&=\begin{pmatrix}\tanh x&1\end{pmatrix}\begin{pmatrix}0&\tanh x-\partial_x\\\tanh x+\partial_x &0\end{pmatrix}\begin{pmatrix}\tanh x\\-1\end{pmatrix}\\I_2&=\left(\begin{pmatrix}0&\tanh x-\partial_x\\\tanh x+\partial_x &0\end{pmatrix}\begin{pmatrix}\tanh x\\1\end{pmatrix}\right)^\dagger\begin{pmatrix}\tanh x\\-1\end{pmatrix}\end{split}$$
The operator is not behaving as we would expect a Hermitian operator to behave. As Qmechanic points out, this is because $\tanh x$ and 1 do not go to zero as $x\to \infty$, so there are surface terms when you do the integration by parts to convert from $I_1$ to $I_2$. 
(added in an edit)
Let's be more precise about this:
$$\begin{split}
I_1&=\int dx\, (\tanh x)\left((\tanh x)-\partial_x(-1)\right)+\left((\tanh x)(\tanh x)+\partial_x(\tanh x)\right)\\
&=\int dx\,(\tanh x)^2-(\tanh x)\partial_x(-1)+(\tanh x)^2+(\partial_x (\tanh x))1\\
&=\left(\int dx\,(\tanh x)^2\right)-\left(\left[(\tanh x)(-1)\right]-\int dx\,(-1)\partial_x\tanh x\right)+\int dx\,(\tanh x)^2+\\
&+[\tanh x]-\int dx\,(\tanh x)\partial_x 1\\
&=\int dx\,(\tanh x)^2+(-1)\partial_x \tanh x+(\tanh x)^2-(\tanh x)\partial_x 1+2[\tanh x]\\
&=I_2+2[\tanh x]
\end{split}$$
However, the more I think about this answer, the more I am convinced it is incomplete. In particular, we may put the system in a box and impose periodic boundary conditions, so that all these issues about integration by parts go away. I'm not mathematically adept enough to solve the boxed problem analytically, but looking at the numerics suggests that the spectrum is very much comparable to the spectrum udrv obtained for the unboxed problem; i.e., "plane waves where one of them is multiplied by $\tanh$," except in the boxed problem they're all orthogonal. I don't have a good explanation for why this no longer holds in the infinite limit.

A: Hamiltonian H is indeed self-adjoint and its eigenvectors must be orthogonal. The reason why a subset of the "eigenbasis" appears to be non-orthogonal has to do with the particular structure assumed for the eigenvectors and with their degeneracy. The very short answer is that their two components carry the wrong normalization. But it is more interesting to see why this is so on algebraic grounds. 
For this, let us try to construct the eigenbasis in a somewhat more systematic way:
1. Preamble (reproducing from comments): 
First rewrite the Hamiltonian in the compact form 
$$
{\hat H} = \tanh(x) {\hat \sigma}_1 + {\hat \sigma}_2 {\hat p}_x 
$$
$$
{\hat p}_x = -i \partial_x,\;\;\;{\hat \sigma}_1 = \left(\begin{array}{cc}0 & 1 \\1 & 0\end{array}\right),\;\;\; {\hat \sigma}_2 = \left(\begin{array}{cc}0 & -i \\i & 0\end{array}\right), \;\;\;{\hat \sigma}_3 = \left(\begin{array}{cc}1 & 0 \\0& -1\end{array}\right)
$$
and cast the eigenvectors as 
$$
{\hat \Psi} = \left(\begin{array}{c}\psi \\ \phi\end{array}\right)
$$
This allows us to note that $\lbrace {\hat H}, {\hat \sigma}_z\rbrace = 0$, with the immediate corollary that if ${\hat \Psi}$ is an eigenvector for eigenvalue $\lambda$, then ${\hat \sigma}_3{\hat \Psi}$ is an eigenvector for $-\lambda$, 
$$
H({\hat \sigma}_3 {\hat \Psi}) = (H{\hat \sigma}_3){\hat \Psi} = -{\hat \sigma}_3 {\hat H}{\hat \Psi} = - \lambda({\hat \sigma}_3{\hat \Psi})
$$ 
So if the spectrum of $H$ is unbounded from above, it is also unbounded from below. 
2. Structure of eigenbasis: 
Consider now a positive eigenvalue $\lambda > 0$ and the related eigenvalue equations
$$
{\hat H}{\hat \Psi}_\lambda = \lambda{\hat \Psi}_\lambda
$$
$$
{\hat H}({\hat \sigma}_3{\hat \Psi}_\lambda) = -  \lambda({\hat \sigma}_3{\hat \Psi}_\lambda)
$$
If we rearrange and separate the orthogonal "up" and "down" components $\frac{{\hat I}+{\hat \sigma}_3}{2}{\hat \Psi}_\lambda = \left(\begin{array}{c}\psi_\lambda \\ 0\end{array}\right)$, $\frac{{\hat I}-{\hat \sigma}_3}{2}{\hat \Psi}_\lambda = \left(\begin{array}{c}0 \\ \phi_\lambda\end{array}\right)$,
$$
{\hat H}\left[\frac{{\hat I}+{\hat \sigma}_3}{2}{\hat \Psi}_\lambda\right] = \lambda\left[\frac{{\hat I}-{\hat \sigma}_3}{2}{\hat \Psi}_\lambda\right] 
$$ 
$$
{\hat H}\left[\frac{{\hat I}-{\hat \sigma}_3}{2}{\hat \Psi}_\lambda\right] = \lambda\left[\frac{{\hat I}+{\hat \sigma}_3}{2}{\hat \Psi}_\lambda\right]
$$
we find that these are in fact degenerate eigenvectors of ${\hat H}^2$ for positive eigenvalue $\lambda^2$:
$$
{\hat H}^2\left[\frac{{\hat I}+{\hat \sigma}_3}{2}{\hat \Psi}_\lambda\right] = \lambda{\hat H}\left[\frac{{\hat I}-{\hat \sigma}_3}{2}{\hat \Psi}_\lambda\right] =  \lambda^2 \left[\frac{{\hat I}+{\hat \sigma}_3}{2}{\hat \Psi}_\lambda\right]
$$
$$
{\hat H}^2\left[\frac{{\hat I}-{\hat \sigma}_3}{2}{\hat \Psi}_\lambda\right] = \lambda{\hat H}\left[\frac{{\hat I}+{\hat \sigma}_3}{2}{\hat \Psi}_\lambda\right] =  \lambda^2 \left[\frac{{\hat I}-{\hat \sigma}_3}{2}{\hat \Psi}_\lambda\right]
$$
The reason why we can separate the components is the form of ${\hat H}^2$ itself,
$$
{\hat H}^2 = \left[ \tanh(x) {\hat \sigma}_1 + {\hat \sigma}_2 {\hat p}_x \right] \left[ \tanh(x) {\hat \sigma}_1 + {\hat \sigma}_2 {\hat p}_x \right] =
$$ 
$$
= {\hat I}\left[ \tanh^2(x) - \partial^2_x\right] - i{\hat \sigma}_1{\hat \sigma}_2 \tanh(x)  \partial_x - i {\hat \sigma}_2 {\hat \sigma}_1 (\partial_x \tanh(x) ) = 
$$
$$
= {\hat I}\left[ \tanh^2(x) - \partial^2_x\right] + {\hat \sigma}_3 \left[ \tanh(x)  \partial_x - (\partial_x \tanh(x) ) \right]= 
$$
$$
= {\hat I}\left[ \tanh^2(x) - \partial^2_x\right] - {\hat \sigma}_3 \frac{1}{\cosh^2(x)} =
$$
$$
= \left[\frac{{\hat I}+{\hat \sigma}_3}{2}\right]\left[ 2\tanh^2(x) -1 - \partial^2_x\right] + \left[\frac{{\hat I}-{\hat \sigma}_3}{2}\right]\left[1 - \partial^2_x\right]
$$
This tells us a few remarkable things at a glance:
1) Components $\psi_\lambda$, $\phi_\lambda$ now appear to satisfy decoupled eigenvalue equations
$$
\left[ 2\tanh^2(x) -1 - \partial^2_x\right] \psi_\lambda =  \lambda^2 \psi_\lambda
$$ 
$$
\left[1 - \partial^2_x\right] \phi_\lambda = \lambda^2 \phi_\lambda
$$
2) The 2nd equation above is a simple Helmholtz equation with plane-wave bounded orthogonal solutions 
$$
\phi_{\lambda, \pm} \sim e^{\pm i k_\lambda x}, \;\;\;k_\lambda >0, 
$$
$$
\langle \phi_{\lambda, \pm}| \phi_{\lambda', \mp}\rangle = 0, \;\;\; \langle \phi_{\lambda, \pm}| \phi_{\lambda', \pm}\rangle = \delta(k_\lambda - k_\lambda')
$$
and eigenvalues 
$$
\lambda^2 = 1 + k^2_\lambda >1
$$
3) With this form for $\lambda^2$ the equation for 1st component becomes 
$$
-  \left[ \frac{2}{\cosh^2(x)} + \partial^2_x\right] \psi_\lambda =  k^2_\lambda \psi_\lambda
$$ 
and can be shown to have bounded orthogonal solutions of the form 
$$
\psi_{\lambda, \pm} \sim \left[ \tanh(x) \mp i k_\lambda \right] e^{\pm i k_\lambda}
$$
$$
\langle \psi_{\lambda, \pm}| \psi_{\lambda', \mp}\rangle = 0, \;\;\; \langle \psi_{\lambda, \pm}| \psi_{\lambda', \pm}\rangle = \delta(k_\lambda - k_\lambda')
$$
Note the individual normalization, in distributions sense, of $\psi_{\lambda, \pm}$ and $\phi_{\lambda, \pm}$. We'll see in a moment that this is actually the fix to the non-orthogonality of the eigenvectors proposed in the OP.
4) The four corresponding eigenvectors of $H^2$, for given $\lambda^2$,
$$
\left(\begin{array}{c} \psi_{\lambda, \pm} \\0\end{array}\right),\;\;\;\left(\begin{array}{c}0 \\  \phi_{\lambda, \pm} \end{array}\right)
$$
provide an orthonormal basis (in the distributions sense), albeit not yet an eigenbasis, for the $\pm \lambda$-eigensubspace of $\hat H$:
$$
{\hat H} \left(\begin{array}{c} \psi_{\lambda, \pm} \\0\end{array}\right) = \lambda \left(\begin{array}{c}0 \\  \phi_{\lambda, \pm} \end{array}\right)
$$
$$
{\hat H} \left(\begin{array}{c}0 \\  \phi_{\lambda, \pm} \end{array}\right) = \lambda \left(\begin{array}{c} \psi_{\lambda, \pm} \\0\end{array}\right)
$$
Since the + and - eigenvectors remain mutually decoupled, the subtended 4-dimensional subspace decomposes further into two disjoint 2-dimensional subspaces. 
3. Answer to the original Question: Diagonalization of ${\hat H}$: In order to finally diagonalize ${\hat H}$, now we only have to consider its matrix representation in the two well-defined bases {$\psi_{\lambda, \pm}$, $\phi_{\lambda, \pm}$}. In both cases it amounts simply to
$$
\left(\begin{array}{cc}0 & \lambda \\ \lambda & 0\end{array}\right) \equiv \lambda {\hat \sigma}_1
$$
and gives the corresponding eigenvalues as $\pm \lambda$, as expected, and the orthonormal eigenvectors as
$$
{\hat \Psi}_{\lambda, \pm} =  \frac{1}{\sqrt2}\left(\begin{array}{c} \psi_{\lambda, \pm} \\  \phi_{\lambda, \pm} \end{array}\right), \;\;\;\text{for}\;\;(+\lambda)
$$
and 
$$
{\hat \sigma}_3 {\hat \Psi}_{\lambda, \pm} = \frac{1}{\sqrt2}\left(\begin{array}{c} \psi_{\lambda, \pm} \\  - \phi_{\lambda, \pm} \end{array}\right) \;\;\;\text{for}\;\;(-\lambda)
$$
At first sight this looks just as in the comments and the OP, with the "tiny" difference that here the $\psi_{\lambda, \pm}$-s and $\phi_{\lambda, \pm}$-s are already individually normalized (and mutually orthogonal). But this turns out to be the crucial point, since unlike before, now we have not only 
$$
\langle {\hat \Psi}_{\lambda, \pm} |{\hat \Psi}_{\lambda', \pm}\rangle = \langle {\hat \sigma}_3 {\hat \Psi}_{\lambda, \pm} | {\hat \sigma}_3 {\hat \Psi}_{\lambda', \pm}\rangle \sim \delta(\lambda - \lambda'), \;\;\;\langle {\hat \Psi}_{\lambda, \pm} |{\hat \sigma}_3|{\hat \Psi}_{\lambda', \mp}\rangle = 0
$$
but also
$$
\langle {\hat \Psi}_{\lambda, \pm} |{\hat \sigma}_3|{\hat \Psi}_{\lambda', \pm}\rangle = \frac{1}{2} \left(\; \langle \psi_{\lambda, \pm} |\;\;\;\langle \phi_{\lambda, \pm}| \;\right) {\hat \sigma}_3 \left(\begin{array}{c} |\psi_{\lambda', \pm}\rangle \\  |\phi_{\lambda', \pm} \rangle \end{array}\right) = 
$$
$$
= \frac{1}{2} \left(\; \langle \psi_{\lambda, \pm} |\;\;\;\langle \phi_{\lambda, \pm}| \;\right)  \left(\begin{array}{c} |\psi_{\lambda', \pm}\rangle \\ - |\phi_{\lambda', \pm} \rangle \end{array}\right) = \frac{1}{2} \left( \; \langle \psi_{\lambda, \pm} |\psi_{\lambda', \pm}\rangle - \langle \phi_{\lambda, \pm}|\phi_{\lambda', \pm} \rangle\; \right) = 0, \;\;\; \forall\;\lambda, \lambda'
$$
The positive and negative energy solutions are at long last mutually orthogonal in a natural way. Problem solved.
Note on the finite-box version: Paradoxically, constraining the problem to a finite box raises much more interesting problems concerning the correct boundary conditions required for Hermiticity and self-adjointness (please note in passing that the boundary conditions at infinity did not enter the above argument in any non-trivial way). This is because it turns out that imposing standard Dirichlet or periodic bc-s on one of the eigenvector components, for instance the one satisfying the Helmholtz equation, does not result in similar bc-s (Dirichlet or periodic) for the other component. This is already visible in the graphs posted by @commutatortot for periodic bc-s: True periodic conditions must match both endpoint values and endpoint derivatives, for both spinor components. But in this problem it can be verified that if one component satisfies periodic bc-s for both endpoint values and endpoint derivatives, then the other component can only satisfy periodic bc-s either for endpoint values or for endpoint derivatives, but not for both. In commutatortot's graphs, the 1st component is truly periodic, while the 2nd one shows matching endpoint values, but not matching derivatives.The issue is beyond the scope of this question, but I think it is worth mentioning: if I didn't miss something important, it seems that Hermiticity requires bc-s that do not seem to have a standard physical justification.      
