I am trying to follow this paper and track the dynamics of vortex motion on a discrete (square) lattice. The idea is to simulate the time evolution of the Gross-Pitaevskii (GP) equation, which reads (in rescaled units)
$$i \partial_t \psi = - \nabla^2 \psi + |\psi|^2 \psi.$$
As initial condition one takes a wavefunction of the form of Eq. (10) in the above-mentioned paper,
$$\psi(z) = \prod_{z_+ \in [Z_+]} \frac{(z-z_+)}{|z-z_+|} \prod_{z_- \in [Z_-]} \frac{(z^*-z_-^*)}{|z-z_-|}$$
where $z=x+iy$ are complex coordinates. This wavefunction describes a set of vortices (V) with coordinates $z_+$ in $Z_+$ and a set of anti-vortices (AV) with coordinates $z_-$ in $Z_-$. This is equivalent to using a wavefunction of the type $e^{i \phi}$ (polar coords.) for vortices and $e^{-i \phi}$ for anti-vortices, with $\phi$ the polar angle. One can calculate the velocity as gradient of the phase of the complex wavefunction. An example of the associated velocity field, for a V at (x,y)=(0,2.5) and AV at (0,-2.5) is shown below.

This should correspond to Fig. 1 in the paper. However, there are several things I do not understand.
First, why do the authors say that they cannot look at the dynamics of just one vortex, but must include an anti-vortex in the initial conditions to satisfy the periodic boundary conditions?
Second, apart from the initial V-AV pair, they also seem to include "a few mirrow images with respect to the boundary". However, they do not specify how many nor where they are positioned. (I assume that the mirrow image of a V is and AV and vice-versa, please correct me if I'm wrong).
Lastly, they mention "evolving the system in imaginary time" in order to "cool it down" (although no temperature is included in the model!) and then "turning on a superflow" and continuing the evolution in real-time. I though all was needed was the discretization of the GP equation and following its evolution given the initial wavefunction and periodic boundary conditions ($\psi(\text{left margin})=\psi(\text{right margin})$ and $\psi(\text{top margin})=\psi(\text{bottom margin})$).
Edit
Following the suggestion of @BebopButUnsteady, I repeated the 'unit' of my system, the V-AV pair, until obtaining the following tiling pattern (vortices are circles, AVs are squares)

Let us first look at the real (left) and imaginary (right) parts of the wavefunctions for the initial V-AV pair:

One can see that, especially in the imaginary part, the values on the boundaries are not equal. Now we start adding copies of the "unit cell" along the x axis, and plot a cut though the imaginary part of the wf:

It appears that, as we increase the number of copies, the values at the ends get closer and closer to 0, which would be the ideal periodic case. Of course this is just numerics, I have no formal proof yet that this procedure actually converges.