Question: "First question: why isn't this defined simply as $p^∗_1M⊗_{O_{X×X}}p^∗_2N$, where a differential operator acts via the Leibniz rule? Second question: why? I tried to prove this with no success."
Answer: "In the Borel book "Algebraic D-modules" there is a section where they construct direct images, inverse images etc for holonomic modules on the Weyl algebra $W_n$. The constructions are explicit and elementary and it should not be too difficult to write down explicit examples for such modules on $W_n$ - have you done this? Then you can check formulas $F1,F2$ explicitly."
In the lecure on the Weyl algebra they define the inverse image $P^*M$ and direct image $P_*M$ of a holonomic module $M$ under polynomial mappings between affine spaces
and prove that $P^*M,P_*M$ are again holonomic. The direct image and inverse image functors $P^*,P_*$ are defined in a purely algebraic manner - there is no "topological inverse image". When working with maps of schemes (or varieties) $f:X \rightarrow Y$ there are two notions: For any $E \in Qcoh(Y)$ there is the topological inverse image $f^{-1}(E)$ of $E$ and the pull back as quasi coherent sheaf $f^*(E)$. By definition
$$f^*(E):= \mathcal{O}_X \otimes_{f^{-1}(\mathcal{O}_Y)} f^{-1}(E).$$
It seems from formula $F1$ ou are using the topological inverse image $p_i^{-1}$. You should clearify what you mean in formula $F1$ and relate it to Borel's book. Maybe this will give more insight. The topological inverse image is not used when studying modules on sheaves of rings of differential operators. A left $D_X$-module $E$ is a left $\mathcal{O}_X$-module equipped with the structure of a $D_X$-module.
Example: When $X$ is regular of finite type over the complex numbers, a $D_X$-module is a quasi-coherent $\mathcal{O}_X$-module $E$ equipped with a flat connection $\nabla$. This is beacuse $D_X$ is locally generated by derivations. When pulling back a $D_X$-module $(E,\nabla)$, you first pull back the $\mathcal{O}_X$-modules $E$ and then you pull back the connection. When you pull back $E$ you use the construction $f^*(E)$ and not the topological inverse image $f^{-1}(E)$.
Example: Formula $F2$ (if you ignore the action of the (ring of) differential operators) says for $A$-modules the following:
Let $p,q:A\rightarrow A\otimes_k A\rightarrow^m A$ where $p,q$ are the caononical maps and $m$ the multiplication map and let $E,F$ be $A$-modules.
It follows
$$ E\boxtimes F \cong (E\otimes_k A)\otimes_{A\otimes_k A}(A \otimes_k F).$$
There are obvious maps between $E\otimes_A F$ and $\Delta_A^*(E\boxtimes F):=(E\boxtimes F)\otimes_{A\otimes_k A}A$. There is one map
$$\rho: E\otimes_A F \rightarrow \Delta_A^*(E\boxtimes F)$$
defined by
$$\rho(x \otimes y):=x \otimes 1 \otimes 1 \otimes y \otimes 1.$$
Use the definitions in [Borel] and check if similar maps exist for holonomic $W_n$-modules.
Question: "The second question boils down to the fact that the inverse image commutes with tensor products, which is well known. As for the first question, I have no calculations. My point is that most books define this object in a very convoluted way and there seems to be a natural definition. I wonder if the natural definition works or not."
Answer: The construction of holonomic modules on the Weyl algebra is elementary and explicit and the pull back and push forward $P^*M,P_*M$ of a holonomic module $M$ are defined explicitly in the Borel book without the use of "sheaf theoretic operations". Before trying to understand the global case you should make sure you understand the case of the Weyl algebra and the affine case. To understand the affine case you will need the PBW theorem mentioned earlier on this site.
Example: If $k$ is a field of characteristic zero and $A:=k[x_1,..,x_n], L:=Der_k(A), X:=Spec(A)$ it follows
$$PBW.\text{ }D_X \cong U(A,L) \cong A\{\partial_{x_1}^{l_1}\cdots \partial_{x_n}^{l_n}: l_i \geq 0\},$$
where $\partial_{x_1},..,\partial_{x_n}$ is partial derivative wrto the $x_i$-variable. The right hand side in $PBW$ means the following: The associative ring $U(A,L)$ is a free left $A$-module on the set
$$B:=\{\partial_{x_1}^{l_1}\cdots \partial_{x_n}^{l_n}: l_i \geq 0\}.$$
Hence any differential operator $D\in Diff_k(A)$ may be written uniquely as
a sum
$$D:=\sum_{I:=(i_1,..,i_n)} a_I \partial_{x_1}^{i_1}\cdots \partial_{x_n}^{i_n}$$
with $a_I \in A$ for all $I$. The ring $D_X\cong U(A,L)$ is isomorphic to the Weyl algebra on $n$ variables. If $i:=i_1+\cdots i_n$ and $\partial:=\partial_{x_1}^{i_1}\cdots \partial_{x_n}^{i_n}$ it follows $\partial \in Diff^i_k(A)$. You get a filtration $D_X^i$ on $D_X$ and $gr(D_X) \cong k[x_i,\partial_{x_i}]$ is isomorphic to a polynomial ring in $2n$ variables. Hence the PBW theorem for $D_X$ says there is an isomorphism of graded $A$-algebras
$$Sym_A^*(L) \cong gr(D_X).$$
Since $L$ is free on the elements $\partial_{x_i}$ it follows $gr(D_X)\cong Sym_A^*(L)$ is a polynomial algebra over $A$ on the element $\partial_{x_i}$.
Note: Some people may try to convince you that you can prove important results in the theory of $D_X$-modules by "diagram chase" and category theory. This is not true in general.
