Finite temperature correlation functions in QFT Suppose that we want to calculate this imaginary time-ordered correlation function for an interacting system (in Heisenberg picture) :
$$\langle \mathscr{T} A(\tau_A)B(\tau_B) \rangle =\frac{1}{Z} Tr\{\mathscr{T}(e^{-\beta H}  e^{H\tau_A}A(0)e^{a-H\tau_A}e^{H\tau_B}B(0)e^{-H\tau_B})\}$$
Assuming that $\tau_A > \tau_B$ we can drop the time-ordering operator $\mathscr{T}$ and using the definition of imaginary time evolution operator in interaction picture :
$$\langle \mathscr{T} A(\tau_A)B(\tau_B) \rangle = \frac{1}{Z}Tr\{ e^{-\beta H_0}S(\beta,\tau_A)A_I(\tau_A)S(\tau_A,\tau_B)B_I(\tau_B)S(\tau_B,0)\}$$
In AGD’s “Methods of quantum field theory”, it's said that we can now write the above relation as :
$$\langle \mathscr{T} A(\tau_A)B(\tau_B) \rangle = \frac{1}{Z}Tr\{e^{-\beta H_0} \mathscr{T} (A_I(\tau_A)B_I(\tau_B)S(\beta,0))\}$$
But it's only true when $\beta > \tau_A,\tau_B$. In fact, if we consider a correlation function of some operators, each one at an arbitrary time, then the finite-temperature version of the Gell-Mann–Low theorem is valid only if the time difference between operators is lower than the characteristic thermal time-scale of the system. 


*

*My question is how we can solve this problem and find the correlation function perturbatively ? 


Or maybe it's not a problem and it's natural that we can't find any correlation between two (or more) quantities in a system in thermal equilibrium when their time separation is so long that the thermal fluctuations kill any correlation between them. 
And even if the latter statement is true, it doesn't say that we can't find the correlation, it just says we must use another method, for example solving the Hamiltonian exactly. Now suppose that we could do this. 


*

*What physical behaviors we expect from it to have (maybe some kind of exponential decrease as the time difference increases with an exponential factor of $e^{-\frac{|\tau_A-\tau_B|}{\beta}}$? Or maybe I'm misled completely because the imaginary time have nothing to do with the real time. 

 A: Although there is already an accepted answer, I'll provide a different point of view here.
It is important to keep in mind that the imaginary time is periodic, which means fields satisfy periodic(or anti-periodic) boundary conditions (the reason is because we use imaginary time path integral to represent the partition function, which is $\mathrm{Tr}\, e^{-\beta H}$). Therefore, the correlation functions also obey the same kind of periodic conditions under the shift $\tau\rightarrow \tau+\beta$, which is why the Matsubara frequencies are $\frac{2\pi n}{\beta}$ for $n\in\mathbb{Z}$ for bosonic fields, and $\frac{(2n+1)\pi}{\beta}$ for fermionic fields. So if $|\tau_A-\tau_B|$ is larger than $\beta$, all you need to do is to fold it back by adding or subtracting several $\beta$'s, just as one would do to any periodic functions. 
EDIT: The periodicity of Green's function can be explicitly derived, without using the path integral representation. After all, the path integral representation is designed to be consistent with the operator formalism.
Let's assume $A$ and $B$ are bosonic, so I do not need to keep track of the fermionic sign. For concreteness, consider $-\beta<\tau<0$ and 
$$
G(\tau)=\langle \mathcal{T} A(\tau)B(0)\rangle=\mathrm{Tr}[e^{-\beta H}B(0)e^{\tau H}A(0)e^{-\tau H}]=\mathrm{Tr}[ e^{\tau H}A(0)e^{-\tau H}e^{-\beta H}B(0)]\\
=\mathrm{Tr}[e^{-\beta H}e^{(\tau+\beta) H}A(0)e^{-(\tau+\beta) H}B(0)]=\langle \mathcal{T} A(\tau+\beta) B(0)\rangle
=G(\tau+\beta)
$$
Here we use the cyclic property of the trace.
A: OP's question is about the long-time behavior of imaginary-time correlation functions in general. But in fact, the correlation function is ill-defined if the time difference $|\tau_A - \tau_B|$ is larger than the inverse temperature $\beta$.
To see this, suppose that the Hamiltonian $H$ has eigenstates $\{ | n \rangle \}$ and associated eigenvalues $\{ E_{n} \}$, and assume $\tau_A > \tau_B$ . Then,
\begin{equation}
\begin{split}
&\mathrm{Tr}\big\{\mathscr{T}[e^{-\beta H}e^{\tau_{A}H} \, A(0) \, e^{-\tau_{A}H}e^{\tau_{B}H} \, B(0) \, e^{-\tau_{B}H}]\big\}\\
&= \sum_{n,n^\prime} \langle n| \, e^{-\beta H}e^{\tau_{A}H} \, A(0) \, |n^\prime\rangle\langle n^\prime|    e^{-\tau_{A}H}e^{\tau_{B}H} \, B(0) \, e^{-\tau_{B}H}| n\rangle\\
&=\sum_{n,n^\prime}e^{-(\beta-\tau_A+\tau_B)E_n} e^{-(\tau_A-\tau_B)E_{n^\prime}}\langle n| \, A(0) \, |n^\prime\rangle\langle n^\prime| \, B(0) \, |n\rangle.
\end{split}
\end{equation}
The energy eigenvalues $\{E_{n}\}$ of a physical system are bounded below and unbounded above. Hence, for both $e^{-(\beta-\tau_A+\tau_B)E_n}$ and $e^{-(\tau_A-\tau_B)E_{n^\prime}}$ to well behave for all $n$ and $n^\prime$, we must have $0\le\tau_A - \tau_B\le\beta$. Similarly, if $\tau_A < \tau_B$, we have $0\le\tau_B - \tau_A\le\beta$.
Combining the two cases $\tau_A > \tau_B$ and $\tau_A < \tau_B$, it follows that the imaginary-time correlation function is well-defined only if
\begin{equation}
|\tau_A - \tau_B| <\beta.
\end{equation}
