To avoid a problem associated with singularities in the discretization near the axis, we follow a method proposed by [@fuk2002].

For each run, we confirm that drag forces, numerically evaluated on the bubble surface, for both the perpendicular and parallel motions are in good agreement with the respective kinetic-energy dissipation rates, numerically integrated over the entire computational domain, normalized by the translational velocities with an error of less than $0.040\%$. Such an agreement between the surface and bulk quantities indicates that the computation is well converged and reaches to the steady state in view of the momentum and kinetic-energy budgets. Further, the drag force for the perpendicular motion shows quantitative agreement with the infinite series solution of [@bar1968] with an error of less than $0.043\%$. The drag force for the parallel motion approaches the wide-gap solution of [@mag2003] with increasing $\varepsilon$. To make sure of numerical stability and accuracy, we set the clearance parameter in a range of $10^{-3}\leq \varepsilon\leq 9$. We performed the convergence tests with varying the size $R_{\rm max}$ of the computational domain, the number $N_\xi$ of nodes in the gap between the bubble and the wall, and the number $N_\eta$ of nodes describing the bubble surface. 