Not a full answer, but some general considerations that should lead to one when followed through.
Denote with $\chi_\rho(\alpha,1)\equiv \operatorname{tr}(D_1(\alpha)\rho)$ the normally-ordered characteristic function of a state $\rho$, with $D_1(\alpha)\equiv e^{\alpha a^\dagger}e^{-\bar\alpha a}$. The $P$ function can be written as the Fourier transform of this characteristic function:
$$P_\rho(\nu) = \int\frac{\mathrm d^2\alpha}{\pi} e^{\nu\bar\alpha-\bar\nu\alpha} \chi_\rho(\alpha,1) \equiv \mathcal F[\chi_\rho(\bullet,1)](\nu),$$
where I used $\mathcal F$ as shorthand for the Fourier transform. Note that with the convention above this operator satisfies $\mathcal F^2=\mathcal F^{-1}=\operatorname{Id}$.
The above generalises to other $s$-ordered quasiprobability distributions, such as Wigner and $Q$, by writing
$$W_\rho(\nu,s) = \mathcal F[\chi_\rho(\bullet,s)](\nu),$$
with $P_\rho(\nu)\equiv W_\rho(\nu,1)$, and $\chi_\rho(\alpha,s)\equiv e^{\frac s2|\alpha|^2}\operatorname{tr}(e^{\alpha a^\dagger-\bar\alpha a} \rho)$.
The standard convolution properties of the Fourier transform then tell us that
$$W_\rho(\bullet,s)\star W_\sigma(\bullet,s) = \mathcal F[\chi_\rho(\bullet,s)\chi_\sigma(\bullet,s)].$$
In words, the convolution of two quasiprobabilities is the Fourier transform of the pointwise product of the characteristic functions (in any fixed ordering $s$).
In particular,
$$P_\rho\star P_\sigma = \mathcal F[\chi_\rho(\bullet,1)\chi_\sigma(\bullet,1)].$$
So the question is equivalent to: what is the state (if any exists) whose characteristic function is the pointwise product of the characteristic functions of $\rho$ and $\sigma$?
Example of coherent states
In the case of coherent states, this is easy enough to compute: for coherent states $|\beta\rangle$ and $|\gamma\rangle$ we have
$$\chi_\beta(\alpha,1) = e^{\alpha\bar\beta-\bar\alpha\beta},
\qquad
\chi_\gamma(\alpha,1) = e^{\alpha\bar\gamma-\bar\alpha\gamma},$$
and thus
$$\chi_\beta(\alpha,1)\chi_\gamma(\alpha,1) = \chi_{|\beta+\gamma\rangle}(\alpha,1).$$
In other words, the convolution of the $P$ functions of two coherent states $|\beta\rangle$ and $|\gamma\rangle$ is the $P$ function of the coherent state $|\beta+\gamma\rangle$.
General case
This isn't complete, but observe that given any channel $\Phi$ we can write
$$\chi_{\Phi(\rho,\sigma)}(\alpha,s) = e^{\frac s2|\alpha|^2}\operatorname{tr}[D(\alpha) \Phi(\rho\otimes\sigma)]
= e^{\frac s2|\alpha|^2}\operatorname{tr}[\Phi^\dagger(D(\alpha)) (\rho\otimes\sigma)],$$
where $\Phi^\dagger$ is the adjoint channel of $\Phi$.
This tells you that if you know how to "evolve back" displacement operators, you know the characteristic function of the evolved state.
And if you can find a channel $\Phi$ such that $\Phi^\dagger(D(\alpha))\simeq D(\alpha_1)\otimes D(\alpha_2)$, then the resulting characteristic function of $\Phi(\rho\otimes\sigma)$ is indeed a pointwise product of the individual characteristic functions, and thus the $P$-function of $\Phi(\rho\otimes\sigma)$ is the convolution of the individual $P$ functions.
Evolution through a beamsplitter
Consider two generic states $\rho\otimes\sigma$ evolving through a balanced beamsplitter. This means in particular that coherent states evolve as
$$|\alpha\rangle\otimes|\beta\rangle\to \left\lvert\frac{\alpha+\beta}{\sqrt2}\right\rangle\otimes \left\lvert\frac{\alpha-\beta}{\sqrt2}\right\rangle.$$
Using this we can derive a formal expression for the result of evolving $\rho\otimes\sigma$ through a beamsplitter, via their $P$ representations:
$$\rho\otimes\sigma=
\int\mathrm d^2\alpha\mathrm d^2\beta
P_\rho(\alpha)P_\sigma(\beta) (\mathbb{P}_\alpha\otimes\mathbb{P}_\beta)
\\
\to
\int\mathrm d^2\alpha\mathrm d^2\beta
P_\rho(\alpha)P_\sigma(\beta) (\mathbb{P}_{(\alpha+\beta)/\sqrt2}\otimes\mathbb{P}_{(\alpha-\beta)/\sqrt2}) ,$$
using the shorthand notation $\mathbb{P}_\alpha\equiv |\alpha\rangle\!\langle\alpha|$.
If we only look at the marginal state on the first output mode, tracing out the second one, we get
$$\rho' =
\int\mathrm d^2\alpha\mathrm d^2\beta
P_\rho(\alpha)P_\sigma(\beta) \mathbb{P}_{(\alpha+\beta)/\sqrt2}
\\
= \int\mathrm d^2\gamma
\left(
\frac{\mathrm d^2\eta}{4}
P_\rho(\tfrac{\gamma+\eta}{\sqrt2})P_\sigma(\tfrac{\gamma-\eta}{\sqrt2})
\right) \mathbb{P}_\gamma
\\
=\int\mathrm d^2\gamma
\left(
\mathrm d^2\theta
P_\rho(\sqrt2\gamma-\theta)P_\sigma(\theta)
\right) \mathbb{P}_\gamma.$$
In words, we found that the state obtained evolving $\rho\otimes\sigma$ through a symmetric beamsplitter and looking only at the first output mode has $P$ representation equal to the convolution of the $P$ representations of $\rho$ and $\sigma$, modulo a rescaling of the input parameter.
