angle difference is to be computed so as to lie within the interval $(-\pi,\pi]$.

To look for vortex lattice ordering we compute the vortex structure function $S({\bf k}_\perp)$ which measures correlations between vortices within the same $xy$ plane, $$S({\bf k}_\perp)={1\over fL^2L_z}\sum_{{\bf r}_\perp,{\bf r}_\perp^\prime,z}{\rm e}^{i{\bf k}_\perp\cdot({\bf r}_\perp-{\bf r}_\perp^\prime)}\langle n_z({\bf r}_\perp,z)n_z({\bf r}_\perp^\prime,z)\rangle
\label{eSk}$$ where ${\bf r}_\perp\equiv (x,y)$ denotes the coordinates in the $xy$ plane and similarly ${\bf k}_\perp \equiv (k_x,k_y)$. The presence of a vortex lattice will be indicated by the appearance of sharp peaks in $S({\bf k}_\perp)$ at reciprocal lattice vectors ${\bf K}$.

To look for a possible vortex glass phase, in which vortex lines are frozen into a disordered configuration and so $S({\bf k}_\perp)$ displays no signal of the ordering, we consider the helicity modulus.[@OlssonVG] To simulate the Hamiltonian of Eq. (\[eH\]) on a finite size grid, we use fluctuating twist boundary conditions, defined by, $$\theta ({\bf r}_i+L_\mu\hat\mu)-\theta({\bf r}_i)=\Delta_\mu\enspace,$$ where $\Delta_\mu$, the total phase angle twist across the system in direction $\hat \mu$, is taken as a thermally fluctuating degree of freedom. We then compute the histogram $P(\Delta_\mu)$ of values that $\Delta_\mu$ takes during the course of the simulation (averaging over the twists in the transverse directions) and define the 