I think there would be no dificulties here. For example, in the Aharonov-Bohm effect studied in a Young Interferometer with a solenoid placed between the slits, we describe the input state (before solenoid) as
$$
|\Psi\rangle = {1\over \sqrt 2} \left(|\Psi_L\rangle + |\Psi_R\rangle\right)
$$
and after the solenoid
$$
|\Psi\rangle = {1\over \sqrt 2} \left(|\Psi_L\rangle + e^{ie\Phi_B/\hbar c}|\Psi_R\rangle\right)
$$
Now, to change this description for density matrices, we just need to consider it through definition
$$
\rho(\Phi_B) = |\Psi\rangle\!\langle \Psi | = \frac 12 \left( |\Psi_L\rangle\!\langle \Psi_L | + e^{ie\Phi_B/\hbar c}|\Psi_L\rangle\!\langle \Psi_R | +e^{-ie\Phi_B/\hbar c} |\Psi_R\rangle\!\langle \Psi_L | + |\Psi_R\rangle\!\langle \Psi_R |\right)
$$
The system described here is equivalent to a two-level system in the slits basis. We will have
$$
\rho(\Phi_B) = \frac 12
\begin{pmatrix}
1 & e^{ie\Phi_B/\hbar c}\\
e^{-ie\Phi_B/\hbar c} & 1
\end{pmatrix}
$$
All interference effects could be obtained by this density matrix
$$
P(x) = \text{tr}(\rho(\Phi_B) |x\rangle\!\langle x|) = \frac 12 \left(|\Psi_L(x)|^2 + |\Psi_R(x)|^2 + 2\text{Re}\left(e^{ie\Phi_B/\hbar c}\Psi_L(x)^*\Psi_R(x)\right) \right)
$$
The interesting feature in this representation is that we can obtain the state of the system in the case where the magnetic flux can't be measured precisely. In this case we have a distribution $g(\Phi_B)$ and the state is obtained through the mean of states by all possible fluxes
$$
\rho = \int_0^\infty g(\Phi_B) \rho(\Phi_B) d\Phi_B
$$
This is necessary, for example, to explain the destruction of interference fringes when the uncertainty about the magnetic flux is big. If we model the random magnetic flux as a gaussian random variable, for example, $g(\phi_B) = \sqrt{1\over 2\pi \Delta \Phi_B^2} e^{-\Phi_B^2/2\Delta\Phi_B^2}$, we get after integration
$$
\rho = \frac 12
\begin{pmatrix}
1 & e^{-e^2 \Delta \Phi_B^2\over 2\hbar^2 c^2}\\
e^{-e^2 \Delta \Phi_B^2\over 2\hbar^2 c^2}& 1
\end{pmatrix}
\stackrel{\Delta \Phi_B \rightarrow \infty}{\rightarrow} \frac 12 \begin{pmatrix}
1 &0\\
0& 1
\end{pmatrix}
$$
