  We present an $L_q(L_{p})$-theory for the equation $$
\partial_{t}^{\alpha}u=\phi(\Delta) u +f, \quad t>0,\, x\in \mathbb{R}^d
\quad\, ;\, u(0,\cdot)=u_0. $$ Here $p,q>1$, $\alpha\in (0,1)$,
$\partial_{t}^{\alpha}$ is the Caputo fractional derivative of order $\alpha$,
and $\phi$ is a Bernstein function satisfying the following: $\exists
\delta_0\in (0,1]$ and $c>0$ such that \begin{equation}
  \label{eqn 8.17.1} c \left(\frac{R}{r}\right)^{\delta_0}\leq
\frac{\phi(R)}{\phi(r)}, \qquad 0<r<R<\infty. \end{equation} We prove
uniqueness and existence results in Sobolev spaces, and obtain maximal
regularity results of the solution. In particular, we prove \begin{align*} \|
|\partial^{\alpha}_t u|+|u|+|\phi(\Delta)u|\|_{L_q([0,T];L_p)}\leq
N(\|f\|_{L_q([0,T];L_p)}+ \|u_0\|_{B_{p,q}^{\phi,2-2/ \alpha q}}), \end{align*}
where $B_{p,q}^{\phi,2-2/\alpha q}$ is a modified Besov space on $\mathbb{R}^d$
related to $\phi$.
  Our approach is based on BMO estimate for $p=q$ and vector-valued
Calder\'on-Zygmund theorem for $p\neq q$. The Littlewood-Paley theory is also
used to treat the non-zero initial data problem. Our proofs rely on the
derivative estimates of the fundamental solution, which are obtained in this
article based on the probability theory.
