Adler, Shiota and van Moerbeke observed that a tau function of the Pfaff lattice is a square root of a tau function of the Toda lattice hierarchy of Ueno and Takasaki. In this paper, we give a representation theoretical explanation for this phenomenon. We consider 2-BKP and two-component 2-KP tau functions. We shall show that a square of a BKP tau function is equal to a certain two-component KP tau function and a square of a 2-BKP tau function is equal to a certain two-component 2-KP tau function.

We introduce a family of models, which we name matrix models associated with children’s drawings—the so-called dessin d’enfant. Dessins d’enfant are graphs of a special kind drawn on a closed connected orientable surface (in the sky). The vertices of such a graph are small disks that we call stars. We attach random matrices to the edges of the graph and get multimatrix models. Additionally, to the stars we attach source matrices. They play the role of free parameters or model coupling constants. The answers for our integrals are expressed through quantities that we call the “spectrum of stars”. The answers may also include some combinatorial numbers, such as Hurwitz numbers or characters from group representation theory.

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.

We study multimatrix models, which may be viewed as integrals of products of tau functions depending on the eigenvalues of products of random matrices. We consider tau functions of the two-component Kadomtsev–Petviashvili (KP) hierarchy (semi-infinite relativistic Toda lattice) and of the B-type KP (BKP) hierarchy introduced by Kac and van de Leur. Such integrals are sometimes tau functions themselves. We consider models that generate Hurwitz numbers H^E,F, where E is the Euler characteristic of the base surface and F is the number of branch points. We show that in the case where the integrands contain the product of n > 2 matrices, the integral generates Hurwitz numbers with E ≤ 2 and F ≤ n+2. Both the numbers E and F depend both on n and on the order of the factors in the matrix product. The Euler characteristic E can be either an even or an odd number, i.e., it can match both orientable and nonorientable (Klein) base surfaces depending on the presence of the tau function of the BKP hierarchy in the integrand. We study two cases, the products of complex and the products of unitary matrices.