  We study the equation $M_\Psi(z+1)=\frac{-z}{\Psi(-z)}M_\Psi(z), M_\Psi(1)=1$
defined on a subset of the imaginary line and where $\Psi$ is a negative
definite functions. Using the Wiener-Hopf method we solve this equation in a
two terms product which consists of functions that extend the classical gamma
function. These functions are in a bijection with Bernstein functions and for
this reason we call them Bernstein-gamma functions. Via a couple of computable
parameters we characterize of these functions as meromorphic functions on a
complex strip. We also establish explicit and universal Stirling type
asymptotic in terms of the constituting Bernstein function. The decay of
$|M_{\Psi}(z)|$ along imaginary lines is computed. Important quantities for
theoretical and applied studies are rendered accessible.
  As an application we investigate the exponential functionals of Levy
Processes whose Mellin transform satisfies the recurrent equation above.
Although these variables have been intensively studied, our new perspective,
based on a combination of probabilistic and complex analytical techniques,
enables us to derive comprehensive and substantial properties and strengthen
several results on the law of these random variables. These include smoothness,
regularity and analytical properties, large and small asymptotic behaviour,
including asymptotic expansions, bounds, and Mellin-Barnes representations for
the density and its successive derivatives. We also study the weak convergence
of exponential functionals on a finite time horizon when the latter expands to
infinity. As a result of new factorizations of the law of the exponential
functional we deliver important intertwining relation between members of the
class of positive self-similar semigroups.
  The derivation of our results relies on a mixture of complex-analytical and
probabilistic techniques.
