I have thought about this apparently simple problem.
Question: Let $\Omega\subset \mathbb{R}_{x,y}^2$ be an open subset. Let's suppose we have $u\in\mathscr{D}'(\Omega)$ that satisfies $$ (\partial_y-\partial_x^2)u=0.$$ Is it true that $u\in C^\infty(\Omega)$?
As a small comment, the same result for the Laplace equation $\Delta u=0$ is true, and this is known as Weyl's Lemma.
So, I am posting an answer with what I have learnt. The answer to the question is yes. The heat operator $L=\partial_t-\Delta_x$ in $\mathbb R^{n+1}$ is in fact hypoelliptic, that is, for any given $u\in \mathscr D’(\Omega)$ and $U\subset \Omega$ open, $$ Lu\in C^\infty(U)\implies u\in C^\infty(U).$$ In other words, a distributional solution to the equation $Lu=f$ is smooth (in space-time) in all the open sets where $f$ is smooth. It is not totally immediate to see, but just to give an idea, this essentially comes from the fact that the fundamental solution of the heat equation $$K(x,t)=\frac{1}{(4\pi t)^\frac{n}{2}}e^{-\frac{|x|^2}{4t}}$$ is smooth away from the origin $(x,t)=(0,0)$. In fact, this condition (having a smooth fundamental solution away from the origin) for constant-coefficient differential operators is equivalent to hypoellipticity.
A good starting point to learn a bit more and see the proofs of these things is Chapter 1, Section F of the book “Introduction to partial differential equations” by G. B. Folland. This section is self-contained and explains essentially all one strictly needs to understand. An equally good or even better reference is the book by Hörmander, the same reference that can be found in the Wikipedia page of Weyl’s lemma about hypoellipticity:
Lars Hörmander, The Analysis of Linear Partial Differential Operators I, 2nd ed., Springer-Verlag (1990), §4.4.
Of course, from this point one can try to generalize this result, for instance thinking about the non-constant coefficient case, or proving regularity results in some (local or global) Sobolev spaces, as well as one does for elliptic operators. Both books are dealing with these subjects (Hörmander’s book is part of a series of books, possibly one has to look at the second book). I found particularly interesting, to have an idea on which results one could expect in the context of the Sobolev spaces $H^s(\mathbb R^n)$, an exercise in another Folland’s book:
G. B. Folland, Real Analysis. Modern techniques and their applications, 2nd ed., (1999), §9.3, exercise 37.
37. The heat operator $\partial_{t}-\Delta$ is not elliptic, but a weakened version of elliptic regularity holds for it. Here we are working on $\mathbb{R}^{n+1}$ with coordinates $(x, t)$ and dual coordinates $(\xi, \tau)$, and $\partial_{t}-\Delta=P(D)$ where $P(\xi, \tau)=2 \pi i \tau+4 \pi^{2}|\xi|^{2}$.
There exist $C, R>0$ such that $|\xi||(\xi, \tau)|^{1 / 2} \leq C|P(\xi, \tau)|$ for $|(\xi, \tau)|>R$. (Consider the regions $|\tau| \leq|\xi|^{2}$ and $|\tau| \geq|\xi|^{2}$ separately.)
If $f \in H_{s}$ and $\left(\partial_{t}-\Delta\right) f \in H_{s}$, then $f \in H_{s+1}$ and $\partial_{x_{i}} f \in H_{s+(1 / 2)}$ for $1 \leq i \leq n .$
If $\zeta \in C_{c}^{\infty}\left(\mathbb{R}^{n+1}\right)$, we have $$ \left[\partial_{t}-\Delta, \zeta\right] f=\left(\partial_{t} \zeta-\Delta \zeta\right) f-2 \sum\left(\partial_{x_{i}} \zeta\right)\left(\partial_{x_{i}} f\right) $$
If $\Omega$ is open in $\mathbb{R}^{n+1}$, $u \in \mathcal{D}^{\prime}(\Omega)$, and $\left(\partial_{t}-\Delta\right) u \in H_{s}^{\text {loc }}(\Omega)$, then $u \in$ $H_{s+1}^{\text {loc }}(\Omega)$.