Hurwitz numbers from Feynman diagrams

To obtain a generating function of the most general form for Hurwitz numbers with an arbitrary base surface and arbitrary ramification profiles, we consider a matrix model constructed according to a graph on an oriented connected surface $\Sigma$ with no boundary. The vertices of this graph, called stars, are small discs, and the graph itself is a clean dessin d'enfants. We insert source matrices in boundary segments of each disc. Their product determines the monodromy matrix for a given star, whose spectrum is called the star spectrum. The surface $\Sigma$ consists of glued maps, and each map corresponds to the product of random matrices and source matrices. Wick pairing corresponds to gluing the set of maps into the surface, and an additional insertion of a special tau function in the integration measure corresponds to gluing in Möbius strips. We calculate the matrix integral as a Feynman power series in which the star spectral data play the role of coupling constants, and the coefficients of this power series are just Hurwitz numbers. They determine the number of coverings of $\Sigma$ (or its extensions to a Klein surface obtained by inserting Möbius strips) for any given set of ramification profiles at the vertices of the graph. We focus on a combinatorial description of the matrix integral. The Hurwitz number is equal to the number of Feynman diagrams of a certain type divided by the order of the automorphism group of the graph.