In nonrelativistic Quantum Mechanics, is the expectation value of a sum of operators always equal to the sum of the expectation values?

Question

Suppose that $\lvert \psi_n \rangle$ are the eigenvectors of a Hamiltonian, $\hat{H}$, which span some Hilbert space $\mathcal{H}$ and satisfy $$\hat{H}\lvert \psi_n \rangle = E_n \lvert \psi_n \rangle.$$

Since all Hilbert spaces are vector spaces and thus linear, I would expect that for any operator $$\hat{O} = \hat{A} + \hat{B},$$ we should have that $$\langle \psi_n \rvert \hat{O} \lvert \psi_n \rangle = \langle \psi_n \rvert \hat{A}+\hat{B} \lvert \psi_n \rangle = \langle \psi_n \rvert \hat{A} \lvert \psi_n \rangle + \langle \psi_n \rvert \hat{B} \lvert \psi_n \rangle.$$

Are there any pathological cases where this simple identity fails? For example, when the Hamiltonian contains some distribution like $\delta(x)$?

Proposed counter-example: this question is motivated by my calculations following this paper on the systems of two identical atoms with s-wave contact interaction in a harmonic trap. The authors provide an analytical solution for this problem which I will employ. All following calculations will be expressed in natural harmonic oscillator units.

After decoupling the center of mass and relative motion, the relative Hamiltonian is given by $$H_\text{rel}= -\frac{1}{2} \nabla_r^2 +\frac{1}{2}r^2 + \sqrt{2} \pi a \delta^3(\vec{r}) \frac{\partial}{\partial r} r = H_\text{osc} + V_\text{pseudo}.$$ The solution (for example) for the ground state with energy $E_0 = -\frac{1}{2}$, corresponding to $\nu = -1$, is $$\psi_0(r) = \frac{e^{\frac{r^2}{2}}}{2 \pi^{3/4} \sqrt{\ln{2}}}\big(\text{ExpIntegralE}(\frac{1}{2},r^2)\big), $$ with $\lvert \vec{r} \rvert = r =\lvert \frac{1}{\sqrt{2}}(\vec{r}_1 - \vec{r}_2) \rvert$ and $\text{ExpIntegralE}(x)$ is the exponential integral function in Mathematica. The scattering length, $a$ associated with this solution is $a = \sqrt{\frac{\pi}{2}}$

It is easy to verify that in fact $$H_{rel} \psi_0(r) = -\frac{1}{2} \psi_0(r),$$ provided that one correctly accounts that the term proportional to $\delta^3(\vec{r})$ is cancelled by a corresponding term coming from the Laplacian. Since $$\psi_0(r) \sim_{r \to 0} \frac{1}{r},$$ we have to treat the Laplacian in a distributional sense such that $$\nabla^2_r \psi_0(r) = \nabla^2_{r \neq 0}\psi_o(r) -4 \pi \delta^3(\vec{r}) \text{Res}(\psi_0(0)).$$ $\text{Res}(f(r))$ is the residue of $f$ at $r$.

We should straight-forwardly have that $$\langle \psi_0 \rvert \hat{H}_{rel} \lvert \psi_0 \rangle = \langle \psi_0 \rvert E_0 \lvert \psi_0 \rangle = E_0.$$ However, when I break apart the Hamiltonian operator, I find that I get additional non-zero terms coming from the terms proportional to $\delta^3(\vec{r})$.

In particular, I find that $$ 4 \pi \int_0^{\infty} \psi_0(r) (-\frac{1}{2} \nabla_{r \neq 0}^2 +\frac{1}{2}r^2) \psi_0(r) r^2 \ dr = E_0,$$ whereas the remaining terms give $$ \int \psi_0(r) \Big(-2 \pi \delta^3(\vec{r})\text{Res}(\psi_0(0)) + \sqrt{2} \pi \sqrt{\frac{\pi}{2}} \delta^3(\vec{r}) \frac{\partial}{\partial r}\big(r \psi_0(r)\big) \Big) \ d^3 r = \frac{\pi}{4 \ln{2}}. $$

I cannot justify ignoring these terms since they were crucial in proving that $\psi_0$ is a solution to the Schrodinger equation with the correct eigenvalue, nor can I find any unaccounted for counter-terms.

Answer 1

Your last integral is $$ J = \int{d^3r \;\psi_0(r) \left[-2\pi\delta(\vec{r})Res(\psi_0(0)) + \sqrt{2}\pi a\delta(\vec{r})\frac{\partial}{\partial r}(r\psi_0(r))\right]} =\\ = 4\pi \int_0^\infty{dr\; r^2 \psi_0(r) \left[-2\pi\frac{\delta(r)}{2\pi r^2}Res(\psi_0(0)) + \sqrt{2}\pi a\frac{\delta(r)}{2\pi r^2}\frac{\partial}{\partial r}(r\psi_0(r))\right]} \\ = 4\pi\int_0^\infty{dr\;\psi_0(r)\left[-\delta(r)Res(\psi_0(0)) + \frac{ a}{\sqrt{2}}\delta(r)\frac{\partial}{\partial r}(r\psi_0(r))\right]} = \\ = 4\pi \int_0^\infty{dr\;\delta(r)\psi_0(r)\left[-Res(\psi_0(0)) + \frac{ a}{\sqrt{2}}\frac{\partial}{\partial r}(r\psi_0(r))\right]} $$ where I used $\delta(\vec{r}) = \frac{\delta(r)}{2\pi r^2}$ (see here). For $\psi_0(r) \approx \frac{\alpha}{r}$ when $r\rightarrow 0$ the last form above becomes $$ J = 4\pi \int_0^\infty{dr\;\delta(r)\left(\frac{\alpha}{r}\right)\left[-Res(\psi_0(0)) + \frac{ a}{\sqrt{2}}\frac{\partial}{\partial r}\alpha\right]} = \\ = - 4\pi \alpha \;Res(\psi_0(0)) \int_0^\infty{dr\;\frac{\delta(r)}{r}} = \\ = 4\pi \alpha \;Res(\psi_0(0)) \int_0^\infty{dr\;\delta'(r)} = 4\pi \alpha \;Res(\psi_0(0)) \frac{d}{dr}(1) = 0 $$ The first expression on the last line uses again a delta-identity, this time $\delta(r) = -r \delta'(r)$. The result obviously holds in distributional sense. This leaves $$ \langle \psi_0|\hat{H}_0|\psi_0\rangle = E_0 $$ as it should. I might have missed a factor of $\frac{1}{2}$ in the last integral over $\delta'$, but it wouldn't matter.

NOTE IN REPLY TO COMMENT:

The above takes into account only the dominant contribution to $\psi_0$ for $r\rightarrow 0$ before proceeding with the integral. It has been objected that this is not a valid procedure since taking into account zero and higher order terms in $\psi_0$ would render the integral nonzero. For instance, using the identity $\delta(r) = -r \delta'(r)$ in the 4th line from the top above yields $$ -4\pi \int_0^\infty{dr\;\delta'(r)\left(r\psi_0(r)\right)\left[-Res(\psi_0(0)) + \frac{ a}{\sqrt{2}}\frac{\partial}{\partial r}(r\psi_0(r))\right]} = \\ = 4\pi \lim_{r\rightarrow 0}\frac{\partial}{\partial r}\left[\left(r\psi_0(r)\right)\left(-Res(\psi_0(0)) + \frac{ a}{\sqrt{2}}\frac{\partial}{\partial r}(r\psi_0(r))\right)\right] \neq 0 $$ since both $r\psi_0$ and $\frac{\partial}{\partial r}r\psi_0$ retain terms linear in $r$ as $r\rightarrow 0$. If this is indeed correct, then $\langle \psi_0|\hat{H}_0|\psi_0\rangle = E_0$ is compromised.

I don't think this is the case, for the following reason. The original equation $\hat{H}_0\psi_0 = E_0\psi_0$ is satisfied provided the laplacian is taken in the distributional sense and the singular pseudo-potential term is taken to mean $$\sqrt{2}\pi a\delta(\vec{r})\frac{\partial (r\psi_0)}{\partial r}(0)$$ This is just the standard interpretation of the singular term. If it is taken into account from the beginning when writing the integral for $\langle\psi_0|\hat{H}_0|\psi_0\rangle$, the integral $J$ discussed above would not appear in the first place and everything would check out trivially. So the terms that make $J$ non-trivial here are the same ones that are wiped out under the singularity in the original differential equation. In view of this, we can say that dismissing higher order contributions in $\psi_0$ when integrating the singularities functions in the same way: it wipes out irrelevant terms accompanying the singularity. On the other hand, including these terms in the integral is no longer as benign as carrying them along in the context of the differential equation. This is simply because they change the nature of the function on which the singular distribution acts. So basically the problem comes down to $$ \int{d^3r \;\psi_0(r) \left[-2\pi\delta(\vec{r})Res(\psi_0(0)) + \sqrt{2}\pi a\delta(\vec{r})\frac{\partial}{\partial r}(r\psi_0(r))\right]} \\ \neq \\ \int{d^3r \;\psi_0(r) \left[-2\pi\delta(\vec{r})Res(\psi_0(0)) + \sqrt{2}\pi a\delta(\vec{r})\frac{\partial(r\psi_0)}{\partial r}(0)\right]} $$