Why are superfluid vortex lattices stable? Both (a) neutral superfluids that are externally rotated, and (b) type-II superconductors (i.e. charged superfluids) under applied magnetic fields between the critical fields $h_{c1}$ and $h_{c2}$, have topological-defect vortex excitations (which are point defects in 2D and line defects in 3D), around which the phase of the order parameter winds by an integer multiple of $2 \pi$.  This integer multiple gives the vortex's topological charge. In the neutral-superfluid case, these vortices carry angular momentum that is quantized to an integer multiple of $2 \pi \hbar$, and in the superconducting case they carry magnetic flux which is quantized to an integer multiple of the magnetic flux quantum.  Topological charges greater than $\pm1$ are possible but unstable, and disintegrate into multiple vortices with charge $\pm 1$, and vortices with opposite charge annihilate, so in equilibrium, all the vortices have the same charge.  Vortices with the same charge repel, so they form a triangular vortex lattice.
My question is, why is this vortex lattice thermodynamically stable? Why don't the repelling vortices all get pushed to the sample boundary?
Intuitively, we imagine that a system with a nonzero net electric charge density does not have a well-defined thermodynamic limit.  For example, consider a classical electric conductor.  Any net charge density that we place on it goes to the boundary, so there is no bulk thermodynamic limit with nonzero charge density. (In fact, since the total charge grows as the volume of the system but the boundary only grows as the area, the surface charge density grows unboundedly with system size, and there's no well-defined boundary thermodynamic limit either).  This intuition can be made more rigorous by considering a gas of particles (either classical or quantum) interacting via the Coulomb repulsion.  Naively one might expect there to never be a thermodynamic limit, since the Coulomb interaction is long-range (its integral over all space diverges).  But it turns out that if the net charge is zero, then the positive and negative charges screen each other and put an exponentially decaying envelope on the effective electric interaction, so there is in fact a well-defined thermodynamic limit.  But if there is a net charge density imbalance, then there is no thermodynamic limit - intuitively, all the excess charge gets pushed to the boundary.  All this is made rigorous in Section V of this review.
It seems to me that the case of a vortex lattice in a superfluid is mathematically analogous to an electric system with a net charge density (physically set by either the rotation or the magnetic field applied to the superfluid).  The vortices interact via a logarithmic effective Coulomb interaction by charge-vortex duality.  And crucially, only the vortices carry topological charge - there's no background "charge" distribution (that I can think of) to screen them, so it seems like they should feel each other's full Coulomb interaction.  So why do they form a thermodynamically stable bulk lattice instead of pushing each other to the boundary, like classical charges in a conductor do?
(Caveat: in a type-II superconductor, the logarithmic interaction does get exponentially screened at distances much longer than the London penetration depth $\lambda$.  This may be enough to thermodynamically stabilize the Abrikhosov vortex lattice, although that's not obvious to me; it seems reasonable at applied fields just above $h_{c1}$, where the lattice is very dilute, but what about at fields near $h_{c2}$, where I believe the inter-vortex spacing should be shorter than $\lambda$ so the screening shouldn't matter?  And in any case, I believe that for a neutral superfluid, there's no long-distance screening and the vortex interaction remains logarithmic, so why don't the vortices push each other to the boundary in the neutral case?)
 A: If you can forgive the hand-waviness: I do not have a proper answer, but I can comment as it follows.
Usually these systems have a sample boundary because they are harmonically trapped. If you consider the ideal case in which the trap is completely lowered, the resulting ground state would consist of an infinite triangular vortex lattice, with vortex density given by Feynman's rule (at least as long as you consider mean field theory and keep the density finite); and this state is certainly stable. 
I believe the introduction of an harmonic potential affects the particle density and not the overall vortex density (Right? I am not too confident on this statement. Certainly the vortex density is affected locally.). I usually think that what prevents the trapped vortices from escaping the sample boundary are other ghost vortices sitting where the density is zero (i.e. outside the sample boundary).
This might give an answer in mean field: it is certainly a different story when considering the microscopic theory.
Does this somehow make sense?
A: "My question is, why is this vortex lattice thermodynamically stable? Why don't the repelling vortices all get pushed to the sample boundary?"
Equal charge vortices are not really repulsive: this is a misunderstanding due to the analogy with electric systems. Rather, they have the tendency to orbit one around the other. This is why the Abrikosov lattice rotates and vortices are not expelled. On the other hand, opposite charge vortices will maintain their distance, but travel in parallel lines.
This is the behaviour of a dissipation-free system. If you introduce dissipation (i.e. the system is at finite temperature and out of equilibrium), then the "equal charge" vortices will spiral in (or out) till a new equilibrium is realized. The opposite charge vortices will tend to decrease their distance or even collide and annihilate in the presence of dissipation.
PS: the vortex lattice thermodynamically stable depending on the experiment you are performing. If you create a vortex state and then you stop the bucket (or change its angular velocity), then the lattice is not thermodynamically stable. If you keep the bucket spinning at a constant angular velocity, then the thermodynamic state that is realized at the end of the initial dissipative relaxation is, by definition, the one that has a certain angular momentum, and this can be achieved only with a certain number of vortices.
