Does the orthogonality of states apply when there is a potential? I know this from a quantum mechanics class:
$$\langle \psi_i | \psi_j\rangle = \delta_{ij}$$
But does it also apply for:
$$\langle \psi_i |V| \psi_j\rangle = \delta_{ij}$$
where $V$ is any arbitrary potential? That is, if I try to compute the expectation value of $V$ using 2 different states $\psi_i$ and $\psi_j$ I get $0$?

The context:
I am trying to compute the eigenvalues of the particle in a box with an added potential $V$. To do this, the strategy I decided on was to compute like so:
$$\langle H_{total} \rangle = \langle -\frac{\hbar ^2}{2m}\frac{\partial ^2}{\partial x^2} + V \rangle$$
$$\langle H_{total} \rangle = \langle -\frac{\hbar ^2}{2m}\frac{\partial ^2}{\partial x^2} + (V_0 + V_{new}) \rangle \quad ; \quad V_0 \text{  is the potential for the particle in a box, which is zero, and  }V_{new} \text{   is the added potential.} $$
$$\langle H_{total} \rangle = \langle H_0 \rangle + \langle V_{new} \rangle$$
These are all matrices. So I'd fill the diagonal of the $H_0$ matrix using the analytic equation for the total energy and then fill the $V_{new}$ matrix, then add the two matrices together. Finally, I would compute the eigenvalues of this final matrix.
So where does my question come into all of this? Well, I know that the only values in $H_0$ are on the diagonal, but I didn't know if that also applies to the $V_{new}$ matrix as well.
 A: No. If you have the regular particle in a one-dimensional box of length $L$, the free Hamiltonian has eigenvalues $|\psi_n\rangle$ with $\langle x | \psi_n \rangle=\sqrt\frac{2}{L} \sin\left(\frac{n \pi x}{L}\right)$ . If, say, $V=V_0 \frac{x}{L}$, then you can do the integral and find $\langle \psi_1|V|\psi_2\rangle=-\frac{16}{9\pi} V_0$.
Some of your equations aren't correct, by the way. The hamiltonian you're interested in is this:
$$ H_{total} = -\frac{\hbar ^2}{2m}\frac{\partial ^2}{\partial x^2} + V$$
In which case the expectation value $\langle H\rangle=\langle \psi|H|\psi\rangle$ is actually:
$$ \langle H_{total}\rangle = \int_{0}^L\left(-\psi^*(x)\frac{\hbar ^2}{2m}\psi''(x) + V(x)\psi^*(x) \psi(x)\right)dx$$
where $\psi(x)$ is some wavefunction.
Furthermore, $H_{total}$ itself isn't a matrix, it's a linear operator which acts on an infinite dimensional space! If you discretize space into $N$ points, your Hamiltonian can be approximated numerically by an $N\times N$ matrix, but as you have encountered, this matrix is not trivial to diagonalize!
If you ask another question I can show you how to do this numerically and give examples using Mathematica and/or javascript. (I choose these two languages specifically because they are easy to visualize, Mathematica has the Eigensystem command (which can diagonalize H no problem), and there is a javascript library called numericjs which also has an eigensystem command).
A: No, it can not be generalized. If your added potential is harmonic then the answer will not be zero. You could write x in terms of annihilation and creation operators and then compute the expectation value of V. You will see it will not be zero generally. In this case, it will further depend on the states upon which you are performing the operations.
