  Einstein-Maxwell-scalar models allow for different classes of black hole
solutions, depending on the non-minimal coupling function $f(\phi)$ employed,
between the scalar field and the Maxwell invariant. Here, we address the linear
mode stability of the black hole solutions obtained recently for a quartic
coupling function, $f(\phi)=1+\alpha\phi^4$ [1]. Besides the bald
Reissner-Nordstr\"om solutions, this coupling allows for two branches of
scalarized black holes, termed cold and hot, respectively. For these three
branches of black holes we calculate the spectrum of quasinormal modes. It
consists of polar scalar-led modes, polar and axial electromagnetic-led modes,
and polar and axial gravitational-led modes. We demonstrate that the only
unstable mode present is the radial scalar-led mode of the cold branch.
Consequently, the bald Reissner-Nordstr\"om branch and the hot scalarized
branch are both mode-stable. The non-trivial scalar field in the scalarized
background solutions leads to the breaking of the degeneracy between axial and
polar modes present for Reissner-Nordstr\"om solutions. This isospectrality is
only slightly broken on the cold branch, but it is strongly broken on the hot
branch.
