We study quartic matrix models with partition function exp(trace). The integral is over the space of Hermitean -matrices, the external matrix E encodes the dynamics, is a scalar coupling constant and the matrix J is used to generate correlation functions. For E not a multiple of the identity matrix, we prove a universal algebraic recursion formula which gives all higher correlation functions in terms of the 2-point function and the distinct eigenvalues of E. The 2-point function itself satisfies a closed non-linear equation which must be solved case by case for given E. These results imply that if the 2-point function of a quartic matrix model is renormalisable by mass and wavefunction renormalisation, then the entire model is renormalisable and has vanishing beta-function. As the main application we prove that Euclidean -quantum field theory on four-dimensional Moyal space with harmonic propagation, taken at its self-duality point and in the infinite volume limit, is exactly solvable and non-trivial. This model is a quartic matrix model, where E has for the same spectrum as the Laplace operator in four dimensions. Using the theory of singular integral equations of Carleman type we compute (for and after renormalisation of ) the free energy density (1/volume) log exactly in terms of the solution of a non-linear integral equation. Existence of a solution is proved via the Schauder fixed point theorem. The derivation of the non-linear integral equation relies on an assumption which in subsequent work is verified for coupling constants .
