In this paper, we work in the setting of Bessel operators and Bessel Laplace equations studied by Weinstein, Huber, and the harmonic function theory in this setting introduced by Muckenhoupt–Stein, especially the generalised Cauchy–Riemann equations and the conjugate harmonic functions. We provide the equivalent characterizations of product Hardy spaces associated with Bessel operators in terms of the Bessel Riesz transforms, non-tangential and radial maximal functions defined via Poisson and heat semigroups, based on the atomic decomposition, the generalised Cauchy–Riemann equations, the extension of Merryfield’s result which connects the product non-tangential maximal function and area function, and on the grand maximal function technique which connects the product non-tangential and radial maximal function. We then obtain directly the decomposition of the product BMO space associated with Bessel operators.