Estimating pairwise interaction effects, i.e., the difference between the joint effect and the sum of marginal effects of two input features, with uncertainty properly quantified, is centrally important in science applications. We propose a non-parametric probabilistic method for detecting interaction effects of unknown form. First, the relationship between the features and the output is modelled using a Bayesian neural network, leveraging on the representation capability of deep neural networks. Second, interaction effects and their uncertainty are estimated from the trained model. For the second step we propose a simple and intuitive global interaction measure: Expected Integrated Hessian (EIH), whose uncertainty can be estimated using the predictive uncertainty. Two important properties of the Bayesian EIH are: 1. interaction estimation error is upper bounded by the prediction error of the neural network, which ensures interaction detection can be improved by training a more accurate model; 2. uncertainty of the Bayesian EIH is well-calibrated provided the prediction uncertainty is calibrated, which is easier to test and guarantee. The method outperforms the available alternatives on simulated and real-world data, and we demonstrate its ability to detect interpretable interactions also between higher-level features (at deeper layers of the neural network).