$
\newcommand{\F}{\mathbb{F}}
\newcommand{\C}{\mathbb{C}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\Sp}{\mathrm{Sp}}
\newcommand{\Cl}{\mathrm{Cl}}
\newcommand{\P}{\mathcal{P}}
$Here's the short version: The Clifford group is a semidirect product if and only if $p$ is an odd prime.
Now for the details. Unfortunately, they will be still very sketchy as a more or less proper treatment would easily fill a few pages. I will add some references below, if you're interested.
Let us first discuss the case when $p$ is an odd prime which is arguably the simpler and mathematically better behaved case.
The Weil representation
For odd characteristic $p$, we have the metaplectic or Weil representation $\mu$ of $\Sp_{2n}(p)$ on $(\C^{p})^{\otimes n}$. Its existence follows from the Stone-von Neumann theorem on the representations of the Heisenberg group (see Ref. 1-2). It is defined by requiring that
$$
\mu(g) W(v) \mu(g)^\dagger = W(gv), \qquad \forall g\in\Sp_{2n}(p),
$$
where $W(v)$ for $v\in\F_p^{2n}$ are the Weyl operators (or generalized Pauli operators, clock and shift operators ...). For completeness, they are defined as
$$
W(z,x) := \omega^{-2^{-1}z\cdot x} Z(z) X(x),
$$
where $\omega = \exp(2\pi i / p)$, $X(x)|y\rangle := |y+x\rangle$, $Z(z)|y\rangle := \omega^{y\cdot x}|y\rangle$, and $x,y,z,\in\F_p^n$.
In general, the above equation only defines a projective representation. For finite fields, however, it can be linearized (not over the reals $\mathbb R$, though!). See Ref. 2 for an explicit construction.
Clearly, every $\mu(g)$ is by definition a Clifford unitary and thus $\mu(\Sp_{2n}(p))$ is a subgroup.
Since we have $\Cl_n(p) / \P_n(p) \simeq \Sp_{2n}(p)$ (at least for equal centres $Z(\Cl_n(p)) = Z(\P_n(p))$), we can conclude that
$$
\Cl_n(p) = \P_n(p) \rtimes \mu(\Sp_{2n}(p)).
$$
The subgroup $\mu(\Sp_{2n}(p))$ is simply the group of "proper" Clifford unitaries, i.e. the qudit versions of the usual Clifford gates such as Hadamard, phase, and controlled-$Z$ gate. In fact, those generate this subgroup. Contrary to the qubit case, the phase gate does not generate any $Z$ operators (at least when properly defined).
The qubit case
The characteristic two case is however more complicated (as always). First of all, the Weil representation does not exist, not even projectively. Of course, that doesn't directly imply that the answer is negative. However, assuming that there is a subgroup of $\Cl_n(2)$ isomorphic to $\Sp_{2n}(2)$ means that there is a representation of $\Sp_{2n}(2)$. One can then use the same particularities that prevent the existence of the Weil representation to derive a contradiction. This has been explicitly shown in Thm. 7 of Ref. 3 (the "Clifford transform group" is what we call "Clifford group" today).
Interestingly, the non-existence of such a factorization is closely related to the fact that no Clifford covariant Wigner function exists for qubit systems, cp. Ref. 4.
So what is the qubit Clifford group, if not a semidirect product?
Let us consider the action of a Clifford unitary $U$ on a Pauli operator $W(v)$. We can write it as
$$
U W(v) U^\dagger = (-1)^{\alpha(v)} W(gv),
$$
where $g\in\Sp_{2n}(2)$ and $\alpha:\, \F_2^{2n}\rightarrow \F_2$ is some phase function.
We can derive some conditions on $\alpha$ by considering the multiplication law of Pauli operators:
$$
W(v)W(u) = i^{\beta(v,u)}W(v+u).
$$
Here $\beta$ is a $\Z_4$-valued function which is given as
$$
\beta(v,u) = (v+u)_z \cdot (v+u)_x - v_z\cdot v_x - u_z\cdot u_x + 2(v_z\cdot u_x) \mod 4.
$$
Then $\alpha$ has to fulfill the following constraint
$$
2\alpha(u+v) = 2\alpha(u) + 2\alpha(v) + \beta(gu, gv) - \beta(u,v) \mod 4.
$$
Hence, we can choose the values of $\alpha$ freely on a basis $e_i$ of $\F_2^{2n}$ (i.e. on generators of the Pauli group), and the above constraint then determines the values on all other elements. I think it is clear that every consistent choice $(g,\alpha)$ defines a Clifford unitary (up to a global phase) and vice versa. However, any such choice of values determines a unique $v\in\F_2^{2n}$ by
$$
\alpha(e_i) = \begin{cases} v_{n+i} & i \leq n, \\ v_i & i > n. \end{cases}
$$
Let $\alpha_0$ be the choice $\alpha_0(e_i)=0$. Then, the Clifford $U$ defined by $(g,\alpha)$ and $U_0$ defined by $(g,\alpha_0)$ are related by
$$
U = U_0 W(v),
$$
again, up to a global phase.
Finally, the projective Clifford group is isomorphic to this set of pairs with the multiplication law
$$
(g,\alpha_g)\bullet (h,\alpha_h) = (g\circ h, \alpha_g\circ h + \alpha_h).
$$
If we want to label the Cliffords by a $g \in \Sp_{2n}(2)$ and a vector $v\in\F_2^{2n}$, we can write this as
$$
(g,v) \bullet (h,u) = (g\circ h, \sum_{i=1}^n (\alpha_g(h e_i) e_{n+i} + \alpha_g(h e_{n+i}) e_i) + u),
$$
where $\alpha_g$ has to be evaluated according to the above constraint (e.g. recursively).
So why is odd $p$ so much easier?
For odd $p$, $\beta$ is a $\F_p$-valued bilinear form which is invariant under $\Sp_{2n}(p)$. Thus, the constraint becomes $\alpha(u+v) = \alpha(u) + \alpha(v)$, i.e. $\alpha$ is independent of $g$ and, moreover, linear.
References
- Gross: "Hudson's Theorem for finite-dimensional quantum systems" (2006)
- Neuhauser: "An Explicit Construction of the Metaplectic Representation over a Finite Field" (2002)
- Bolt, Room, and Wall: "On the Clifford collineation, transform and similarity groups. II." (1960)
- Zhu: "Permutation Symmetry Determines the Discrete Wigner Function" (2016)
- [Gottesman: "Fault-Tolerant Quantum Computation with Higher-Dimensional Systems" (1999)][5]
Remark: Everything in this answer continues to hold if we replace the finite field $\F_p$ by some finite field of characteristic $p$, say $\F_{p^k}$. In practice, this case is not so relevant, thus I omitted it here to keep the answer as simple as possible.
