Let $$ \mathcal{L}\equiv\sum_{i,j}^{n}a_{i,j}\frac{\partial^{2}}{\partial x_{i}\,\partial x_{j}}+\sum_{i}^{n}b_{i}\frac{\partial}{\partial x_{i}}+c $$ be uniformly elliptic on $\Omega\equiv\mathbb{R}^{n}\times\left(0,T\right)$. Let $$ \mathcal{M\equiv L}-\frac{\partial}{\partial t} $$ Consider the parabolic equality $$ \mathcal{M}u=f\text{ on }\Omega. $$ Lastly, suppose $u\left(\mathbf{x},0\right)$ is monotone increasing in $\mathbf{x}$ and that $f\left(\mathbf{x},t\right)$ is monotone decreasing in $\mathbf{x}$ for all times $t$. A function $g:\mathbb{R}^n \rightarrow \mathbb{R}$ is monotone increasing if for all $\mathbf{x},\mathbf{y}\in\mathbb{R}^n$, $g\left(\mathbf{x}\right)\leq g\left(\mathbf{y}\right)$ whenever $\mathbf{x}\leq \mathbf{y}$. Monotone decreasing is defined similarly.
First, consider the case of $a_{i,j}$, $b_{i}$, $c$ being functions of $t$ only. Define $v\left(\mathbf{x},t\right)=u\left(\mathbf{x}-\epsilon,t\right)$ (where $\epsilon \in \mathbb{R}^n$ with $\epsilon > 0$) and note that $$ \mathcal{M}v=f\left(\mathbf{x}-\epsilon,t\right)\text{ on }\Omega. $$ Let $w=u-v$ so that $$ \mathcal{M}w=f\left(\mathbf{x},t\right)-f\left(\mathbf{x}-\epsilon,t\right)\leq0\text{ on }\Omega. $$ Since $u$ is monotone at time zero, we have $$w\left(\mathbf{x},0\right)=u\left(\mathbf{x},0\right)-v\left(\mathbf{x},0\right)=u\left(\mathbf{x},0\right)-u\left(\mathbf{x}-\epsilon,0\right)\geq0 \text{ on } \mathbb{R}^n.$$ These last two inequalities yield $w\geq0$ (and hence $u$ is monotone in $\mathbf{x}$) everywhere in $\Omega$.
However, this argument is rather restrictive, as it requires the assumption
that the coefficients $a_{i,j}$, $b_{i}$, and $c$ are constant in space (i.e. w.r.t. $\mathbf{x}$). Does anyone know of any relaxations to the above? (see edit below)
EDIT: Some things were forgotten that are important: We need to assume the boundedness of the coefficients in the operator, along with $$\liminf_{\left|\mathbf{x}\right|\rightarrow\infty} \mathbf{u} \geq 0$$ for all times.
EDIT: I have found a potentially interesting relaxation. Suppose $\mathbf{u}$ is convex in $\mathbf{x}$ on $\Omega$. Let $A\in\mathbb{R}^{n\times n}$ be the matrix with $a_{j,k}$ in the $j^{\text{th}}$ row, $k^{\text{th}}$ column. Since $\mathcal{L}$ is assumed to be uniformly elliptic, $A$ is positive semidefinite. Denote this by $A\succeq0$. Now, relax the requirement on the coefficients $a_{j,k}$ as follows: $A\left(\mathbf{x}+\epsilon,t\right)-A\left(\mathbf{x},t\right)\succeq0$ on $\Omega$ for all $\epsilon>0$. Note that in the one-dimensional case (i.e. $n=1$), this requirement becomes $a\equiv a_{1,1}$ is monotone increasing in space. Noting that $$ \left(A\left(\mathbf{x}+\epsilon,t\right)-A\left(\mathbf{x},t\right)\right)\circ\left[\nabla_{\mathbf{x}}u_{i}\right]\left(\mathbf{x},t\right)\succeq0 $$ where $\nabla$ denotes the Hessian operator and $\circ$ the Hadamard (entrywise) product, we get $$ \sum_{i,j}^{n}a_{i,j}\left(\mathbf{x},t\right)\frac{\partial^{2}u}{\partial x_{i}\,\partial x_{j}}\leq\sum_{i,j}^{n}a_{i,j}\left(\mathbf{x}+\epsilon,t\right)\frac{\partial^{2}u}{\partial x_{i}\,\partial x_{j}}. $$ Using the above, we arrive once again at $\mathcal{M}w\leq 0$. This brings me to a more interesting question...
Are there better relaxations than the above when $u\left(\mathbf{x},t\right)$ is convex in $\mathbf{x}$?