Solving for the density operator in the quantum Brownian motion master equation

I want to solve for the density operator in the quantum Brownian motion master equation,

\begin{align} \begin{aligned} \frac{d\rho_S(t)}{dt}=&-\left(\frac{i}{\hbar}\right)\Big[H_S+\frac{1}{2}M\widetilde{\Omega}^2X^2,\rho_S(t)\Big]-i\gamma[X,[P,\rho_S(t)]_+]\\ &-D[X,[X,\rho_S(t)]]-f[X,[P,\rho_S(t)]]\,, \end{aligned} \end{align} where \begin{align} \begin{aligned} H_S&=\frac{1}{2M}P^2+\frac{1}{2}M\Omega^2X^2\,,\\ \widetilde{\Omega}^2&=\frac{2}{M\hbar}\int_0^{\infty}\eta(\tau)\cos(\Omega\tau)d\tau\,,\\ \gamma&=\frac{1}{M\Omega\hbar^2}\int_0^{\infty}\eta(\tau)\sin(\Omega\tau)d\tau\,,\\ D&=\frac{1}{\hbar^2}\int_0^{\infty}\nu(\tau)\cos(\Omega\tau)d\tau\,,\\ f&=-\frac{1}{M\Omega\hbar^2}\int_0^{\infty}\nu(\tau)\sin(\Omega\tau)d\tau\,. \end{aligned} \end{align}

The only operators are $X$ and $P$ - everything else is a constant. If I want to write this master equation in the momentum basis, I know that I need to sandwich both sides of the equation with $\langle{p}|$ and $|{p}\rangle$. If I just focus upon $\Big[H_S,\rho_S(t)\Big]$, ignoring constants, I have

\begin{align} \langle{p}|\Big[P^2+X^2,\rho_S(t)\Big]|p\rangle&=\Big[p^2+(i\hbar)^2\frac{d^2}{dp^2},\rho_S(p,t)\Big]\,,\\ &=\Big[(i\hbar)^2\frac{d^2}{dp^2},\rho_S(p,t)\Big]\,. \end{align} But now I'm not sure what to do. This final commutator gives me a derivative that's acting on nothing. How do I interpret this when solving for the density operator?

Since I am going to eventually solve for the density operator in the momentum basis, is there an easier way to do so than sandwiching both sides of the master equation with the bra and ket of the momentum?

How come you took an expectation value of an operator and still got an operator instead of a c-number? In order to get a c-number equation, you will get it for the matrix elements of the density matrix, so in position space this would be $$\left\langle x | \rho | x'\right\rangle \equiv \rho(x,x')$$ I'll give you an example for the last operator and you can do the rest (I set $\hbar=1$).
$$\left\langle p | \left[ X^2 + P^2,\rho\right] | p' \right\rangle = \left\langle p | \left[ X^2 ,\rho\right] | p' \right\rangle + \left(p^2 - p'^2\right)\left\langle p |\rho|p'\right\rangle$$ $$=\left(-\frac{\partial^2}{\partial p^2}\left\langle p \right| \right) \rho \left| p' \right\rangle + \left\langle p\right|\rho\left(-\frac{\partial^2}{\partial p'^2}\left| p' \right\rangle \right) + \left(p^2 - p'^2\right)\left\langle p |\rho|p'\right\rangle$$ $$=\mathscr{F}_{pp'}\left[ \left(x'^2 - x^2\right)\rho(x,x') + \frac{\partial ^2}{\partial x^2}\rho(x,x')-\frac{\partial ^2}{\partial x'^2}\rho(x,x')\right]$$
Where $\mathscr{F}$ stands for the double fourier transform.
By the way, you taking expectation value between the same state, you are only finding the differential equation for the diagonal elements. Using different $p$ and $p'$ gives you the general matrix element.