We present novel statistical tools to cross-correlate frequency cleaned thermal Sunyaev–Zel'dovich (tSZ) maps and tomographic weak lensing (wl) convergence maps. Moving beyond the lowest order cross-correlation, we introduce a hierarchy of mixed higher order statistics, the cumulants and cumulant correlators, to analyse non-Gaussianity in real space, as well as corresponding polyspectra in the harmonic domain. Using these moments, we derive analytical expressions for the joint two-point probability distribution function for smoothed tSZ (y) and convergence (κ) maps. The presence of tomographic information allows us to study the evolution of higher order mixed tSZ–wl statistics with redshift. We express the joint PDFs pκy(κ, y) in terms of individual one-point PDFs [pκ(κ), py(y)] and the relevant bias functions [bκ(κ), by(y)]. Analytical results for two different regimes are presented that correspond to the small and large angular smoothing scales. Results are also obtained for corresponding hotspots in the tSZ and convergence maps. In addition to results based on hierarchical techniques and perturbative methods, we present results of calculations based on the lognormal approximation. The analytical expressions derived here are generic and applicable to cross-correlation studies of arbitrary tracers of large-scale structure including, e.g., that of tSZ and soft X-ray background. We provide detailed comparison of our analytical results against state of the art Millennium Gas Simulations with and without non-gravitational effects such as pre-heating and cooling. Comparison of these results with gravity only simulations, shows reasonable agreement and can be used to isolate effect of non-gravitational physics from observational data.