The three-dimensional flow past two identical spherical bubbles moving side by side in a viscous fluid is studied numerically by solving the full Navier–Stokes equations. The bubble surface is assumed to be clean so that the outer flow obeys a zero-shear-stress condition. The present study describes the interaction between the two bubbles over a wide range of Reynolds number ($0.02\,{\le}\,Re\,{\le}\,500$, $Re$ being based on the bubble diameter and rise velocity), and separation $S$ ($2.25\,{\le}\,S\,{\le}\,20$, S being the distance between the bubble centres normalized by the bubble radius). The flow structure, the vorticity field, the sign of the interaction force and the magnitude of the drag and lift forces are analysed; in particular the latter are compared with analytical expressions available in the potential flow limit and in the limit of low-but-finite Reynolds number. This study sheds light on the role of the vorticity generated at the bubble surface in the interaction process. When vorticity remains confined in a boundary layer whose thickness is small compared to the distance between the two bubbles, the interaction is dominated by the irrotational mechanism that results in an attractive transverse force. In contrast, when viscous effects are sufficiently large, the vorticity field about each bubble interacts with that existing about the other bubble, resulting in a repulsive transverse force. Computational results combined with available high-Reynolds-number theory provide empirical expressions for the drag and lift forces in the moderate-to-large Reynolds number regime. They show that the transverse force changes sign for a critical Reynolds number whose value depends on the separation. Using these computational results it is shown that, depending on their initial separation, freely moving bubbles may either reach a stable equilibrium separation or move apart from each other up to infinity.