You may notice that they are re-deriving the Dyson equation (or whatever it is called): G(b,a)=G0(b,a)+∫dcG(b,c)V(c)G0(c,a), where a,b,c denote complete sets of indices, say in this case (r,τ), G(b,a)=−⟨Tψ(b)ψ†(a)⟩, and V is an external potential. This equation is one that can be used to derive Feynman's rules and it can be easily derived just by equation of motion of the Hamiltonian operator. (If you are interested, check Chapter 11, Many-body quantum theory in condensed matter physics. ) Anyway, let me try to answer your questions:
(1) Yes, Gn(r,r′)=⟨⟨c†n(r)cn(r′)⟩⟩imp is a scalar for fixed r and r′, but if we treat r and r′ as indices, Gn(r,r′) can also be regarded an infinite-dimensional matrix (an operator if you prefer, but in a different sense, depending on which space it acts.) For clarification, let me show all indices explicitly in the relation: ˆG(r,r′)=1iωnδ(r−r′)−ˆH0(r,r′)−ˆHimp(r,r′)=1ˆG−10(r,r″)[1−ˆG0(r,r″)ˆHimp(r″,r′)]. Summation over repeated indices is assumed. The expansion should also carry these indices with caution. After that, you can get the expression for ˆG(r,r′).
(2) Notice that there are always two common ways to evaluate Green's functions (correlation functions, propagators...), namely, using operators and using functional integration (path integral). If you do not know the difference, compare Eq(4.31) with Eq(9.18) in Peskin and Schroeder. In path integral, these fields are scalars.
(3) Yes. You are right. It seems to me that the ''scattering lines'' should connect to a single point. Hope this can help :~)