"Natural" proof of $P\left(\frac{d}{dx}\right)\bigl(e^{xy}Q(x)\bigr)=Q\left(\frac{d}{dy}\right)\bigl(e^{xy}P(y)\bigr)$. In the context of linear differential equations, I've stumbled upon the following identity for an arbitrary pair of polynomials $P$ and $Q$ with real or complex coefficients:
$$
  P\left(\frac{d}{dx}\right)\bigl(e^{xy}Q(x)\bigr)
  =\sum_{n=0}^\infty\frac{P^{(n)}(y)e^{xy}Q^{(n)}(x)}{n!}
  = Q\left(\frac{d}{dy}\right)\bigl(e^{xy}P(y)\bigr).
$$
This can be more or less easily checked by using Taylor expansions of $P\bigl(\frac{d}{dx}\bigr)$ at $y$ and of $Q\bigl(\frac{d}{dy}\bigr)$ at $x$:
$$
  P\left(\frac{d}{dx}\right)
  =\sum_{n=0}^\infty\frac{P^{(n)}(y)}{n!}\left(\frac{d}{dx} - y\right)^n,
  \quad
  Q\left(\frac{d}{dy}\right)
  =\sum_{n=0}^\infty\frac{Q^{(n)}(x)}{n!}\left(\frac{d}{dy} - x\right)^n.
$$
Is there any easy way to "see" that
$P\bigl(\frac{d}{dx}\bigr)\bigl(e^{xy}Q(x)) = Q\bigl(\frac{d}{dy}\bigr)\bigl(e^{xy}P(y)\bigr)$
without "getting hands dirty"?
Is this identity a part of some general theory?
It makes me think of Fourier analysis, but I do not know much about it.
 A: Here I post a proof based on the comment by KCd since he has not done it himself (yet):
$$
P\left(\frac{\partial}{\partial x}\right)\bigl(e^{xy}Q(x)\bigr)
= P\left(\frac{\partial}{\partial x}\right)Q\left(\frac{\partial}{\partial y}\right)e^{xy}
= Q\left(\frac{\partial}{\partial y}\right)P\left(\frac{\partial}{\partial x}\right)e^{xy}
= Q\left(\frac{\partial}{\partial y}\right)\bigl(e^{xy}P(y)\bigr).
$$
A: The equation is bilinear in $P$ and $Q$ on both sides, thus it is sufficient to show it for $P(x)=x^m$ and $Q(x)=x^n$. Using the Leibniz rule for the derivatives of a product,
\begin{align}
\frac{∂^m}{∂x^m}\left(e^{xy}x^n\right)
&=\sum_{k=0}^m\binom{m}{k}\frac{∂^k}{∂x^k}\left(e^{xy}\right)\frac{∂^{m-k}}{∂x^{m-k}
}\left(x^n\right)
\\
&=\sum_{k=0}^m\binom{m}{k}y^ke^{xy}\cdot n(n-1)\cdots(n-(m-k)+1)x^{n-(m-k)}
\\
&=\sum_{k=0}^mm(m-1)\cdots (k+1)y^ke^{xy}\cdot \binom{n}{m-k}x^{n-(m-k)}
\\
&=\sum_{k=0}^mm(m-1)\cdots (m-k+1)y^{m-k}e^{xy}\cdot \binom{n}{k}x^{n-k}
\\
&=\sum_{k=0}^m\binom{n}{k}\frac{∂^k}{∂y^k}\left(y^m\right)\cdot \frac{∂^{n-k}}{∂y^{n-k}}\left(e^{xy}\right)
\\
&=\frac{∂^n}{∂y^n}\left(e^{xy}y^m\right)
\end{align}
Note that $\binom{m}{k}=0$ for $k>m$ and $\frac{d^k}{dy^k}\left(y^m\right)=0$ for $k>m$, so that the upper bound of the sums can be adjusted as needed.
This establishes the identity in the basis cases and thus also in general.

Or shorter as proposes in a comment by @KCd, refined by OP @Alexey, using that the partial $x$ and $y$ derivatives commute:
$$
\frac{∂^m}{∂x^m}\left(e^{xy}x^n\right)
=\frac{∂^m}{∂x^m}\left(\frac{∂^n}{∂y^n}e^{xy}\right)
=\frac{∂^n}{∂y^n}\left(\frac{∂^m}{∂x^m}e^{xy}\right)
=\frac{∂^n}{∂y^n}\left(e^{xy}y^m\right).
$$
A: You can in fact generalise the identity for arbitrary functions $P,Q$, using Fourier analysis. In general, a function of a derivative is defined as
$$
(f(D)g)(x):=\int_{\mathbb R}f(2\pi ik)\mathrm e^{2\pi ikx}g(k)\mathrm dk
$$
where $g(k)$ is the Fourier transform of $g(x)$:
$$
g(k)=\int_{\mathbb R}\mathrm e^{-2\pi ikx}g(x)\mathrm dx
$$
Now to your identity: the l.h.s. is, by definition
$$
P\left(\frac{d}{dx}\right)\bigl(e^{xy}Q(x)\bigr)\equiv\int_{\mathbb R}P(2\pi ik)\mathrm e^{2\pi ikx}\left[\int_{\mathbb R}\mathrm e^{-2\pi ikx+xy}Q(x)\mathrm dx\right]\mathrm dk
$$
while the r.h.s. is
$$
Q\left(\frac{d}{dy}\right)\bigl(e^{xy}P(y)\bigr)\equiv\int_{\mathbb R}Q(2\pi ik)\mathrm e^{2\pi iky}\left[\int_{\mathbb R}\mathrm e^{-2\pi iky+xy}P(y)\mathrm dx\right]\mathrm dk
$$
Both sides are manifestly equal, under a trivial change of variables.
A: In this answer, I will take your suggestion and articulate it in abstract terms. 

Let $\mathsf X := C^\infty(\mathbb R^2)$ be the algebra of smooth functions on $\mathbb R^2$ w.r.t. pointwise function multiplication $\cdot$, and define the linear operators $\partial_1, \partial_2 : \mathsf X \to \mathsf X$ as the partial derivatives w.r.t. $x$ and $y$ respectively. These two operators commute thanks to Schwarz's lemma:
$$[ \partial_1, \partial_2 ] = \partial_1\circ \partial_2 - \partial_2\circ\partial_1 \equiv 0. $$
Notice that for the specific function $\varphi \in \mathsf X$ given by $\varphi(x,y) = e^{xy}$, these operators act in a likewise specific way:
$$\partial_1 \varphi = \pi_2 \cdot \varphi, \qquad \partial_2 \varphi = \pi_1 \cdot \varphi,$$
where $\pi_i \in \mathsf X$, $i=1,2$, is the projection to the $i$-th entry, so that e.g. $\pi_1(x,y) = x$. By the linearity of the operators $\partial_i$ and the bilinearity of pointwise multiplication, we have that for any choice of real polynomials $P,Q \in \mathbb R[\mathtt x]$ in the indeterminate $\mathtt x$, 
$$P(\partial_1)\varphi = P(\pi_2) \cdot \varphi, \qquad Q(\partial_2)\varphi=Q(\pi_1) \cdot \varphi, $$
(with the obvious understanding that powers of the differential operators are to be intended in terms of operator composition, while powers of elements of $\mathsf X$ must be viewed in terms of pointwise multiplication).
The fact that $\partial_1\circ \partial_2 = \partial_2\circ \partial_1$ yields, by linearity of $\partial_1$, that $P(\partial_1) \circ \partial_2 = \partial_2 \circ P(\partial_1)$ for any arbitrary $P \in \mathbb R[\mathtt x]$; similarly for $Q \in \mathbb R[\mathtt x]$ we find, by linearity of $\partial_2$, that $\partial_1 \circ Q(\partial_2) = Q(\partial_2) \circ \partial_1$. These two facts together imply that $P(\partial_1)$ and $Q(\partial_2)$ commute too:
$$[P(\partial_1), Q(\partial_2)] = P(\partial_1) \circ Q(\partial_2) - Q(\partial_2)\circ P(\partial_1) \equiv 0 .$$
Hence we may conclude the thesis:
$$P(\partial_1) \big( Q(\pi_1) \cdot \varphi \big) = P(\partial_1) \big( Q(\partial_2)\varphi \big) = Q(\partial_2) \big( P(\partial_1) \varphi) = Q(\partial_2) \big( P(\pi_2) \cdot \varphi).  $$
