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_{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_{osc} + V_{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.

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.