In this paper, we propose two new representation formulas for the conditional marginal probability density of the multi-factor Heston model. The two formulas express the marginal density as a convolution with suitable Gaussian kernels whose variances are related to the conditional moments of price returns. Via asymptotic expansion of the non-Gaussian function in the convolutions, we derive explicit formulas for European-style option prices and implied volatility. The European option prices can be expressed as Black–Scholes style terms plus corrections at higher orders in the volatilities of volatilities, given by the Black–Scholes Greeks. The explicit formula for the implied volatility clearly identifies the effect of the higher moments of the price on the implied volatility surface. Further, we derive the relationship between the VIX index and the variances of the two Gaussian kernels. As a byproduct, we provide an explanation of the bias between the VIX and the volatility of total returns, which offers support to recently proposed methods to compute the variance risk premium. Via a series of numerical exercises, we analyse the accuracy of our pricing formula under different parameter settings for the one- and two-factor models applied to index options on the S&P500 and show that our approximation substantially reduces the computational time compared to that of alternative closed-form solution methods. In addition, we propose a simple approach to calibrate the parameters of the multi-factor Heston model based on the VIX index, and we apply the approach to the double and triple Heston models.