# Reduced density matrix time evolution

So, I'm trying to study the interaction between a harmonic oscillator-like particle (with variables $$X$$ and $$P$$) and a heat bath, also modeled by a compound of $$N$$ harmonic oscillators, with variables $$x_n$$ and $$p_n$$. I have derived the following hamiltonian:

$$H=\frac{P^2}{2M}+\frac{1}{2}M\Omega_0X^2+\sum_{n=1}^N \left( \frac{p_n^2}{2m}+\frac{1}{2}m\omega x_n^2 \right)+gX\sum_n x_n =H_S+H_B+H_{int}\ \ ,$$ where $$g$$ acts as a coupling constant (the same for all pairs of oscillators). In this fashion, I want to find a master equation for the reduced density operator of the central oscillator, $$\rho_S$$. I have thought the following: the combined system of central oscillator plus bath is a closed system, therefore the time evolution of the total density operator is unitary: $$\rho_T (t)=\mathcal{U}(t)\rho_T(0)\mathcal{U}^{\dagger}(t) \ \ ,$$

with $$\mathcal{U}$$ being the time evolution operator $$\mathcal{U}(t)=\exp\left( -\frac{i}{\hbar}Ht\right)$$. The definition of the reduced density matrix is clear, as we only have to trace over the bath coordinates $$\left\{x_n,p_n \right\}$$ which in the continuum limit transforms into integrals:

$$Tr_B \longrightarrow \int \prod_n dp_ndx_n$$

The last part of this approach is assigning a density operator at time zero. Assuming uncorrelated variables, the total density operator can be expressed as the outer product of the reduced density operators both from the system (central oscillator) and the bath $$\rho_T (0)=\rho_S (0) \otimes \rho_B(0) \ \ ,$$ where the latter can be assumed to be a Gibbs state in thermodynamic equilibrium at temperature $$T$$, namely $$\rho_B(0)=\frac{\exp\left( -\beta H_B\right)}{Z_B} \ \ ,$$ where $$Z_B$$ is the initial partition function of the bath and $$H_B$$ is the part of the hamiltonian corresponding to the bath (without the interaction part) . Now, my reasoning has been the following, I just need to derivate with respect to time the expression of $$\rho_S$$ as a partial trace of $$\rho_T$$: $$\frac{d}{dt}\rho_S=\frac{d}{dt}\left[ Tr_B \left( \mathcal{U}(t)\rho_T(0)\mathcal{U}^{\dagger}(t) \right) \right]=\frac{d}{dt}\int \prod_n dp_n dx_n \mathcal{U}(t)\left( \rho_S(0) \otimes \rho_B(0)\right)\mathcal{U}^{\dagger}(t)$$

But I'm kinda stuck here, as this gets algebraically difficult for me, and also taking into account that I don't have any information on the initial system's density operator, $$\rho_S(0)$$. I have tried a couple of things, which include:

1. Assuming that $$\rho_S (0)$$ can also be described as a thermodynamic equilibrium distribution; this turned out to be useless, as plugging the same $$\beta$$ just throws off the trivial result that the system doesn't evolve (of course), and plugging a different $$\beta$$ just doesn't help.

2. Expanding the time evolution operator as an infinite series, but I got lost in the algebraic mess that emerged.

Any help or tips are appreciated.

Note: Reading through a vast lot of literature I've read about the Caldeira-Leggett and Rubin models for heat baths, but I just see everything done as "trivial" and want to develop the reasoning as a whole.

Edit: I have thought about using von-Neumann's equation, since the trace over the bath coordinates can be swapped with the time derivatives to yield the time derivative of $$\rho_S$$ and since the hamiltonian and $$\rho_T$$ don't generally commute, maybe something can be worked out. Is this right?