  In this paper we establish the product Hardy spaces associated with the
Bessel Schr\"odinger operator introduced by Muckenhoupt and Stein, and provide
equivalent characterizations in terms of the Bessel Riesz transforms,
non-tangential and radial maximal functions, and Littlewood--Paley theory,
which are consistent with the classical product Hardy space theory developed by
Chang and Fefferman. Moreover, in this specific setting, we also provide
another characterization via the Telyakovski\'i transform, which further
implies that the product Hardy space associated with this Bessel Schr\"odinger
operator is isomorphic to the subspace of suitable "odd functions" in the
standard Chang--Fefferman product Hardy space. Based on the characterizations
of these product Hardy spaces, we study the boundedness of the iterated
commutator of the Bessel Riesz transforms and functions in the product BMO
space associated with Bessel Schr\"odinger operator. We show that this iterated
commutator is bounded above, but does not have a lower bound.
