What is the physical/deeper reason for the second order shift of the ground state energy in time independent perturbation theory to be always down?
I know that it follows from the formula quite straightforwardly, but I could not find a deeper reason for it even if if I took up examples like the linear Stark effect.
Would you please tell me why?
I) The lowering of the ground state energy is a special case of the more general phenomenon of level repulsion (because the excited energy levels by definition must be larger than the ground state energy).
II) Level repulsion is not just a quantum phenomenon. It also happens for purely classical systems, e.g. two coupled oscillators, as mentioned in the link.
III) Mathematically, the non-zero off-diagonal elements ("the interactions") of an $n\times n$ Hermitian matrix $(H_{ij})_{1\leq i,j \leq n}$ ("the Hamiltonian") cause the real eigenvalues $E_1$, $\ldots$, $E_n$, (counted with multiplicity) to, loosely speaking, spread out more than the distribution of diagonal elements $H_{11}$, $\ldots$, $H_{nn}$. This eigenvalue repulsion effect is e.g. encoded in the Schur-Horn Theorem.
IV) Perhaps a more physical understanding of level repulsion is as follows. When considering a perturbation problem $H=H^{(0)}+V$, we would like to find the true energy eigenstates $|1 \rangle$, $\ldots$, $|n \rangle$, with energy eigenvalues $E_1$, $\ldots$, $E_n$. A priori we only know the unperturbed energy eigenstates $|1^{(0)} \rangle$, $\ldots$, $|n^{(0)} \rangle$, with energy eigenvalues $H_{11}$, $\ldots$, $H_{nn}$. (We have for the sake of simplicity assumed that the interaction part $V$ has no diagonal part. This is always possible via reorganizing $H^{(0)} \leftrightarrow V$.) Thus an unperturbed energy eigenstates
$$|i^{(0)} \rangle ~=~ \sum_{j=1}^n |j \rangle~ \langle j |i^{(0)} \rangle$$
is really a linear combination of the true energy eigenstates $|1 \rangle$, $\ldots$, $|n \rangle$. The square of the overlaps $\langle j |i^{(0)} \rangle$ has a probabilistic interpretation
$$ \sum_{j=1}^n |\langle j |i^{(0)}\rangle|^2 ~=~\langle i^{(0)}|i^{(0)}\rangle ~=~1. $$
Hence the unperturbed energy eigenvalue
$$H_{ii}~=~\langle i^{(0)}|H|i^{(0)}\rangle ~=~\sum_{j=1}^n E_j|\langle j |i^{(0)}\rangle|^2 $$
is a quantum average of the true energy eigenvalues $E_1$, $\ldots$, $E_n$ of the system. Intuitively, the distribution of unperturbed energy eigenvalues $H_{11}$, $\ldots$, $H_{nn}$, must therefore, loosely speaking, be closer to each other than the distribution of true energy eigenvalues $E_1$, $\ldots$, $E_n$.
V) Finally, let us mention that eigenvalue repulsion plays an important role in random matrix theory. Integrating out "angular" d.o.f. leads to a Vandermonde measure factor
$$ \prod_{1\leq i<j \leq n} |E_j-E_i|^{\beta}, $$
so that the partition function ${\cal Z}$ favors that the eigenvalues $E_1$, $\ldots$, $E_n$ are different. Here the power $\beta$ is traditionally $1$, $2$, or $4$, depending on the random matrix ensemble. For Hermitian matrices $\beta=2$. See also this post.
The true reason is the "perturbation" operator $\hat{V}$ being accidentally Hermitian (self-conjugate) in the linear space of non-perturbed Hamiltonian $\hat{H}_0$ in many practical cases: $V_{mn}=V_{nm}^*.$
Generally, a perturbation theory develops a function $f$ in a Taylor series in powers of a small parameter, let us say, $\varepsilon$, and this function is not obliged to have its second derivative of a negative sign at the expansion point $\varepsilon = 0$. The function $f(\varepsilon)$ may be concave or convex at the expansion point.
In terms of operators, the total Hamiltonian $\hat{H}=\hat{H}_0+\hat{V}$ should be Hermitian in the linear space of its exacts eigenfunctions $\varphi_n$, not in the linear space of the non perturbed Hamiltonian $\hat{H}_0$. So, the perturbation operator $\hat{V}$ may be non Hermitian in the space of non perturbed eigenfunctions $\varphi_n^{(0)}$. Then, the second-order correction may take a positive sign.
I encountered such an example in my practice. See it here, Chapter 2.1.
