Two-fluid plasma flow equations describe the plasma flow of ions and electrons with different densities, velocities, and pressures. The system of equations has flux consisting of three independent components, one each for ions and electrons flows and a linear Maxwell`s equation flux for the electromagnetic fields. The coupling of these components is via source terms.

This article introduces second-order finite difference schemes aimed at maintaining the consistent evolution of divergence constraints on both electric and magnetic fields. The core concept involves developing a numerical solver for Maxwell`s equations, utilizing a multidimensional Riemann solver at vertices to ensure discrete divergence constraints. To handle the fluid components, we employ an entropy-stable discretization approach. Our proposed schemes are collocated, boasting second-order accuracy, entropy stability, and ensuring the divergence-free evolution of the magnetic field. Time discretization is accomplished through explicit and Implicit-Explicit (IMEX) schemes. To demonstrate the accuracy, stability, and divergence constraint-preserving ability of the proposed schemes, we present several test cases in one and two dimensions. We also compare numerical results with schemes with no divergence cleaning and perfectly hyperbolic Maxwell (PHM) equations based divergence cleaning schemes for Maxwell`s equations.