Is it true that for each real $p\ge1$ there is some real $C_p$ such that for all smooth real-valued functions $u$ compactly supported on $S:=(0,1)^2$ one has $$\|D_1D_2u\|_p\le C_p(\|D_1^2u\|_p+\|D_2^2u\|_p),$$ where $D_j$ denotes the operator of the partial differentiation with respect to the $j$th argument ($j=1,2$) and $\|\cdot\|_p$ is the $L^p(S)$ norm?
The special case when $p=2$ is follows easily from the Plancherel isometry.
Since you tagged reference-request:
In the PDE/harmonic analysis literature this is a consequence of the Calderon-Zygmund Inequality, it is one of the main tools for studying elliptic regularity theory in $L^p$. It originated in Calderon and Zygmund's 1952 paper (see Chapter 3; given the age of the paper it is not surprising they state it in a way that is harder to parse using modern language, and the emphasis is on something different so the inequality you seek has to be chained together from several of the results proven). The $L^2$ version is in fact slightly earlier, due to Mikhlin (1948).
For a more modern summary, this is presented in Chapter 9 (see Corollary 9.10) in Gilbarg and Trudinger's Elliptic PDE of second order. Specifically one has the estimate that for any open bounded domain $\Omega\subseteq \mathbb{R}^n$ and any $p\in (1,\infty)$ there exists a constant (dependent on $\Omega$) $C_p$ (with $C_2= 1$) such that any $u\in W^{2,p}_0(\Omega)$ satisfies $$ \| D^2 u\|_p \leq C_p \|\Delta u \|_p $$ (here $\Delta$ is the Laplacian).
Note that when the CZ inequality holds we have that $$ \sum_i \| D_i D_i u\|_p \leq \|D^2 u\|_p \leq C_p \| \Delta u\|_p \leq C_p \sum_i \|D_i D_i u\|_p $$ and so all three quantities are mutually comparable.
For the endpoint cases $p = 1$ and $\infty$, it turns out that the weaker inequality $$ \| D^2 u\| \lesssim \sum_i \|D_i D_i u\| $$ already fails, and hence the CZ-inequality fails also.
For the case $p= \infty$, this can be shown following this answer (which just rehashes a well-known exercise in Gilbarg and Trudinger); taking $P(x_1, x_2) = x_1x_2$ as suggested, the construction gives a function $u$ that is smooth away from the origin, whose second partials $\partial^2_{11}u$ and $\partial^2_{22}u$ extend continuously over the origin, but whose second partial $\partial^2_{12}u$ blows up near the origin.
For the case $p = 1$, this was a result of Ornstein, which says:
Thm Given any set of linearly independent multi-indices in $n$-variables $\alpha_0, \ldots, \alpha_m$, with $|\alpha_i| = |\alpha_j|$, for any $K > 0$ one can find a smooth function $u$ on $\mathbb{R}^n$, vanishing outside $[-1,1]^n$ such that $\| D^{\alpha_0} u\|_{L^1} \geq K$ while $\|D^{\alpha_j} u\|_{L^1} \leq 1$ for all $i > 0$.
In the paper he provides a as an example a fairly simple construction in the case of $n = 2$, and the differential operators being $D^{\alpha_0} = \partial^2_{xy}$, $D^{\alpha_1} = \partial^2_{xx}$ and $D^{\alpha_2} = \partial^2_{yy}$. (Note however there are quite a few printing errors/typos in the published paper.)
It should be true for $p>1$. For a function $\varphi$ in the Schwartz class it holds that \begin{equation} D_1D_2 \varphi(x) = -R_1 R_2 \Delta \varphi(x), \end{equation} where $R_1, R_2$ are the Riesz transforms (See Grafakos, Classical Fourier analysis, Proposition 5.1.17). Therefore \begin{equation} \Vert D_1 D_2 \varphi \Vert_p = \Vert R_1 R_2 \Delta \varphi \Vert_p \leq C_{p} \Vert \Delta \varphi \Vert_p \leq C_{p} (\Vert D_1^2 \varphi \Vert_p + \Vert D_2^2 \varphi \Vert_p), \end{equation} by the $L^p$ boundedness of the Riesz transforms. For $p=1$ maybe one should look for a counterexample.
