Three-dimensional (3D) numerical simulations of the flow around a cylinder undergoing vortex-induced vibration (VIV) at Re = 3900 (defined as Re = U∞D/ν, where U∞ is the inlet flow velocity, D is the diameter of the cylinder, and ν is the fluid kinematic viscosity) are performed using large eddy simulations (LES). Detailed mesh and time-step convergence studies are carried out for the flow past a stationary cylinder to obtain the optimal mesh and time-step. The validation studies are performed by comparing the present results with the previous published experimental data and numerical results to confirm feasibility of the present numerical model for predicting the wake characteristics. The numerical model is then applied to investigate the flow around a self-excited cylinder vibrating in the cross-flow direction at three different reduced velocities. The effects of the reduced velocity are analyzed based on the cross-flow vibration amplitudes, the power spectral densities of lift and drag coefficients, and the wake flow coherent structures. Different vortex shedding modes are identified and categorized using the iso-surfaces of the pressure coefficient and the Q-criterion. The dominant 3D wake flow features are extracted and discussed by carrying out the proper orthogonal decomposition (POD) for flow velocities on multiple two-dimensional (2D) sampling planes in the wake region. The modes on the XY planes display a traveling-wave behavior and the modes on the XZ planes show cells of streamwise vortices, which display the three-dimensionality of the wake flow.