Let's say that I am trying to find the solution to a complicated linear PDE on an infinite 3D domain for a scalar field $\Phi$. $$ \mathcal{L}_1[\Phi]=0 $$ I can integrate over the volume on both sides to get: $$ \int \mathcal{L}_1[\Phi] dV=0. $$ Am I now allowed to add a total divergence inside the integral and manipulate the differential operator such as: $$ \int \mathcal{L}_1[\Phi] +\nabla\cdot\mathcal{F}[\Phi]dV=\int \mathcal{L}_2[\Phi]dV=0 $$ to then solve the resulting PDE $$ \mathcal{L}_2[\Phi]=0? $$ where $\mathcal{L}_2$ is a different linear differential operator.
Many thanks.
No, you are not allowed to do this.
It's easiest to illustrate using 1D examples, where every function $f$ is the "divergence" of some $F$, i.e. $f = F'$. If you could do what you are suggesting, this would mean $$ \int_{\mathbb{R}} f = \int_{\mathbb{R}} F' = 0 $$ for every $f$, which we know is false.
To see when you are allowed to do this, notice that it is true that $$ \int_a^b f = F(b) - F(a), $$ so if $F(b) - F(-b) \rightarrow 0$ as $ b \rightarrow \infty$, we are good to go ... except not really. If e.g. $f(x) = x$, then $\int_{-b}^b f = 0$, but the integral is not absolutely convergent, i.e. the convergence depends on a choice of $a = -b$ above, and other choices lead to different things. So typically the requirements you will find to do this will ask that $$ \lim_{R\rightarrow \infty} \int_{\partial B_R} |F| = 0 $$ and $\int |f| < \infty$, or something along those lines (for $f = \text{div}F$).
It's also worth noting that your last line doesn't follow unless e.g. $\mathcal{L}'[\Phi]\geq 0$. You can do your computation against test functions instead, which also avoids this integrability issue.