Let $\Omega\subset \mathbb{R}^d$: open, connected.
Suppose $u\in L^2(\Omega)$ and $\Delta^n u\in L^2(\Omega)$ for $n=1,\dotsc,N$, where $\Delta$ is in the weak sense.
Let $\zeta\in C_c^{\infty}(\Omega)$.
Can we say that $\Delta^N(u\zeta)\in L^2(\Omega)$?
We know if $v\in H^2(\Omega)$ and $\zeta\in C_c^{\infty}(\Omega)$, then $v\zeta\in H^2(\Omega)$, but here we don't know if $u$ is in Sobolev (we don't know how $u$ would behave near the boundary).
I first considered the case $N=1$, $\Delta (u\zeta)$, used the product rule for the Laplacian assuming $u$ is smooth, but got cross term $u_{xy}\zeta_{xy}$, $u_{yz}\zeta_{yz}$, $u_{zx}\zeta_{zx}$ and now I am not sure if it holds because we don't really know if they are square integrable. (Do we?)
But we have $\Delta^n u\in L^2(\Omega)$ and just multiplying a smooth function. Is $\Delta^N(u\zeta)\in L^2(\Omega)$ too good to be true?
If the weak Laplacian is in $L^2$, then the function is in $H^2$ locally. Indeed, the Poisson equation with $L^2$ data admits an $H^2_{\rm loc}$ solution $w$ (e.g., here. Then the difference $u-w$ satisfies the Laplace equation in the weak sense, and therefore (by Weyl's Lemma) is a classical harmonic function, in particular $C^\infty$. (I think you can also find the above result in PDE by Evans somewhere under "interior regularity for Poisson's equation".)
Applying the above repeatedly, we get $u\in H^{2N}$ locally. The boundary behavior is irrelevant because you multiply by a compactly supported $\zeta$. So, $\zeta u\in H^{2N}(\Omega)$ under your assumptions. In particular, $\Delta^N(\zeta u)\in L^2(\Omega)$.