  In this paper, we consider the numerical approximations for the commonly used
binary fluid-surfactant phase field model that consists two nonlinearly coupled
Cahn-Hilliard equations. The main challenge in solving the system numerically
is how to develop easy-to-implement time stepping schemes while preserving the
unconditional energy stability. We solve this issue by developing two linear
and decoupled, first order and a second order time-stepping schemes using the
so-called "Invariant Energy Quadratization" approach for the double well
potentials and a subtle explicit-implicit technique for the nonlinear coupling
potential. Moreover, the resulting linear system is well-posed and the linear
operator is symmetric positive definite. We rigorously prove the first order
scheme is unconditionally energy stable. Various numerical simulations are
presented to demonstrate the stability and the accuracy thereafter.
