We aim to investigate the aerodynamic and acoustics characteristics of a twin square jet using an implicit Large Eddy Simulation (ILES). A screeching cold jet condition, a nozzle pressure ratio (NPR) of 3.0, is considered to simulate a coupled twin-jet. A second-order central scheme with a modified version of Jameson’s artificial dissipation is adopted to damp numerical oscillations and to mimic the effect of small-scale turbulence without an explicit subgrid-scale (SGS) model. Numerical results show that the overall trends of time-averaged streamwise velocity profiles are similar to the experimental data, with the largest differences observed at locations associated with the presence of the shock-cell structures. A detailed investigation of the flow fluctuations in jet shear layers is performed. The amplitude of the velocity fluctuations is highly dependent on the location of the shear layers with respect to the twin-jet configuration (upper, lateral, or inner). The coupling mode of twin jets associated with the screech tone is determined as a symmetrical flapping mode be a two-points spacetime cross-correlation analysis. The overall trends of near-field pressure fluctuation spectra by LES agree well with the experimental results in both upstream and downstream regions. Near-field pressure fluctuation spectra by ILES agree well with the experimentally obtained spectra at different locations in the nozzle exit plane as well as at several downstream locations in the near-field acoustic region. The highest screech tone is observed at the inter-nozzle region where superposition of in-phase waves and standing waves are found. Fourier phase and amplitude fields at the fundamental frequency also confirm the symmetrical flapping mode of the twin jets by showing in-phase relations of hydrodynamic/acoustic waves and noise directivities.