In this paper, we consider numerical approximations of the binary surfactant phase-field model on complex surfaces. Consisting of two nonlinearly coupled Cahn-Hilliard type equations, the system is solved by a fully discrete numerical scheme with the properties of linearity, decoupling, unconditional energy stability, and second-order time accuracy. The IGA approach based on Loop subdivision is used for the spatial discretizations, where the basis functions consist of the quartic box-splines corresponding to the hierarchic subdivided surface control meshes. The time discretization is based on the so-called explicit-IEQ method, which enables one to solve a few decoupled elliptic constant-coefficient equations at each time step. We then provide a detailed proof of unconditional energy stability along with implementation details, and successfully demonstrate the advantages of this hybrid strategy by implementing various numerical experiments on complex surfaces.(c) 2023 Elsevier B.V. All rights reserved.