Why is the Matsubara Green function $\mathscr{G}(i\omega_n)$ equal to the retarded Green function (also the linear response susceptibility) $\chi(\omega+i\epsilon)$ under the substitution $i\omega_n \mapsto \omega+i\epsilon$.
I understand that you can compute the spectral representation and find that the 2 are equal. However, this seems rather unsatisfying, since there should be some motivation (hopefully) of defining the Masturbara Green function as $$ \mathscr{G}(\tau)=-\langle TA(\tau)A\rangle $$ which ultimately results in the fact that it is equal to $$ \chi(t)=-i\theta(t)\langle [A(t),A]\rangle $$ after Fourier transform.
Of course, it is also possible that the thought process is actually done in reverse, i.e., we first calculate $\chi(z)$ where $z \in \mathbb{C}^+$ is in the upper half-plane, and think: you know what? What if we tried computing $\chi(i\omega_n)$ and see what happens when we apply the inverse Fourier transform to get $\chi(\tau)$. What do you know, it happens to be a desired form of time-ordering correlation function in imaginary time.