Presentation of the Problem
I am making a measurement which yields me the complex quantity:
\begin{equation} S_k = e^{i \cdot \phi} + X_k + i \cdot Y_k \quad \text{with} \quad X_k, Y_k \stackrel{iid}{\sim} \mathcal{N}\left(0, \sigma_c^2\right) \end{equation}
which is affected by a bivariate, independent, Gaussian noise. My main objective here is to compute the MSE that may be reached on the estimate of $\varphi$ - let us call it $\theta$ - as a function of $\sigma_c$. This MSE is defined as:
\begin{equation} MSE = \mathbf{E}\left((\varphi - \theta)^2\right) \end{equation}
Each ${S_k}$ can be seen as a point in the complex plane of coordinates:
\begin{equation}\label{eq:xy_distri} \begin{cases} x = \Re{\left({S_k}\right)} = cos\ \varphi + X_k \sim \mathcal{N}(cos \varphi, \sigma_c^2)\\ y = \Im{\left({S_k}\right)} = sin\ \varphi + Y_k \sim \mathcal{N}(sin \varphi, \sigma_c^2) \end{cases} \end{equation}
And we can write $f(x,y)$ as the bivariate normal distribution:
\begin{equation} f(x,y) = \frac{1}{2 \cdot \pi \cdot \sigma_c^2} \cdot e^{-\frac{1}{2 \cdot \sigma_c^2}\cdot ((x - cos \varphi)^2+(y-sin \varphi)^2)} \end{equation}
We can use polar coordinate as:
\begin{equation} g_\varphi(\theta) = \int_{0}^{+\infty}f(\theta, r) rdr \quad \text{with} \quad \begin{cases} x = r \cdot cos\ \theta\\ y = r \cdot sin\ \theta \end{cases} \end{equation}
Yielding ultimately to:
\begin{equation} \boxed{ g(\theta) = \frac{1}{2\cdot \pi} \cdot e^{-\frac{1}{2 \sigma_c^2}} + \frac{1}{2 \cdot \sqrt{2 \cdot \pi} \cdot \sigma_c} \cdot cos(\theta-\varphi) \cdot e^{-\frac{sin^2(\theta - \varphi)}{2\cdot \sigma_c^2}} \cdot erfc\left(-\frac{cos(\theta - \varphi)}{\sigma_c \cdot \sqrt{2}}\right) } \end{equation}
I checked the latter equation in simulations as can be seen in the Figure below:
Wherein the blue histogram is a simulation of many ($\sim 10^7$) $S_k$ from which the angle was retrieved.
Question
Given the afore-mentioned $g(\theta)$, how can I compute the MSE on $\theta$ with respect to $\varphi$ estimation? I think that the usual expected value and variance formulae do not apply here given the "circular" nature of $\theta$ i.e. the fact that the $\left[0;2\cdot\pi\right]$ folds on itself.
I would greatly appreciate some help on this issue, even only as keywords or tracks to follow. Do not hesitate if you think that further details are requested either. I tried to keep it short and did not develop all the calculations here, but I may provide additional details upon request.