The interaction of qubits with quantized modes of electromagnetic fields has been largely addressed in the quantum optics literature under the rotating wave approximation (RWA), where rapid oscillating terms in the qubit-mode interaction picture Hamiltonian can be neglected. At the same time, it is generally accepted that provided the interaction is sufficiently strong or for long times, the RWA tends to describe physical phenomena incorrectly. In this work, we extend the investigation of the validity of the RWA to a more involved setup where two qubit-mode subsystems are brought to interaction through their harmonic coordinates. Our treatment is all analytic thanks to a sequence of carefully chosen unitary transformations which allows us to diagonlize the Hamiltonian within and without the RWA. By also considering qubit dephasing, we find that the purity of the two-qubit state presents non-Markovian features which become more pronounced as the coupling between the modes gets stronger and the RWA loses its validity. In the same regime, there occurs fast generation of entanglement between the qubits which is also not correctly described under the RWA. The setup and results presented here clearly show the limitations of the RWA in a scenario amenable to exact description and free from numerical uncertainties. Consequently, it may be of interest for the community working with cavity or circuit quantum electrodynamic systems in the strong coupling regime.