The covariance of a stationary process $X$ is diagonalized by a Fourier transform. It does not take into account the complex Fourier phase and defines Gaussian maximum entropy models. We introduce a general family of phase harmonic covariance moments, which rely on complex phases to capture non-Gaussian properties. They are defined as the covariance of $\hat{H} (L X)$, where $L$ is a complex linear operator and $\hat{H} $ is a non-linear phase harmonic operator which multiplies the phase of each complex coefficient by integers. The operator $\hat{H} (L X)$ can also be calculated from rectifiers, which relates $\hat{H} (L X)$ to neural network coefficients. If $L$ is a Fourier transform then the covariance is a sparse matrix whose non-zero off-diagonal coefficients capture dependencies between frequencies. These coefficients have similarities with high order moment, but smaller statistical variabilities because $\hat{H} (L X)$ is Lipschitz. If $L$ is a complex wavelet transform then off-diagonal coefficients reveal dependencies across scales, which specify the geometry of local coherent structures. We introduce maximum entropy models conditioned by these wavelet phase harmonic covariances. The precision of these models is numerically evaluated to synthesize images of turbulent flows and other stationary processes.