We investigate the prospects for the discovery of a neutral Higgs boson produced in association with a $b$ quark, followed by the Higgs decay into a pair of bottom quarks, $pp to bphi^0 to b bbar{b} +X$, at the CERN Large Hadron Collider (LHC) within the framework of unified supersymmetric models. The Higgs boson $phi^0$ can be a heavy scalar $H^0$ or a pseudoscalar $A^0$. Furthermore, this direct discovery channel is compared with the indirect Higgs searches in the rare decay $B_s to mu^+mu^-$ at hadron colliders. Promising results are found for the minimal supergravity (mSUGRA) model, the anomaly mediated supersymmetry breaking (AMSB) model, and the gauge mediated supersymmetry breaking (GMSB) model. We find that the indirect search for $B(B_s to mu^+mu^-) ge 5times 10^{-9}$ is complementary to the direct search for $bphi^0 to bbbar{b}$ with $sqrt{s} = 14$ TeV and an integrated luminosity ($L$) of 300 fb$^{-1}$. In the AMSB and GMSB models, $bphi^0 to bbbar{b}$ with $L = 300$ fb$^{-1}$ covers a larger area in the parameter space than $B(B_s to mu^+mu^-) ge 5times 10^{-9}$. In addition, we present constraints from $b to sgamma$ and muon anomalous dipole moment ($Delta a_mu$) on the parameter space.