The constant is $C = \frac{\pi}{4}$.
The $2N$-dimensional vector composed of the real and imaginary parts of the coefficients of a normalized state vector is uniformly distributed on an $S^{2N-1}$. Thus the needed expectation value is:
$$\frac{C}{N} = \frac{\int_{S^{2N-1}}|(c_{1x} + i c_{1y}) (c_{2x} - i c_{2y})| d\mu(S^{2N-1})}{\mathrm{Vol}(S^{2N-1})}$$
Integrals of homogeneous functions over spheres can be traded by Gaussian integrals over the ambient Euclidean space the with proper normalizations; since the angular integrations are the same and we need only to normalize the different radial integrals. In our case we get:
$$\frac{C}{N} = \frac{\int_{\mathbb{C}^N}|(c_{1x} + i c_{1y}) (c_{2x} - i c_{2y})| e^{-\sum_j|c_j|^2}\prod_j d\mathrm{Re}(c_j)d\mathrm{Im}(c_j)}{\int_{\mathbb{C}^N}(\sum_j|c_j|^2e^{-\sum_j|c_j|^2}\prod_j d\mathrm{Re}(c_j)d\mathrm{Im}(c_j)}$$
The denominator is the unique symmetric polynomial of the same homogeneity degree as the numerator, thus gives rise to the same radial integration value.
The Gaussian integral in the numerator is decomposable to two types of integrals and in the one in the denominator denominator is a sum of N similar complex gaussian integrals, thus we get:
$$\frac{C}{N} =\frac{\left(\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \sqrt{c_x^2 + c_y^2} e^{-c_x^2 - c_y^2}\right)^2 \left(\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} e^{-c_x^2 - c_y^2}\right)^{N-2}}{N \left(\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} (c_x^2 + c_y^2) e^{-c_x^2 - c_y^2}\right) \left(\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} e^{-c_x^2 - c_y^2}\right)^{N-1}}= \frac{\pi}{4N}$$
(In the above, the integrals can be solved in polar and cartesian coordinates yielding Gamma functions of intgral or half integral arumements)
Remark:
The equivalence of the integration formulas is based on the following Dirac delta function limit representation:
$$ \delta(y) = \lim_{\epsilon \rightarrow 0} \frac{1}{\sqrt{4\pi \epsilon}} e^{-\frac{-y^2}{4\epsilon}}$$
For a homogeneous function $f(x): \mathbb{R}^{2N} \rightarrow \mathbb{R}$:
$$f(ax) = a^{\alpha} f(x)$$
Denoting: $ r^2 = \sum_{j=1}^{2N} x_j^2$, the ratio of integrals:
$$\frac{\int_{\mathbb{R}^{2N}} f(x) \frac{e^{-\frac{r^2-1}{4\epsilon}}}{\sqrt{4\pi \epsilon}}d\mu_L(\mathbb{R}^{2N})}{\int_{\mathbb{R}^{2N}} r^{\alpha} \frac{e^{-\frac{r^2-1}{4\epsilon}}}{\sqrt{4\pi \epsilon}}d\mu_L(\mathbb{R}^{2N})}$$
is invariant under the scaling $x_j \rightarrow \sqrt{\frac{\epsilon}{\epsilon’}}x_j$, and the integral conserves its form with $\epsilon$ replaced by $\epsilon’$.
Therefore, the ratio of integrals is independent of $\epsilon$. By taking the limit $\epsilon = \frac{1}{4}$, we get our Gaussian integral. By taking the limit $\epsilon \rightarrow 0$, we get a delta function concentrating the measure to the spherical shell:
$$r^2=1$$
in the Euclidean space: $\mathbb{R}^{2N}$, which the integral over the $2N-1$ dimensional sphere.