I'm following Evans here. The statement is as follows: take a subsolution $u\in C^2(U)\cap C^1(\overline{U})$ of a 2nd order elliptic operator $L$ in non-divergence form. Suppose also $U$ satisfies the interior ball condition. If, for example, the constant term of $L$ is identically $0$, then the normal derivative at a maximum point $x_o\in \partial U$ is strictly positive.
The proof, roughly speaking, is to perform a perturbation of $u$ by a subsolution $v$ (essentially the Gaussian, but modified so that it vanishes on $\partial B$) such that $v|_{\partial B}\equiv 0$ and $\frac{\partial v}{\partial N}|_{x_o}<0$ for $N$ the unit outward normal. Here $B=B_x(R)$, for $x\in U$ such that $x_o\in \partial B$. Applying weak maximum principle essentially gets you the conclusion. A technical point here, however, is that actually we cannot take $v$ to be defined on all of $B$, but rather the outer annulus $A:=B_x(R)-B_x(\frac{R}{2})$. The reason this is the case is that Evans can only guarantee that $v$ remains a subsolution on $A$ as he can make it behave like a second order polynomial with negative leading coefficient, but only away from a neighborhood of the origin. Therefore, it is conceivable to me that one may choose another subsolution $v$ on all of $B$, instead of an annular neighborhood away from the origin. Is this possible, or is there some fundamental obstruction in finding such a $v$? Another question I have that comes from this proof is that it's surprising to me that (essentially) the Gaussian is a subsolution to any 2nd order elliptic operator $L$. Is there any reason I should've expected this to be the case a priori? I feel like there may be a connection between the second question and the Fourier transform, and I would appreciate any insight.
First question: Yes, there is a fundamental obstruction: the weak maximum principle! As part of the proof you require that $v=0$ on $\partial B_r(x)$ (or possibly $v\leqslant 0$ on $\partial B_r(x)$) so that you can conclude that $u \geqslant v$ on $\partial B_r(x)$ in order to satisfy the assumptions of the weak maximum principle; however, if $v$ is a subsolution in all of $B_r(x)$ then $v\leqslant 0$ by the weak maximum principle which implies that $\frac{\partial v}{\partial \nu } \geqslant 0$ on $\partial B_r(x)$.
Second question: To be honest, I don’t think there is anything particularly special about the Gaussian per se, but rather the convexity properties of the Gaussian.
For example, let’s consider the simplest case $n=1$ and $B_r(x)=(-1,1)$. We want some $v$ such that $-v’’ \leqslant 0$ in $(-1,-1/2) \cup (1/2,1)$, $v(-1)=v(1)=0$, and $v>0$ in $(-1,-1/2)$ (which, using that $v(-1)=0$ gives that $v’(-1)=\frac{\partial v}{\partial \nu}(-1)<0$). Let $\varphi \in C^2(\mathbb R)$ be an even function to be chosen later, $\lambda >0$, and $v_\lambda (x) := \varphi(\lambda x ) - \varphi(\lambda)$. If $Lv=-a(x)v’’+b(x)v’+c(x)$ with $a>0$, we have that \begin{align*} Lv_\lambda (x) &= -a\lambda^2 \varphi’’(\lambda x ) + b \lambda \varphi’(x)+c(\varphi(\lambda x)-\varphi(\lambda)). \end{align*} So long as $\varphi’’(\lambda x)\geqslant k >0$ for $x \in (-1,-1/2)$, $\lambda $ large and $k$ independent of $\lambda$ (which equivalent to being strictly convex in this region), and $\varphi,\varphi’$ are bounded then you can conclude that $Lv_\lambda \leqslant 0$ in $(-1,-1/2)\cup (1/2,1)$. Then you would just need to find such a $\varphi$ that is nonnegative which there would be plenty. Another example of such a function is $\varphi(x)=\frac 1 {1+x^2}$ which looks kind of like a Gaussian but isn’t. In $\mathbb R^n$, you could make the exact same argument except with radial functions.