We consider a low complexity polynomial expansion (PE) detector in a massive multiple-input multiple-output (MIMO) uplink channel. In contrast to most massive MIMO systems in the literature, where single antenna user equipments (UEs) are assumed, multiple antenna UEs are employed in this paper. Moreover, the channel between a base station (BS) and a UE is a jointly correlated Rician fading channel. The PE detector reduces the computational complexity of the minimum mean square error (MMSE) detector by replacing the matrix inversion with an approximate matrix polynomial. The coefficients of the approximate matrix polynomial are computed from the deterministic equivalents of the moments of the channel Gram matrix. We use operator-valued free probability, which is a more general version of free probability, to derive the deterministic equivalents. In particular, we use the operator-valued moment-cumulant formula. The proposed low complexity PE detector is easy to compute. Simulation results show that the proposed detector can achieve performance close to the MMSE detector.