Unfortunately, this is one of those cases where you can get an expression for the probability mass function, but it is quite a complicated expression that does not simplify. To see this, let's consider the more general case where you have independent random vectors $\mathbf{X}_1 \sim \text{Mu}(n, \mathbf{P}_1)$ and $\mathbf{X}_2 \sim \text{Mu}(n, \mathbf{P}_2)$ and you set $\mathbf{Y} \equiv \mathbf{X}_1-\mathbf{X}_2$. Then for any possible outcomes $\mathbf{y}$ we get the probability mass function:
$$\begin{align}
p_\mathbf{Y}(\mathbf{y})
&\equiv \mathbb{P}(\mathbf{Y}=\mathbf{y}) \\[12pt]
&= \sum_{\mathbf{x}} \text{Mu}(\mathbf{x}+\mathbf{y}|n, \mathbf{p}_1) \cdot \text{Mu}(\mathbf{x}|n, \mathbf{p}_2) \\[6pt]
&= \sum_{\mathbf{x} \in \mathscr{X}(\mathbf{y})} {n \choose \mathbf{x}+\mathbf{y}} {n \choose \mathbf{x}} \prod_{i=1}^k p_{1,i}^{x_i+y_i} p_{2,i}^{x_i}. \\[6pt]
\end{align}$$
where $\mathscr{X}(\mathbf{y}) \equiv \{ \mathbf{x} \in \{ 0,...,n \}^k | \min (x_i+y_i) \geqslant 0, \sum x_i = n \}$ is the set of all possible values of $\mathbf{x}$ consistent with the outcome $\mathbf{y}$. Now, this is technically a closed-form expression for the mass function, since it is a finite sum of closed form expressions. Unfortunately, it does not simplify any further. With a bit of effort you can program this function in computational software so that you automate the computation, and that is about as good as you can do here.
Of course, if you don't mind an approximation, another approach that is useful for large $n$ is to approximate the multinomial distribution by the multivariate normal distribution, so that the difference vector also have a multivariate normal distribution. You could get a reasonably approximating expression that uses the multivariate normal and then rounds the values to integers and then imposes the support constraints. That would be another way to deal with it.
