  We study glueball $G$ production in gluonic penguin decay $B\to G + X_s$,
using the next-to-leading order $b\to s g^*$ gluonic penguin interaction and
effective couplings of a glueball to two perturbative gluons. Subsequent decays
of a scalar glueball are described by using techniques of effective chiral
Lagrangian to incorporate the interaction between a glueball and pseudoscalar
mesons. Mixing effects between the pure glueball with other mesons are
considered. Identifying the $f_0(1710)$ to be a scalar glueball, we find that
both the top and charm penguin to be important and obtain a sizable branching
ratio for $B\to f_0(1710) + X_s$ of order $1.3\times 10^{-4}
(f/0.07\mbox{GeV}^{-1})^2$, where the effective coupling strength $f$ is
estimated to be $0.07$ GeV$^{-1}$ using experimental data for the branching
ratio of $f_0(1710) \to K \overline K$ based on chiral Lagrangian estimate. An
alternative perturbative QCD based estimation of $f$ is a factor of 20 larger,
which would imply a much enhanced branching ratio. Glueball production from
this rare semi-inclusive $B$ decay can be probed at the LHCb and Belle II to
narrow down the allowed parameter space. Similar branching ratio is expected
for the pseudoscalar glueball. We also briefly comment on the case of vector
and tensor glueballs.
