Advection-diffusion with periodic boundary conditions and tilt Context: Consider the advection-diffusion equation with periodic boundary conditions (PBC) over a flat square domain $L \times L$.
The scalar density $\rho $ is transported by a prescribed field $\mathbf{v}=-\nabla U$, where $U(\mathbf{x})$ is a scalar potential that has the periodicity imposed by the PBC. The density $\rho$ evolves as
$$
\partial_t \rho(\mathbf{x},t) = -\nabla \cdot [ \mathbf{v}(\mathbf{x}) \rho(\mathbf{x},t) -   \nabla \rho(\mathbf{x},t) ] = 0
$$
The steady-state solution is found by imposing $\partial_t \rho(\mathbf{x},t) =0$ and has the usual Gibbs form:
$$
\rho(\mathbf{x}) \,  \propto  \, e^{-U(\mathbf{x}) } 
$$
The problem: I am wondering how to find the steady-state in a slightly more general case, where
$$\mathbf{v} = -\nabla U + \mathbf{q}$$
The potential $U$ has the periodicity imposed by the PBC and $\mathbf{q} =(q_x,q_y)$ is a constant vector field (the constant \mathbf{q} defines the so-called "tilt" of the potential $U$).
Hence, the equation we have to solve is
$$
\nabla \cdot [ \, \rho(x,y) \, \mathbf{q} -  \rho(x,y) \nabla U(x,y) - \nabla \rho(x,y) \, ] = 0
$$
with the periodic conditions
$\rho(0,y) = \rho(L,y)$,
$\rho(x,0) = \rho(x,L)$,
$U(0,y) = U(L,y)$,
$U(x,0) = U(x,L)$. For simplicity, I tried to consider the case $\mathbf{q}=(q,0)$, but the problem still seems non-trivial.
Question: Any idea or reference about the diffusion-advection equation in periodic boundary conditions (in particular about the steady-state)?  Which is the "Gibbs-like solution" in this case?
Further considerations: I have the feeling that finding a solution is not easy because the "tilt" potential that generates the constant field $\mathbf{q}$ is $-\mathbf{x}\cdot \mathbf{q}$. This "tilt" contribution males the total potential $U-\mathbf{x}\cdot \mathbf{q}$ not periodic (i.e. it does not satisfy the PBC).
Moreover, define the total current in the steady-state as
$$
\mathbf{J}(x,y) = \rho(x,y) \, [\mathbf{q} -  \nabla U(x,y)] - \nabla \rho(x,y) \, ,
$$
so that we have to find the $\mathbf{J}$ such that
$$
\nabla \cdot \mathbf{J} = 0  \quad \Rightarrow \quad \mathbf{J} = R \nabla g
$$
where $R$ is a 90-degrees rotation and $g$ is an unknown scalar potential. Note that $g$ does not have to respect the PBC, but $\mathbf{J}$ does: (probably) the most general form of $g$ is
$$
g(x,y) = G(x,y) + a x + b y
$$
where  $G$ respects the PBC and $a$ and $b$ are constants. Even though this problem is more likely to be studied by physicists, I have the feeling that the problem is intimately related to the topology of the 2D torus, so I posted also a similar question on math SE.
 A: Any $\rho$ which solves the equation on the whole torus must also be a solution locally on every subset. In particular, it must be solution on the (non-toroidal) open $L \times L $ square. Since solutions on the torus are a subset of the solutions on the square, the question becomes: Do there exist solutions on the square which happen to match at the boundaries?
On this square, we can define $V = U - \mathbb{x} \cdot \mathbb{q}$, and we have an ordinary advection-diffusion equation. We know there exist solutions of the form $\alpha e^{-V(\mathbb{x})}$. We also know that $U$ is periodic, so $V$ can only be periodic if $\mathbb{q} = \mathbb{0}$. However $e^{-V}$ could still be periodic if $\mathbb{q}$ is imaginary. Specifically, we have periodic solutions for $\mathbb{q} = \frac{2\pi i}{L}\mathbb{n}, \mathbb{n} \in \mathbb{Z}^2$.
For other $\mathbb{q}$, solutions proportional to $ e^{-V(\mathbb{x})} $ cannot extend to solutions on the whole torus. The remaining question: Are such solutions the whole solution space?
Now,  Matthew Kvalheim points to Zeeman, 1988. Theorem 3 reads

Let $U$ be a vector field on a compact manifold $X$ without boundary, and
let $\epsilon$ > 0. Then the Fokker-Planck equation for $U$ with $\epsilon$-diffusion has
a unique steady state, and all solutions tend to that steady state.

The torus is a compact manifold without boundary, Zeeman's $U$ is our $-\nabla V$, and we have $\epsilon = 1$, so the theorem tells us a solution $\rho$ must exist and is unique (up to an overall scalar). Unfortunately, this proof is not constructive.
In one dimension, variation of parameters gives the solution
$$\rho = C_1 e^{-V}\left(C_2 + \int_0^x e^V\right)$$
and the requirement $\rho(0) = \rho(L)$ fixes $C_2$. We can try to extend this to two dimensions as follows:
Assume $\rho$ is of the form $\alpha(x)e^{-V}$. Then the equation becomes
$$ \nabla \cdot [-\nabla V \alpha(x)e^{-V} - \nabla (\alpha(x) e^{-V})] = 0 $$
which simplifies to
$$ \nabla \cdot (\nabla\alpha(x) e^{-V}) = 0 $$
Solutions are
$$ \nabla\alpha(x) e^{-V} = \nabla \times \mathbf{\psi} $$
for $\mathbf{\psi} = \mathbf{e}_z $ and $g$ some scalar function. Then
$$ \nabla\alpha(x)  = e^V(\nabla \times \mathbf{\psi}) $$
If
$$ \nabla \times (e^V(\nabla \times \mathbf{\psi})) = 0 $$
then this has solution
$$ \alpha(x,y) = C + \left(\int_0^x -e^V g_y dx\right) + \left(\int_0^y e^V g_x dy\right) $$
The requirement of periodic boundary conditions picks out some unique $g$, $C$ up to an overall constant. We need
$$ \alpha(x,0) = \alpha(x,L)e^{-Lq_y} $$
or
$$ C + \left(\int_0^x -e^V g_y dx\right) = Ce^{-Lq_y} + \left(\int_0^x -e^V g_y dx\right)e^{-Lq_y} + \left(\int_0^L e^V g_x dy\right)e^{-Lq_y} $$
At $x = 0$ this simplifies to
$$ C =  \frac{1}{e^{Lq_y} - 1}\int_0^L e^V g_x(0,y) dy $$
It remains to find $g$.
I'm not sure that there is a nice expression for the solution in general. Some miscellaneous thoughts:

*

*When $U = 0$, $\rho = C$ is a solution, which corresponds to $\alpha = e^V, g = xq_y-yq_x$. This shows that $g$ may be defined only on the square, not on the torus.

*When $\nabla U \gg q$ or $q \gg \nabla U $, we can start with the nearby known solution and series expand.

