I'm reading this paper. But I haven't read anything on how to calculate the density operator in a QFT or how to calculate its trace. Now I can't follow this part of the paper. Can anyone clarify?

Well to answer your question, it seems one first needs to learn CFT, I have the book of Di Franchesko but I still didn't read it, so I guess I cannot help here, what is your knowledge of CFT?

I know some CFT but actually its not necessary for this question. My question is about calculating the density matrix of any QFT. Maybe @samalkhaiat can help.

Density functional is different from density matrix(operator). In QM, density operator is calculated by the formula ## \rho=\sum_n p_n |\phi_n\rangle \langle \phi_n |## where ## \{|\phi_n\rangle\} ## is a set of not necessarily orthonormal states. I want to know how to calculate the corresponding quantity in QFT.

Tadashi Takayanagi is good young physicist, I was impressed by his lectures at CERN Winter School on Supergravity 2012. He presented the subject very well, I suggest (if you have not done so) you look at them. As you might know, to explain this properly, one needs to draw pictures. Unfortunately I cannot do this, but this does not prevent us from making formal analogy with QM. In statistical physics, the density matrix in a thermal state at (inverse) temperature [itex]\beta[/itex] is given by [tex]\rho (x,y) \propto \langle x | e^{-\beta H} | y \rangle .[/tex] The proportionality constant which ensures [tex]\mbox{Tr}( \rho ) = \int dx \ \rho (x,x) = 1 ,[/tex] is given by [tex]\frac{1}{Z} = \frac{1}{\mbox{Tr}( e^{-\beta H})} .[/tex] One can (if one wants to) express [itex]\rho[/itex] as path integral over the variables [itex]x(\tau)[/itex] which take “boundary” values [itex]y[/itex] at [itex]\tau = 0[/itex], and [itex]x[/itex] at [itex]\tau = \beta[/itex]. So, you should have no problem with path integral expression [tex]\rho (x , y) = \frac{\langle x | e^{-\beta H} | y \rangle }{Z} = (1/Z) \int [dx(\tau)] \ \delta (x(0) - y) \delta (x(\beta) - x) \ e^{- \int^{\beta}_{0} \ L d\tau} . \ \ \ (1)[/tex] Taking the trace, i.e., setting [itex]x = y[/itex] and integrating, has, therefore, the effect of “gluing the edges” along [itex]\tau = 0[/itex] and [itex]\tau = \beta[/itex] and forming a "cylinder" of circumference [itex]\beta[/itex].
Okay, let us formally generalize the above to a bosonic field theory in (1+1)-dimension: the operator [itex]\hat{x}[/itex] will be “replaced” by a complete set of local observables [itex]\{ \hat{\varphi}(x) \}[/itex], and [itex]\hat{x}|x\rangle = x |x\rangle[/itex] with the eigenvalue equations [tex]\{ \hat{\varphi}(x) \} \left( \otimes_{x}| \{ \varphi (x) \} \rangle \right) = \{ \varphi (x) \} \left( \otimes_{x}| \{ \varphi (x) \} \rangle \right) .[/tex] So, instead of (1), you will have
[tex]
\begin{align*}
\rho \left( \varphi^{''}(x^{''}) , \varphi^{'}(x^{'}) \right) &= \frac{ \langle \varphi^{''}(x^{''}) | e^{-\beta H} | \varphi^{'}(x^{'}) \rangle}{Z} \\
&= \frac{1}{Z} \int [d\varphi (x , \tau)] \prod_{x}\left[ \delta \left( \varphi (x,0) - \varphi^{'}(x^{'}) \right) \delta \left( \varphi (x ,\beta) - \varphi^{''}(x^{''}) \right) \right] e^{-\int_{0}^{\beta} L_{E} d\tau} \ \ \ \ (2)
\end{align*}
[/tex]
Again, the delta functions are there to enforce the boundary values [itex]\varphi^{'}(x)[/itex] at [itex]\tau = 0[/itex], and [itex]\varphi^{''}(x)[/itex] at [itex]\tau = \beta[/itex]. And the trace is found by setting [itex]\varphi^{''}(x) = \varphi^{'}(x)[/itex] and integrating over these variables except, in this case, you really have edges to glue along [itex]\tau = 0[/itex] and [itex]\tau = \beta[/itex].
So, if you want to calculate [itex]\mbox{Tr}(\rho^{n})[/itex], you will need to take the products of n path integral of the form (2) which will be equivalent to path integral over n-sheeted Reimann surface.

For an intro to relativistic quantum many-body theory (QFT), I can recommend three textbooks

J. I. Kapusta and C. Gale, Finite-Temperature Field Theory; Principles and Applications, Cambridge University Press, 2 ed., 2006.
M. LeBellac, Thermal Field Theory, Cambridge University Press, Cambridge, New York, Melbourne, 1996.

It seems to me you're talking about something else!

1) The calculations in sec.3.1 of the paper are done for a zero temperature CFT, i.e. ##-\infty<t_E<\infty##.

2) The density matrix for the field in all points is ## [\rho]_{\phi_0\phi_0'}=\Psi(\phi_0)\overline \Psi(\phi'_0) ##. Since ## \displaystyle \Psi(\phi_0)\propto \int_{t_E=-\infty}^{\phi(t_E=0,x)=\phi_0(x)} D \phi e^{-S[\phi]}## and ## \overline\Psi(\phi'_0)\propto \displaystyle \int^{t_E=\infty}_{\phi(t_E=0,x)=\phi'_0(x)} D \phi e^{-S[\phi]} ##, the density matrix should be ## \displaystyle [\rho]_{\phi_0\phi'_0}\propto \int_{t_E=-\infty}^{\phi(t_E=0,x)=\phi_0(x)} \int^{t_E=\infty}_{\phi'(t_E=0,x)=\phi'_0(x)} D \phi D \phi' e^{-S[\phi]} e^{-S[\phi']}##.

3) The path integral with those dirac deltas, is actually the reduced density matrix which is calculated by integrating out all the degrees of freedom outside of region A and the density matrix for all points of space should have no dirac delta! This means that in the density matrix above, we should put ## \phi(x)=\phi'(x) ## and ## \phi_0(x)=\phi'_0(x) ## except for ## x\in A ##. Taking a full trace gives ## \displaystyle Tr \rho=\int_{t_E=-\infty}^{\phi(t_E=0,x)=\phi_0(x)} \int^{t_E=\infty}_{\phi'(t_E=0,x)=\phi'_0(x)} D \phi D \phi' e^{-S[\phi]} e^{-S[\phi']} \prod_{x} \delta( \phi(t_E,x)-\phi'(t_E,x) )=\int_{t_E=-\infty}^{t_E=\infty}D \phi e^{-2S[\phi]}##. But because we don't want the full trace and just want to integrate out degrees of freedom outside of A and have a reduced density matrix for A, we should avoid integrating over all fields for ## x\in A ##. We also need two boundary conditions at zero(## \phi_{\pm}(x) ##). So we put those dirac deltas in the formula for the full trace and get:
## \displaystyle [\rho_A]_{\phi_+\phi_-}\propto \int_{t_E=-\infty}^{t_E=\infty}D \phi e^{-2S[\phi]} \prod_{x\in A} \delta(\phi(+0,x)-\phi_+(x))\delta(\phi(-0,x)-\phi_-(x)) ##.

But as you can see, instead of ## e^{-S[\phi]} ##, I have ## e^{-2S[\phi]} ##. Also I'm not sure about all these because its just my guts filling the blanks!!!
I also don't quite understand how ## \displaystyle \int_{t_E=-\infty}^{\phi(t_E=0,x)=\phi_0(x)} \int^{t_E=\infty}_{\phi'(t_E=0,x)=\phi'_0(x)} \to \int_{t_E=-\infty}^{t_E=\infty} ## because the second integral should simply give 1 and shouldn't extend the integration interval of the other integral. But it seems that should be done to get the desired result!

You integrate "over all paths", where in QFT the paths are in field-configuration space. Your first expression describes the integration over all paths first involving only the part where ##t<0## and at ##t=0## the field is fixed as ##\phi_0(t=0,\vec{x})=\phi_0(\vec{x})## and then a part where ##t>0## and ##\phi_0(t=0,\vec{x})=\phi_0'(\vec{x})##, leading to a functional of ##\phi_0## and ##\phi_0'##.

Now if you want the trace (i.e., the partition sum) you set ##\phi_0=\phi_0'##, i.e., you integrate first for ##t<0## and ##t>0## with fixed field configuration ##\phi_0## for both parts. Finally you integrate also over ##\phi_0##. So altogether you integrate in this final step over all field configurations without restrictions, and that's the meaning of the trace.

Thanks guys. I now understand it clearly.
But I just realized there is another problem. In the paper in the OP, the ground state wave functional is given by ## \displaystyle \Psi(\phi_0(x))=\int_{t_E=-\infty}^{\phi(t_E=0,x)=\phi_0(x)} D\phi e^{-S[\phi]} ##. As you can see the argument of the wave functional appears as the ## t_E=0 ## boundary condition of the path integral. But the calculations mentioned in this thread, give a wave functional whose argument appears as the ## \pm \infty ## boundary condition of the path integral. Where does this discrepancy come from?

EDIT:
Is ## t=it_E ## only a convention for Wick rotation or is it the only right way? Because if I'm allowed to do the Wick rotation with ## t=-it_E ##, the above problem would be solved!