Let $ξ$ be a Lévy process and $I_ξ(t):=\int_{0}^te^{-ξ_s}\mathrm{d} s$, $t\geq 0,$ be the exponential functional of Lévy processes on deterministic horizon. Given that $\lim_{t\to \infty}ξ_t=-\infty$ we evaluate for general functions $F$ an upper bound on the rate of decay of $\mathbb{E}\left(F(I_ξ(t))\right)$ based on an explicit integral criterion. When $\mathbb{E}\left(ξ_1\right)\in\left(-\infty,0\right)$ and $\mathbb{P}\left(ξ_1>t\right)$ is regularly varying of index $α>1$ at infinity, we show that the law of $I_ξ(t)$, suitably normed and rescaled, converges weakly to a probability measure stemming from a new generalisation of the product factorisation of classical exponential functionals. These results substantially improve upon the existing literature and are obtained via a novel combination between Mellin inversion of the Laplace transform of $\mathbb{E}\left(I^{-a}_ξ(t)\mathbf{1}_{\left\{I_ξ(t)\leq x\right\}}\right)$, $a\in (0,1)$, $x\in(0,\infty],$ and Tauberian theory augmented for integer-valued $α$ by a suitable application of the one-large jump principle in the context of the de Haan theory. The methodology rests upon the representation of the aforementioned Mellin transform in terms of the recently introduced bivariate Bernstein-gamma functions for which we develop the following new results of independent interest (for general $ξ$): we link these functions to the $q$-potentials of $ξ$; we show that their derivatives at zero are finite upon the finiteness of the aforementioned integral criterion; we offer neat estimates of those derivatives along complex lines. These results are useful in various applications of the exponential functionals themselves and in different contexts where properties of bivariate Bernstein-gamma functions are needed. $ξ$ need not be non-lattice.