Functional integration for the order parameter in $XY$ model In the continuum limit the Hamiltonian of the classical XY model is given by, ignoring the inessential constant:
$$H=\int d\vec{r}\ (\nabla\theta)^2$$
and the x-component of the order parameter is given by 
$$\langle S_x\rangle=\langle \cos\theta(\vec{r}) \rangle=\langle \cos\theta(0) \rangle=\frac{\int D\theta \cos\theta(0)e^{-\beta H}}{Z}=\operatorname{Re}\frac{\int D\theta e^{-\beta H+i\theta(0)}}{Z}$$
Going to Fourier space:
$$\theta(0)=\int \frac{d\vec{k}}{(2\pi)^d}\theta(\vec{k})$$
$$\int d\vec{r}\ (\nabla\theta)^2=\int \frac{d\vec{k}}{(2\pi)^d}k^2\theta(\vec{k})\theta(-\vec{k})=\int \frac{d\vec{k}}{(2\pi)^d}k^2\theta(\vec{k})\theta^*(\vec{k})$$
We have, using $\theta(\vec{k})=\theta_R+i\theta_I$:
\begin{equation}
\begin{aligned}
\int D\theta e^{-\beta H+i\theta(0)}&=\int D\theta_R D\theta_Ie^{-\int\frac{d\vec{k}}{(2\pi)^d}(\beta k^2\theta^2_R-i\theta_R)+(\beta k^2\theta^2_I+\theta_I)}\\
&=\int D\theta_R D\theta_Ie^{-\int\frac{d\vec{k}}{(2\pi)^d}\beta k^2(\theta_R-\frac{i}{2\beta k^2})^2+1/(4\beta k^2)+\beta k^2(\theta_I+\frac{1}{2\beta k^2})^2-1/(4\beta k^2)}\\
&=\int D\theta_R D\theta_Ie^{-\int\frac{d\vec{k}}{(2\pi)^d}\beta k^2(\theta_R-\frac{i}{2\beta k^2})^2+\beta k^2(\theta_I+\frac{1}{2\beta k^2})^2}\\
&=\int D\theta_R D\theta_Ie^{-\int\frac{d\vec{k}}{(2\pi)^d}\beta k^2(\theta_R)^2+\beta k^2(\theta_I)^2}\\
&=Z
\end{aligned}
\end{equation}
which is not correct. The above steps basically involve completing the square in Fourier space followed by shifting the functional variables. The $1/(4\beta k^2)$ shouldn't have been canceled out, but I don't see any sign problem here. Note that the $k$ integration is only over half of the $k$-space, because both $\theta_R$ and $\theta_I$ are only independently defined over half of the $k$-space due to the reality of $\theta(x)$. What I have been trying to do is to fill the gap between Eq.(19) and (20) in this XY model notes.
======================= Edit ==================================
Thanks to Shane for drawing my attention to the integration measure. In my original calculation, I directly used $\theta_R$ and $\theta_I$ as the two independent variables, which is not exactly right here. The reason is the following:
Due to the reality of $\theta(x)$, we have $\theta^*(\vec{k})=\theta(-\vec{k})$, which means $\theta_R(\vec{k})-i\theta_I(\vec{k})=\theta_R(-\vec{k})+i\theta_I(-\vec{k})$, which again means $\theta_R$ is symmetric and $\theta_I$ is anti-symmetric. Using the anti-symmetry of $\theta_I$, the $\int d\vec{k}\theta_I$ part of my above integration readily vanishes, i.e.
$$\int D\theta e^{-\beta H+i\theta(0)}=\int D\theta_R D\theta_Ie^{-\int\frac{d\vec{k}}{(2\pi)^d}(\beta k^2\theta^2_R-i\theta_R)+(\beta k^2\theta^2_I)}$$ 
and this is why there would be no completing of square for the $\theta_I$ part and extra term $1/(4\beta k^2)$ remains. In the end will have
$$\int D\theta e^{-\beta H+i\theta(0)}=Z\mathrm{exp}\left(-\int \frac{d\vec{k}}{(2\pi)^d}1/(4\beta k^2)\right)$$
where the exponent $\sim -\int^{\pi/a}_{\pi/L}dk k^{d-3}$.
For $d=1,2$, the exponent $\to -\infty$ when system size $L\to \infty$, meaning $S_x\to 0$. Whereas, for $d>2$ the order parameter takes a non-zero value.  
 A: So I think there is an issue with pulling the $Re$ of $Re(e^{i\theta})$ out of the path integral because the measure $\mathcal{D}$ has a bunch of factors of i in it after Fourier transforming.  I don't know for sure but I think this messes with your ability to do: $\theta(-k)=\theta(k)^*$(i.e. assume the field is real).
I might be completely wrong though but my reason for believing this is completing the square as follows gets the right result:
$
\int_{\infty}^{\infty} \frac{d\vec{k}}{(2\pi)^d} \beta k^2\theta(\vec{k})\theta(-\vec{k}) +i\theta(\vec{k})=
$
$
2\int_{0}^{\infty} \frac{d\vec{k}}{(2\pi)^d} \beta k^2\theta_1(\vec{k})\theta_2(\vec{k}) +i\theta_1(\vec{k})+i\theta_2(\vec{k})=
$
$
\int_{0}^{\infty} \frac{d\vec{k}}{(2\pi)^d} 2\beta k^2(\theta_1(\vec{k})+i/2\beta k^2)(\theta_2(\vec{k})+i/2\beta k^2) +1/\beta k^2
$
With $\theta_1(k)=\theta(k)$ and $\theta_2(k)=\theta(-k)$.  Assuming you can shift $\theta_1$ and $\theta_2$ independently.  Things now become:
$
2\int_{0}^{\infty} \frac{d\vec{k}}{(2\pi)^d} \beta k^2\theta_1(\vec{k})\theta_2(\vec{k}) +1/4\beta k^2
$
You can now do the integration on the second term and pull it out of the path integral.  The remain first term is the same as the $-\beta H$ so the remaining path integral becomes the the partition function $Z$ after functional integration.
You can  do the rotation:
$\theta_1=\theta_+ +\theta_-$,
$\theta_2=\theta_+ -\theta_-$ To diagonalize and do the functional integration.
Note to preform the shift independently, the following constraint must be relaxed $\theta_1=\theta_2^*$.  This why I think there is something funny about the complex analysis.  It might be related to symmetry breaking but I don't see that clearly.
A: I don't know what the author of the notes has in mind, but in any dimension, and in the absence of any direction breaking field, we must have 
$\langle \cos\theta \rangle=0$ because all directions of the order parameter are equally likely. To get spontanous symmetry breaking you need to impose a small symmetry breaking perturbation  that makes  $\langle \cos\theta \rangle$ non-zero, and then take the perturbation to zero to see if $\langle \cos\theta \rangle$ has a non-zero limit.  
I suspect that the author has ignored the effect of the singular ${\bf k}=0$ mode in doing his momentum space calculation.   
