To investigate the characteristics of the bi-stable flow at subcritical Reynolds numbers, large eddy simulation is adopted to simulate the crossflow around two tandem circular cylinders at Re = 3900. The reattachment/co-shedding bi-stability is observed in the simulations with spacing ratios (L/D, L is the center-to-center cylinder spacing and D is the diameter) of 4.5 and 4.7. Statistical analyses are performed on the hydrodynamic coefficients, time-averaged flow fields, three-dimensional characteristics, wake pattern, and vortex shedding frequencies at different spacing ratio and time period. In addition, a detailed analysis and explanation were conducted on the secondary vortices identified in the reattachment flow regime, revealing that the secondary vortices, generated from the instability of the shear layer, significantly influence the variation in vortex shedding frequency over time. The reduced-order variational mode decomposition method is employed to decompose the flow field during the flow regime transitions, unveiling their spatial and temporal features. It is revealed that the shear layer instability and the low-frequency modulation behavior are the predominant factors triggering the bi-stable phenomenon at subcritical Reynolds numbers. This study aims to uncover triggering mechanisms underlying the bi-stable phenomenon in the flow around two tandem cylinders and provides valuable insight for relevant engineering applications.