  In this paper, we consider the numerical approximation of a time-fractional
stochastic Cahn--Hilliard equation driven by an additive fractionally
integrated Gaussian noise. The model involves a Caputo fractional derivative in
time of order $\alpha\in(0,1)$ and a fractional time-integral noise of order
$\gamma\in[0,1]$. The numerical scheme approximates the model by a piecewise
linear finite element method in space and a convolution quadrature in time (for
both time-fractional operators), along with the $L^2$-projection for the noise.
We carefully investigate the spatially semidiscrete and fully discrete schemes,
and obtain strong convergence rates by using clever energy arguments. The
temporal H\"older continuity property of the solution played a key role in the
error analysis. Unlike the stochastic Allen--Cahn equation, the presence of the
unbounded elliptic operator in front of the cubic nonlinearity in the
underlying model adds complexity and challenges to the error analysis. To
overcome these difficulties, several new techniques and error estimates are
developed. The study concludes with numerical examples that validate the
theoretical findings.
