We prove that the grand canonical Gibbs state of an interacting quantum Bose gas converges to the Gibbs measure of a nonlinear Schrödinger equation in the mean-field limit, where the density of the gas becomes large and the interaction strength is proportional to the inverse density. Our results hold in dimensions d \leq 3. For d>1 the Gibbs measure is supported on distributions of negative regularity and we have to renormalize the interaction. More precisely, we prove the convergence of the relative partition function and of the reduced density matrices in the L^r-norm with optimal exponent r. Moreover, we prove the convergence in the L^\infty-norm of Wick-ordered reduced density matrices, which allows us to control correlations of Wick-ordered particle densities as well as the asymptotic distribution of the particle number. Our proof is based on a functional integral representation of the grand canonical Gibbs state, in which convergence to the mean-field limit follows formally from an infinite-dimensional stationary phase argument for ill-defined non-Gaussian measures. We make this argument rigorous by introducing a white-noise-type auxiliary field, through which the functional integral is expressed in terms of propagators of heat equations driven by time-dependent periodic random potentials and can, in turn, be expressed as a gas of interacting Brownian loops and paths. When the gas is confined by an external trapping potential, we control the decay of the reduced density matrices using excursion probabilities of Brownian bridges.
