Wake of two side-by-side square cylinders was investigated through direct numerical simulation at low Reynolds numbers (16~200). The gap between the two cylinders varied from 0-10D, where D is the dimension of the square cylinder (edge length). 9 different wake patterns and their dependency on both Reynolds number and gap spacing were identified and analysed. A system classification map, demarcated by Reynolds number and gap ratio g* (g/D, where g is the gap spacing between 2 cylinders), was derived for these 9 wake modes. Steady-state wake (mode I) was observed when Reynolds number is lower than the critical Reynolds number, which depends on g*. For gap ratio less than 0.7, only single vortex street was observed. The single vortex street wake can be either symmetric and periodic (mode II), or asymmetric and periodic (mode III), or irregular (mode IV). In this gap ratio range (less than 0.7), shedding frequency decreases with gap ratio due to the damping role of the gap flow. For gap ratio larger than 0.7, two vortex streets were also observed. For gap ratio larger than 1, only two vortex streets were observed. Vortex shedding can be either synchronized and in-phase (mode V), synchronized and anti-phase (mode VI), in-phase dominated with low frequency modulation (mode VII), anti-phase dominated with low frequency modulation (mode VIII), asymmetric synchronized anti-phase (mode IX), or irregular (mode IV). For gap ratio larger than 4, only synchronized anti-phase mode was observed under the conditions of this study. In the two vortex streets regime, shedding frequency is higher than that of a single cylinder, due to stronger gap flow than that in the freestream side. The impact of gap ratio and Reynolds number on the drag and lift forces was also studied.2