If, $u : \Omega (\underset{\text{open}}{\subset} \mathbb{R}^n) \to \mathbb{R}$ satisfies $\Delta u + v \cdot \nabla u = 0 \tag{1}$
Since, $$\displaystyle \Delta \left(u(x)e^{\frac{1}{2}v\cdot x}\right) = e^{\frac{1}{2}v\cdot x}\Delta u(x) + e^{\frac{1}{2}v\cdot x}v\cdot \nabla u(x) + \frac{|v|^2}{4}e^{\frac{1}{2}v\cdot x}u(x)$$ writing, $\displaystyle w(x) = u(x)e^{\frac{1}{2}v\cdot x}$, the equation $(1)$ can be written as: $$\Delta w(x) - c^2w(x) = 0 \tag{$1^{*}$}$$ where, $\displaystyle c = |v|/2 > 0$.
If, $B(x_0,r) \subset\subset \Omega$, then for a solution with $ w \in C^{2m}(\Omega)$, expanding in Taylor series leads to,
\begin{align*}w(x) = \sum\limits_{0 \le |\alpha| \le 2m} \frac{1}{\alpha !}D^{\alpha}w(x_0)(x - x_0)^{\alpha} + O(|x-x_0|^{2m})\end{align*}
Now, call a multi-index $\alpha = (\alpha_1, \cdots , \alpha_n)$ even if each coordinate is even.
Then using $$\int_{B(x_0,r)} (x-x_0)^{\alpha}\,dx = \begin{cases} \dfrac{\Gamma\left(\frac{\alpha_1 +1}{2}\right) \cdots \Gamma\left(\frac{\alpha_n +1}{2}\right)r^{|\alpha|+n}}{\Gamma\left(\frac{|\alpha| + n}{2}+1\right)} &\text{ when, } \alpha \text{ (even) } \\ 0 &\text{ otherwise }\end{cases} \tag{2}$$ integrating both sides of Taylor expansion we have,
\begin{align}\int_{B(x_0,r)} w(x)\,dx &= \sum\limits_{0 \le |\alpha| \le 2m} \frac{D^{\alpha}w(x_0)}{\alpha !}\int_{B(x_0,r)} (x-x_0)^{\alpha}\,dx + O(r^{2m}) \\&= \sum\limits_{0 \le |\alpha| \le m} \frac{D^{2\alpha}w(x_0)}{(2\alpha) !}\int_{B(x_0,r)} (x-x_0)^{2\alpha}\,dx + O(r^{2m})\tag{3}\\&= \sum\limits_{0 \le |\alpha| \le m} \frac{D^{2\alpha}w(x_0)}{(2\alpha) !}\frac{\left(\alpha-\frac{1}{2}\right)!}{\left(|\alpha|+\frac{n}{2}\right)!}r^{2|\alpha|+n} + O(r^{2m}) \tag{4}\\&= \sum\limits_{0 \le |\alpha| \le m} \frac{D^{2\alpha}w(x_0)}{\alpha !}\frac{\pi^{n/2}r^n}{\left(|\alpha|+\frac{n}{2}\right)!}\left(\frac{r}{2}\right)^{2|\alpha|} + O(r^{2m}) \tag{5}\\&= \sum\limits_{0 \le k \le m} \frac{\Delta^{k} w(x_0)}{k!}\frac{\pi^{n/2}r^n}{\left(k+\frac{n}{2}\right)!}\left(\frac{r}{2}\right)^{2k} + O(r^{2m}) \tag{6}\\&= |B(x_0,r)|\left(\frac{n}{2}\right)!\sum\limits_{0 \le k \le m} \frac{\Delta^{k} w(x_0)}{k!\left(k+\frac{n}{2}\right)!}\left(\frac{r}{2}\right)^{2k} + O(r^{2m}) \tag{7}\end{align}
where, In line $(3)$ the sum only survives over even indices. In line $(4)$ we used the Legendre Duplication identity, $$(2\alpha)! = \prod\limits_{j=1}^{n} \Gamma\left(2\alpha_j+1\right) = \prod\limits_{j=1}^{n} \frac{2^{2\alpha_j}\Gamma\left(\alpha_j+\frac{1}{2}\right)\Gamma\left(\alpha_j+1\right)}{\sqrt{\pi}} = \frac{2^{2|\alpha|}\alpha!\left(\alpha-\frac{1}{2}\right)!}{\pi^{n/2}}$$. In line $(5)$ we used the identity, $\displaystyle \Delta^{k}w = \sum\limits_{|\alpha| = k} \frac{k!}{\alpha!}D^{2\alpha}w$. Finally in line $(6)$ we used the volume of ball $\displaystyle |B(x_0,r)| = \frac{\pi^{n/2}r^n}{(n/2)!}$.
Now, if $w \in C^{\infty} (\Omega)$, from $(1^{*})$ we have, $\displaystyle \Delta^{k}w = c^{2k}w$ in $\Omega$ and hence from $(7)$ we have, $$\frac{1}{|B(x_0,r)|}\int_{B(x_0,r)} w(x)\,dx = w(x_0)\left(\frac{n}{2}\right)!\sum\limits_{k=0}^{\infty} \frac{(cr/2)^{2k}}{k!\left(k+\frac{n}{2}\right)!} = \frac{\left(\frac{n}{2}\right)!I_{n/2}\left(cr\right)}{(cr/2)^{n/2}}w(x_0) \tag{8}$$
where, $I_{n/2}$ is the Modified Bessel Function of first kind.
Now, $u(x) = 1$ also satisfies $(1)$, hence, we may rewrite $(8)$ in terms of the mean value formula,
$$u(x_0)\int_{B(x_0,r)} e^{\frac{1}{2}v.x}\,dx = \int_{B(x_0,r)} u(x)e^{\frac{1}{2}v.x}\,dx \tag{9}$$
Differentiating both sides w.r.t., $r$ we get the corresponding mean value formula for spheres,
$$u(x_0)\int_{\partial B(x_0,r)} e^{\frac{1}{2}v.x}\,dS(x) = \int_{\partial B(x_0,r)} u(x)e^{\frac{1}{2}v.x}\,dS(x)\tag{10}$$
Similar to the Laplace equation the converse to the mean value property should also hold in this case (I'll try and add more details if needed).
