wave-function must be single valued function at each point of the condensate. On writing $
\psi = \psi_0 e^{i\chi}, 
$ the gauge invariant superfluid velocities can be expressed through the gradient of the phase of superfluid order parameter $\chi$ and the vector potential, ${\bf A}$ $$\label{eq:ASv}
{\bf v}_{\tau}=\frac{\hbar}{2m_{\tau}}{\bf \nabla }
\chi _{\tau}-\frac{e_{\tau}}{m_{\tau}c} {\bf A},$$ where $ e_{\tau} \equiv (e, 0 )$ is the electric charge of protons and neutrons respectively, $m_{\tau}$ is their bare mass; $\tau \in \{n, p\}$, where $n$ and $p$ stand for neutrons and protons respectively. Applying the curl operator to Eq. (\[eq:ASv\]) and implementing quantization of the circulation (the phase of the superfluid order parameter changes by $2\pi$ around a closed path) we find $$\begin{aligned}
\label{eq:AS:curlv}
\mathop{\rm curl} {\bf v_{\tau}} = \frac{\pi \hbar}{m_{\tau}}
{\bf  \nu}_{\tau} \sum_{j} \delta^{(2)} ({{\bf x}}-{{\bf x}}_{\tau j})
- \frac{ e_{\tau}}{m_{\tau}c}{\bf B} \equiv {\bf  \omega}_{\tau},\end{aligned}$$ where $\pi \hbar/m_{\tau}$ is the circulation quantum, ${\bf  \nu}_{\tau}
\equiv {\bf  \omega}_{\tau}/\omega_{\tau}$ is a unit vector along the vortex lines, ${{\bf x}}_{\tau j}$ defines the position of a vortex line in the plane orthogonal to the vector ${\bf  \nu}_{\tau}$, $\delta^{(2)}$ is a two–dimensional delta function in this plane and ${\bf B}=\mathop{\rm curl} \, {\bf A}$ is the magnetic 