  Modular graph functions associate to a graph an $SL(2,Z)$-invariant function
on the upper half plane. We obtain the Fourier series of modular graph
functions of arbitrary weight $w$ and two-loop order. The motivation for this
work is to develop a deeper understanding of the origin of the algebraic
identities between modular graph functions which have been discovered recently,
and of the relation between the existence of these identities and the
occurrence of cusp forms. We show that the constant Fourier mode, as a function
of the modulus $\tau$, consists of a Laurent polynomial in $y = \pi \, {\rm Im}
(\tau)$ of degree $(w,1-w)$, plus a contribution which decays exponentially as
$y \to \infty$. The Laurent polynomial is a linear combination with rational
coefficients of the top term $y^w$, and lower order terms $\zeta (2k+1)
y^{w-2k-1}$ for $1\leq k \leq w-1$, as well as terms $\zeta (2w-2\ell-3) \zeta
(2\ell+1)y^{2-w}$ for $1 \leq \ell \leq w-3$. The exponential contribution is a
linear combination of exponentials of $y$ and incomplete $\Gamma$-functions
whose coefficients are Laurent polynomials in $y$ with rational coefficients.
