Abstract: Hypervolume improvement (HVI) is commonly employed in multi-objective Bayesian optimization algorithms to define the acquisition function due to its Pareto-compliant property, e.g., the expected hypervolume improvement. Instead of computing a specific probability moment of this quantity, this work aims to provide the exact expression of the probability distribution of HVI for bi-objective scenarios. We consider a bi-variate Gaussian distribution in the objective space resulting from the Gaussian process modeling. We derive the probability distribution of HVI based on a cell partition of the objective space induced by the approximation points to the Pareto front. Our exact expression is superior in numerical accuracy and computation efficiency compared to the Monte Carlo approximation.Moreover, we propose a novel acquisition function - $\varepsilon$-probability of hypervolume improvement ($\varepsilon$-PoHVI), which directly utilizes HVI's distribution function. Experimentally, we show that on many widely-applied bi-objective test problems, $\varepsilon$-PoHVI significantly improves over other closely related ones, $\varepsilon$-PoI and expected hypervolume improvement, demonstrating the benefits of using the exact computation of HVI's distribution in Bayesian optimization.