We present novel statistical tools to cross-correlate frequency cleaned thermal Sunyaev-Zeldovich (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 analyze 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 (2PDF) for smoothed tSZ (y_s) and convergence (kappa_s) maps. The presence of tomographic information allows us to study the evolution of higher order {em mixed} tSZ-weak lensing statistics with redshift. We express the joint PDFs p_{kappa y}(kappa_s,y_s) in terms of individual one-point PDFs (p_{kappa}(kappa_s), p_y(y_s)) and the relevant bias functions (b_{kappa}(kappa_s), b_y(y_s)). Analytical results for two different regimes are presented that correspond to the small and large angular smoothing scales. Results are also obtained for corresponding {em hot spots} 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.