Path integral formulation of the density matrix ρ In Feynman's Statistical Mechanics - A Set of Lectures, upon the introduction of the path integral, a series of approximations are made in order to calculate integrals. I am unsure how exactly to get to the following important approximation.
Section 3.1 Path Integral Formulation of the Density Matrix:

For low $\epsilon$, because  $\rho_\text{free}$ is a very localized
Gaussian, most of the contribution in the integral over $x''$ occurs
near $x''=x_0$ with
$$x_0=\frac{ux+(\epsilon-u)x'}{\epsilon}$$
So we can, for small $\epsilon$ , write
$$\rho(x,x';\epsilon)\approx -\int\limits_{0}^{\epsilon}V(x_0)\underbrace{\rho_\text{free}(x,x';\epsilon)}_{…?}\ \text{d}u$$

with the density matrix of the free particle
$$\rho_0(x,x',\epsilon)=\rho_{\text{free}}(x,x',\epsilon)=\sqrt{\frac{m}{2m\hbar \epsilon}}\exp\bigg(\frac{-m}{2\hbar \epsilon}(x-x')^2\bigg)$$
and the integral in question
$$\rho(x,x';\epsilon)=-\int\limits_{0}^{\epsilon}\int\limits_{-\infty}^{+\infty}\rho_0(x,x'';\epsilon-u)V(x'')\rho_0(x'',x';u)  \ \text{d}u\text{d}x''$$

I thought, that with
$$ \int\limits_{-\infty}^{+\infty}\rho_0(x,x'';\epsilon-u)\rho_0(x'',x';u)  \ \text{d}x''=\int\limits_{-\infty}^{+\infty}<x|\rho_0(\epsilon-u)\underbrace{|x''><x''|}_{\mathbb{1}}\rho_0(u)|x'>  \ \text{d}x''$$
I might get further, that still leaves me with
$$ \rho_0(\epsilon-u)\cdot\rho_0(u)$$
 A: You have more-or-less got the answer already. As you noted, the integral consists of a product of narrow Gaussian functions of the intermediate coordinate $x''$, multiplied by the potential energy function $V(x'')$. The first step is to calculate where the product of Gaussians has its maximum value: this gives the formula for $x_0$. Then $V$ is evaluated at $x''=x_0$ and taken out of the integral over $x''$.
Your concern is with the remaining integral,
which now involves just the product of
free particle density matrices.
It is the convolution identity for density matrices, which applies to any density matrix (not just the free particle one),
and holds for any value of $\epsilon$, not just small $\epsilon$. It comes from
\begin{align*}
\exp(-\epsilon \hat{H}) &= \exp(-[\epsilon-u] \hat{H}) \, \exp(-u \hat{H}) 
\\
\Rightarrow\qquad
\langle x |\, \exp(-\epsilon \hat{H}) \, | x'\rangle &=
\langle x |\, \exp(-[\epsilon-u] \hat{H}) \, \exp(-u \hat{H}) \, |x'\rangle
\\
&=
\int dx''
\langle x | \exp(-[\epsilon-u] \hat{H}) | x''\rangle \,
\langle x'' | \exp(-u \hat{H}) | x'\rangle
\end{align*}
where an integral $1=\int dx'' |x''\rangle\langle x''|$ over intermediate states has been inserted.
This is essentially what you have already!
The same identity is typically used to decompose the density matrix
into a product of $M$ density matrices all evaluated at a temperature
$M$ times higher than the true temperature  (i.e. inverse temperature $\beta\rightarrow\beta/M$),
integrated over all the intermediate points.
