The XY-model Hamiltonian is the following,
$${\cal H}~=~-J\sum_{\langle i,j\rangle} \cos (\theta_i -\theta_j).$$
The Goldstone mode corresponds to term $(\nabla \theta)^2$ in the effective Lagrangian.
Then what's the form of effective Lagrangian that produces both Goldstone mode and topological excitation (vortices)?
And how to derive the effective Lagrangian from the XY-Hamiltonian?
This is quite a standard textbook example. Good treatments can, for example, be found in Altland (Condensed Matter Field Theory) or Nagaosa (Quantum Field Theory in Condensed Matter Physics). However, the basic reasoning can be understood quite easily:
Consider the change in angle produced as we perform a closed loop around some point. Naively we obtain \begin{align} \Delta \theta = \oint d \vec r \cdot \nabla \theta(\vec r) = \int d^2 r \, \vec e_z \cdot \vec \nabla \times \nabla \theta(\vec r) = 0, \end{align} where Stokes theorem was used in the second step and the fact that $\vec \nabla \times \nabla f = 0$ for any scalar function f. However, we know that due to the compact nature of our angle variable $\theta(\vec r)$, we don't have to get back exactly to the same angle, but have more freedom and are allowed to have $\Delta \theta = 2\pi n$ with $n \in \mathbb{Z}$. When Taylor-expanding the original cosine, we have lost the information about the compact nature of the $\theta$-variables and thus have to reinclude it manually.
To include this additional freedom, we have to generalize the phase difference field $\vec u$ in the form $\vec u = \nabla \theta(\vec r) + \vec \nabla \times \vec A(\vec r)$, where the vector field $\vec A(\vec r) = \psi(\vec r) \vec e_z$ can have only an $z$-component in order for $\vec u$ to lie within the xy-plane. Then \begin{align} 2 \pi n = \oint d \vec r \cdot \vec u(\vec r) = \int d^2 r \, \vec e_z \cdot \vec \nabla \times \vec u(\vec r) = -\int d^2 r \, \vec \nabla^2 \psi(\vec r) \end{align} where $\vec \nabla^2 = (\vec \nabla \cdot \vec \nabla)$ is the Laplace-operator and the vector identity $\vec \nabla \times (\vec \nabla \times \vec A) = \nabla (\vec \nabla \cdot \vec A) - (\vec \nabla \cdot \vec \nabla) \vec A$ was used whose first term vanishes for $\vec A = \psi(\vec r) \vec e_z$ since $\psi(\vec r)$ has no $z$-dependence. Consequently, one finds \begin{align} \vec \nabla^2 \psi(\vec r) = -2\pi \sum_i n_i \delta(\vec r - \vec r_i) \end{align} and we interpret this by saying that the field $\psi(\vec r)$ describes vortices with winding numbers $n_i$ centered at positions $\vec r_i$. The above equation for $\psi(\vec r)$ is the defining equation of linear combinations of fundamental solutions of the Laplace equation (in 2D) which are given by \begin{align} \psi(\vec r) = - 2 \pi \sum_i n_i \log(\vert \vec r - \vec r_i\vert). \end{align} Plugging the phase difference field into the effective Hamiltonian, one finds \begin{align} H = -\frac{J}{2} \int d^2 r \, \vec u^2 = -\frac{J}{2} \int d^2 r \, (\nabla \theta(\vec r))^2 + (\vec \nabla \times \psi(\vec r) \vec e_z)^2, \end{align} where the mixed term involving both the $\theta$ and the $\psi$ field vanished after a partial integration. The first part corresponds to the spin-wave excitations you mentioned, while the second encodes the vortices. It is straight-forward to show that $ (\vec \nabla \times \psi(\vec r) \vec e_z)^2 = (\nabla \psi)^2 = -\psi (\vec \nabla \cdot \vec \nabla) \psi$ after a partial integration, such that the topological part one ends up with is given by \begin{align} H_\text{top} = \frac{J}{2} \sum_{ij} n_i n_j \frac{\log(\vert \vec r_i - \vec r_j\vert)}{2\pi}. \end{align} The passage to the Lagrangian formulation of the problem is trivial since no canonical momenta are involved.
