We develop a coupled-cluster theory for bosonic mixtures of binary species in external traps, providing a promising theoretical approach to demonstrate highly accurately the many-body physics of mixtures of Bose–Einstein condensates. The coupled-cluster wavefunction for the binary species is obtained when an exponential cluster operator eT, where T = T(1) + T(2) + T(12) and T(1) accounts for excitations in species-1, T(2) for excitations in species-2, and T(12) for combined excitations in both species, acts on the ground state configuration prepared by accumulating all bosons in a single orbital for each species. We have explicitly derived the working equations for bosonic mixtures by truncating the cluster operator up to the single and double excitations and using arbitrary sets of orthonormal orbitals for each of the species. Furthermore, the comparatively simplified version of the working equations are formulated using the Fock-like operators. Finally, using an exactly solvable many-body model for bosonic mixtures that exists in the literature allows us to implement and test the performance and accuracy of the coupled-cluster theory for situations with balanced as well as imbalanced boson numbers and for weak to moderately strong intra- and interspecies interaction strengths. The comparison between our computed results using coupled-cluster theory with the respective analytical exact results displays remarkable agreement exhibiting excellent success of the coupled-cluster theory for bosonic mixtures. All in all, the correlation exhaustive coupled-cluster theory shows encouraging results and could be a promising approach in paving the way for high-accuracy modeling of various bosonic mixture systems.