Let $f\colon\mathbb{R}^n\to\mathbb{R}^n$ be a smooth, bounded vector field. Further, let $u\colon\mathbb{R}^n\to\mathbb{R}$ satisfy $$-\Delta u=\operatorname{div}(fu).$$ If $u\in L^1(\mathbb{R}^n)$, then $u$ has one sign, i.e., either $u>0$ everywhere, or $u<0$ everywhere, or $u=0$ everywhere.

I have a direct proof of this for $n=1$. For $n>1$, I have a proof using the theory of parabolic equations (see below). My question: Is there a direct proof using only the theory of elliptic PDEs?

(Edited to assume $f$ bounded and fix the case $n=1$ below.)

For $n=1$, my proof goes as follows. The equation is $-u''=(fu)'$, which integrates to $u'+fu=A$ for some constant $A$. If $u$ changes sign then we may without loss of generality take $u(0)=0$. Thus $$u(x)=A\int_0^x e^{F(t)-F(x)}\,dt$$ where $F'=f$. If $|f|\le c$ then $F(t)-F(x)\ge c(t-x)$ for $t<x$, so $u(x)\ge Ac^{-1}(1-e^{-cx})$ when $x>0$, and hence $u\notin L^1$ (unless $A=0$).

For $n>1$, my only proof is much more involved. Here is a brief outline. Assume the conclusion is wrong, so we can write $u=u_+-u_-$ with $u_\pm\ge0$ everywhere and neither identically zero, and $u_+u_-=0$ everywhere.

Now let $v_\pm$ solve $$\frac{\partial v_\pm}{\partial t}=\Delta v_\pm+\operatorname{div}(fv_\pm)$$ for $t>0$, with initial conditions $v_\pm(0,x)=u_\pm(x)$. By uniqueness for this equation (with suitable growth conditions at infinity), $v_+(t,x)-v_-(t,x)=u(x)$ for $t>0$ and $x\in\mathbb{R}^n$. Also, for $t>0$ we find $v_\pm>0$ everywhere, and also $$\int_{\mathbb{R}^n} v_\pm(t,x)\,dx=\int_{\mathbb{R}^n} u_\pm(x)\,dx$$ since the equation is on divergence form. We conclude $$\int_{\mathbb{R}^n}|u(x)|\,dx=\int_{\mathbb{R}^n}(u_+(x)+u_-(x))\,dx=\int_{\mathbb{R}^n}(v_+(t,x)+v_-(t,x))\,dx>\int_{\mathbb{R}^n}|v_+(t,x)-v_-(t,x)|\,dx,$$ which is a contradiction.