Existence of unique critical points to second order elliptic PDEs

Question

Let $\Omega\subset\mathbb{R}^2$ is bounded and convex and $\partial\Omega$ be smooth. Consider the second order elliptic PDE (1) $$ \begin{cases} Lu = f &\text{ on } \Omega\,,\\ u=0 &\text{ on } \partial\Omega\,, \end{cases} $$ where $f\in C^\infty(\Omega)$ satisfies $f(x) > 0$ and $f$ has no isolated critical points on $\Omega$. Then, does the solution $u$ of the above have only one critical point? I am looking for a result similar to Theorem 2 of this paper(ScienceDirect):

Theorem 2: Let $\Omega$ be a bounded, strictly convex domain in $\mathbb{R}^2$ and $u\in C^3(\Omega)\cap C^1(\bar{\Omega})$ a solution to the boundary value problem \begin{align*} \Delta u &= f(u,\nabla u)\quad\text{ in }\Omega\,,\\ u=\text{const}\,,& \quad\nabla u\neq0\,,\quad\text{ on }\partial\Omega\,, \end{align*} where $f\in C^1$, $f_u\geq0$. Then $u$ has exactly one critical point in $\bar{\Omega}$ and there $\det(D^2u)>0$ holds (i.e., a global maximum or minimum).

The paper claims that Theorem 2 is true even if we replace Laplace's operator by another elliptic operator. This means, for example, that if $f$ is constant, then the solution to (1) has a unique critical point.

Edit 1: Here is the specific PDE I am interested in: \begin{equation} (1+h^2y^2)\partial_{xx}u +(1+h^2x^2)\partial_{yy}u -2xyh^2\partial_{xy}u-\frac{h^2x(3+h^2\rho^2)}{1+h^2\rho^2}\partial_{x}u-\frac{h^2y(3+h^2\rho^2)}{1+h^2\rho^2}\partial_{y}u = \frac{1}{1+h^2\rho^2} \end{equation} where $h\in (0,1]$, $\rho^2=x^2+y^2$ and the domain $\Omega$ is convex, bounded and does not contain the origin. I have verified that this PDE is elliptic for $h\in(0,1]$.

Edit 2: Due to Mateusz's counterexample, I am adding the condition that if $a(x,y)$ is a coefficient of $L$, then on any simply connected subset $D$ of $\Omega$, $a$ attains its maximum and minimum on $\partial D$. The point of this is to avoid any 'humps' in the coefficient functions (since this is the case for my Pde in edit 1).

Answer 1

(This is an extended comment, I suppose).

What do we assume about $L$? In complete generality this cannot be true, as shown by the example that follows. Even if we assume, as in the paper that you linked, that $L$ is in non-divergence form without the first-order term, the result seems to be false. The article gives a reference to a dissertation in German; I suppose further details are given there.

Take $\Omega$ to be a rectangle $[-4, 4] \times [-1, 1]$. Pick an auxiliary bump function $\phi$ which is smooth, non-negative, equal to zero outside $[-2, 2]$ and equal to some large number $\lambda$ in $[-1, 1]$. Define $$ L u(x, y) = \partial_{xx} u(x, y) + (1 + \phi(x)) \partial_{yy} u(x, y). $$ Set $f$ to be equal to $1$ everywhere. Then I am rather convinced that the solution of $L u = f$ will have (at least) two local minima for $\lambda$ large enough.

Here is a probabilistic explanation: $-u(x, y)$ is the mean exit time from $\Omega$ of a diffusion generated by $L$. This diffusion moves as the usual Brownian motion, except that it is sped up by a factor of $(1 + \phi(x))$ in the $y$ direction. It will therefore exit $\Omega$ almost immediately with very high probability when started in $(0, y)$ for $y \times (-1, 1)$. On the other hand, it will usually take some time before it exits $\Omega$ when started at $(\pm 3, 0)$. That will make $-u$ small near the vertical axis, but of constant order near $(\pm 3, 0)$.

If you like, I can try to give a more rigorous (and purely analytic) argument.