Table Of ContentStructure and Magnetization of Two-Dimensional Vortex Arrays in the Presence of
Periodic Pinning
Toby Joseph∗ and Chandan Dasgupta†
Centre for Condensed Matter Theory, Department of Physics, Indian Institute of Science, Bangalore 560012, India
and
3
0 Condensed Matter Theory Unit, JNCASR, Bangalore 560064, India
0
2
n magnetic induction that determines the areal density ρ0
Ground-state properties of a two-dimensional system of
a of vortices (ρ0 = B/Φ0). The filling ratio, n, defined as
J superconducting vortices in the presence of a periodic array n B/B ,measuresthe commensurabilityofthe vortex
ofstrongpinningcentersarestudiedanalytically andnumer- ≡ φ
3 systemwiththeunderlyingpinlattice. Theinterplaybe-
ically. The ground states of the vortex system at different
1 tween the lattice constant of the pin array (determined
filling ratios are obtained usinga simple geometric argument
] under the assumption that the penetration depth is much byBφ)andtheintervortexseparation(determinedbyB)
n canleadtoavarietyofinterestingeffectsinsuchsystems.
smaller than the spacing of the pin lattice. The results of
o Some of these effects have been observed in recent
this calculation are confirmed by numerical studies in which
c experiments. Imaging experiments using various tech-
simulated annealing is used to locate the ground states of
- niques such as Lorentz microscopy8 and scanning Hall
r the vortex system. The zero-temperature equilibrium mag-
p probe microscopy5,7 have shown the formation of or-
netization as a function of the applied field is obtained by
u numerically calculating the energy of the ground state for a deredstructuresofthevortexsystematlowtemperatures
s
. largenumberofcloselyspacedfillingratios. Theresultsshow for commensurate values of n. Magnetization measure-
at interesting commensurability effects such as plateaus in the ments3,4,6 in the irreversible (vortex solid) regime have
m B H diagram at simple fractional filling ratios. demonstratedtheoccurrenceofanomaliesatcertainhar-
−
monics of B . The effectiveness of pinning at integral
- φ
d 74.78.-w,74.25.Qt,74.25.Ha values of n has been found6,10–13 to produce regularly
n spacedsharpminima in the resistivityversusfieldcurve.
o
A pinning-induced reconfiguration of the vortex lattice
c has been observed13 in a thin-film superconductor with
[ I. INTRODUCTION
a rectangular array of magnetic dots. Some of these ef-
2 fectshavebeenstudiedtheoretically,usinganalytic14and
v In the mixed phase of type-II superconductors, mag- numerical15–20 methods. Experimental realizations of a
1 netic fluxpenetratesthe sampleinthe formofquantized
system of interacting “particles” in the presence of an
3 vortex lines1. The amountof flux carriedby each vortex
external periodic potential are also obtained in colloidal
5
2 line is equal to the basic flux quantum Φ0 = hc/(2e) = suspensions in interfering laser fields21, and in periodic
1 2.07 10−7G.cm2. Thesevortexlinesformaspecialphys- arrays of optical traps22.
×
2 ical system knownas “vortex matter”. In the absence of In this paper, we have used analytic and numerical
0 any pinning sites in the material, the vortexlines form a methods to analyze the zero-temperature structure of
/ triangular lattice known as the Abrikosov lattice2.
t vortex arrays in the presence of periodic pinning. We
a Equilibrium and transport properties of the mixed have also carried out a numerical study of the zero-
m
phaseoftype-IIsuperconductorsarestronglyaffectedby temperature equilibrium magnetization of a supercon-
- thepresenceofpinningcenters,eitherintrinsictothesys-
ducting film with a square array of pinning centers as
d
tem or artificially generated. Understanding the effects
n a function of the applied field. In section II, we consider
of pinning in these systems is very important for practi-
o the ground states of a vortex system in a square array
c calapplicationsbecause the presenceofpinning strongly of pinning centers for fillings less than unity. We look
: influences the value of the critical current in the mixed
v ata classof structures that areBravaislattices with one
i phase. vortex per basis if the filling n is of the form 1/q, and
X Inrecentyears,avarietyofnanofabricationtechniques
with p vortices per basis if n=p/q (p and q are integers
r have been used to create periodic arrays of pinning cen- greater than unity, with p < q). The structure with the
a ters in thin-film superconductors 3–13. Such arrays may
lowest energy in this class can be obtained rather eas-
consist of microholes (“antidots”)3–7, defects produced ily. We find that the “ground-state” structure obtained
bythebombardmentofion8 orelectron9 beams,ormag- this way matches those obtained from experiments8 and
netic dots10–13. These pinning centers are “strong” in simulations15,18 for a large number of simple fractional
the sense that each pinning site can trap one or more
values of n. The results obtained in this section can
vorticesatlowtemperatures. Theeffectsofperiodicpin- also be used to predict the ground-state structures for
ningdependstronglyonthe relativevaluesofB andB,
φ 1 < n < 2. In section III, we consider the ground-state
where Bφ =ρpΦ0 (ρp is the areal density of the pinning structures for fillings greater than two. In these calcu-
centers) is the so-called “matching field”, and B is the
1
lations, we use simple geometric arguments to arrive at Consider now fillings of the form n = 1/q, q being
the ground states. This analysis is performed under the an integer greater than one. Let us look at Bravais lat-
assumptionthatthe rangeofthe intervortexinteraction, tices that can be formed for a specific n by distributing
which is set by the penetration depth, is much smaller the vortices on the square pin lattice with one vortex
than the spacing between the pinning sites. We show per basis. The motivationfor consideringsuchlattices is
that the ground-state structures obtained from this sim- thatthis will automaticallyensure that there is no shear
ple analysis match the ones obtainedfrom simulated an- of the vortices with respect to the pin lattice, since the
nealing. This analysis is extended to rectangular and forces on a vortex due to other vortices will add up ex-
triangular pin lattices in section IV. In section V, the actly to zero. The unit cell area of these structures has
zero-temperature equilibrium magnetization of a vortex to be q.d2. So the possible unit cells can be obtained by
system in a square array of pinning centers is obtained factorizing q into products of the form r.s (r and s are
byfirstcalculatingtheground-stateenergyasafunction integers),arrangingthe vorticesatthe cornersofrectan-
of the magnetic induction and then finding the applied gles of dimension rd sd, and then sliding the parallel
×
field from a numerical differentiation of the data. The sides relative to each other. This procedure produces a
ground-stateenergiesfordifferentvaluesofthe magnetic large number of structures depending on the value of n
induction are obtained using a simulated annealing pro- and we have to pick the one that minimizes the energy.
cedure. The calculated B H curve exhibits interesting For small values of q, this can be done by hand, but as
−
commensurability effects, manifested as plateaus occur- q becomes large and highly factorizable, the number of
ringatsimplerationalvaluesofthefillingfractionn. The possible structures increases rapidly. In such cases, we
main results of our study are summarized in section VI. have resorted to computers to generate these structures
and compare their energies.
The structures so obtained for fillings 1/2 and 1/4
II. GROUND STATES FOR A SQUARE PIN match those found in the imaging experiment8. Also
ARRAY WITH FILLING RATIO LESS THAN for fillings 1/2,1/3,1/4,1/5,1/8,1/9,1/10and 1/15, we
ONE find the same structures as those obtained by solving
the “greedy lattice gas model”23 exactly. This is un-
We consider a superconducting film that has a square derstandable because when the intervortex interaction
array of pinning sites with lattice constant d. The mag- falls off rapidly as the distance is increased above the
netic field is assumed to be perpendicular to the sur- defect spacing d, the ground state can be attained by
face of the film. The “matching field” BΦ is then given finding the lattice that maximizes the shortest distance
by Bφ Φ0/d2, and the filling fraction n is given by between vortex pairs. If two structures have same value
n=B/B≡φ =Bd2/Φ0whereBisthemagneticinduction. andnumber of shortestdistances,then the next shortest
We assume that the pinning potential is much stronger distance should be maximized, and so on. For fillings
thantheintervortexinteraction,butisofextremelyshort 1/2,1/3,1/4, and 1/5, our analysis also yields the same
range. Thelargestrengthofthepinningpotentialimplies structures as those found in the large Ue (Ue is the en-
thatthevorticesmustoccupypinningsitesaslongasthe ergy of on-site repulsion between two electrons) limit of
numberofvorticesdoesnotexceedthenumberofpinning the neutral Falicov-Kimball model24. In Fig. 1(a,b,c),
sites. We alsoassumethatapinning site cannotaccom- we have shown the structures so obtained for a few fill-
modate more than one vortex. If the pinning centers in ings of the form n = 1/q. The ground-state structure
the film are microholes, then this assumption amounts shown in Fig. 1(b) for n = 1/5 is different from that
to the requirement that the radius of each hole is close found in Ref. 16 from a simulated annealing calculation.
to two times the coherence length ξ of the superconduc- This difference is probably due to the use of a different
tor. These assumptions ensure that interstitial vortices (logarithmic) intervortex potential in Ref. 16.
appear only when the filling fraction n is greater than Let us now compare the energies of the nearest and
unity. The assumed short range of the pinning potential next-nearest neighbors in one of these lattices. The in-
can be justified if the defect diameter is small compared teractionenergybetweentwovorticesseparatedbyadis-
tothedefectspacingd. Anotherassumptionthatwewill tance r is given by the expression,
make in most of our calculations is that the intervortex Φ2t r
0
interactionfallsoffrapidlywithdistance. Thisisensured U(r)= 8π2λ2K0 λ (1)
if the penetration depth λ is much smaller that the pin-
(cid:16) (cid:17)
lattice spacing d. In our calculations, we take the ratio where K0 is the zeroth order Hankel function of imagi-
λ/dtobe10. Thisvalueisappropriateforthepinlattice nary argument and t is the film thickness. For n = 1/2,
ofRef.8. Weconsidertemperaturesthatarelowenough thenearestneighbordistanceis√2dandthenext-nearest
toneglecteffectsofdepinningandvortex-latticemelting. neighbor distance is 2d. So the interaction energies are,
Theproblemoffindingthestructureofthevortexsystem to within a constant prefactor, given by
thenreducestolocatingthegroundstateinthepresence
√2d
of the pinning potential. Un K0 =0.2 10−6,
∝ λ ! ×
2
Unn K0 2d =0.6 10−10. in a basis. We have shown in Fig. 1(d,e,f) some of the
∝ λ × ground-state structures obtained this way. These struc-
(cid:18) (cid:19)
tures match those obtained from our simulated anneal-
One can see here that there is orders of magnitude dif- ing calculation. These ground states show the “stripe”
ference in these energies which cannot be compensated structure predicted by Watson23 and Kennedy24 in ap-
bydifferencesarisingfrominteractionswithmoredistant propriate density ranges.
neighbors. This differenceisgoingto bemoreprominent
atlowerdensities. Thistells usthatthe maximizationof
the shortest intervortex distance in a lattice for a given III. GROUND STATES FOR A SQUARE PIN
ARRAY WITH FILLING RATIO GREATER
(a) (d) THAN TWO
If n is greater that 1 but smaller than 2, then the
ground-state structures are similar to the ones for the
case of n less than one. The only difference is that
the pinning sites are all occupied and the centers of the
squareunitcellsofthepinlatticenowactasnewpinning
centers for interstitial vortices. But things look different
(b) (e)
whenthefillinggoesabovethevaluen=2. Forsuchval-
uesofn,wecannolongerplacetheinterstitialvorticesat
the centers of the squares and look for simple structures
obtained this way. Also, we now have to start looking
into the stability of the structures since the square sym-
metry would not be present.
(c) (f)
A. The ground state for n=5/2
Here we are faced with the task of placing more than
one vortex in a square. Before going to the problem of
finding the ground state for n = 5/2, let us ask a more
basic question: given a single unit cell of the square pin
lattice with each corner occupied by a vortex. how can
FIG.1. The ground state structures for a few filling frac- we arrange two more vortices inside this square so as
tions n < 1. The different filling fractions are (a) 1/2, (b) to minimize the energy. Since “greedy lattice gas” has
1/5, (c) 1/9, (d) 2/3, (e) 3/5, and (f) 2/7. The dots in the beenagoodapproximationforthepreviouscases,wetry
figures represent the pinning sites and the circles represent to tackle this problem by “maximizing the shortest dis-
thevortices. tance”method. Inordertostabilizeaninterstitialvortex
bymaximizingtheshortestdistance,itsdistancefromat
filling would lead to the ground states, provided the lat-
leastthreenearestvorticesmustbetheshortestdistance.
tice spacing is large compared to the penetration depth
It is also requiredthat these nearest set of vortices must
of the film. This, in fact, is exactly the definition of the
be spread in such a way that if we draw straight lines
greedylatticegas. However,onehastobecautiousabout
from the vortex in question to these neighbors, the an-
this method because, as noted in Ref. 23, the structures
glesformedbyadjacentlinesmustbelessthan180o. The
can be strongly dependent on the form of the potential
proof of this statement is given in Appendix A.
in certain ranges of n and we can even have aperiodic
It can be seen fromthe symmetry of the problem that
structuresasgroundstates. The ground-statestructures
we have to place the two vortices on the lines joining
shown here have been cross checked with simulations to
the centers of the sides to meet the condition mentioned
ensure that they are indeed the lowest energy configura-
above. This leaves us with only two possible ways of
tions. Togiveanexampleofacasewhere this treatment
doing it, which are shown in Fig. 2. In (a) the shortest
doesnotleadtothe truegroundstate,we foundthatfor
distance D can be obtained by solving the equation
ma
filling 1/16, the energy per vortex for the structure with
minimumenergyobtainedthiswaywasgreaterthanthat 2
1 d
2 2
for filling 1/15, implying that the structure obtained for 4(d−Dma) + 4 =Dma. (2)
n=1/16 was not the ground state.
When the filling fraction is of the form p/q with p not On solving this equation, we get AB = BC = BE =
equaltoone,onecanlookforgroundstatesinasubsetof D = (√7 1)d/3. In (b) the vortices A′,B′ and C′
ma
structureswheretheunitcellhassizeqd2 withpvortices −
3
A D A’ D’
θ m n
C’ p
φ
θ
B
t q
C B’
φ
(a) (b)
E E’
(a) (b)
FIG.2. Puttingtwovorticesinasquare. Theconfiguration
in (a) is the global energy minimum, and the configuration
in (b) can at best be a local minimum of the energy. The
angles θ=15o and φ=24.49o. The squareis drawn for easy (c) (d)
visualizationandthedottedlinesarethebisectorsofthesides.
FIG.3. Possible 2d 2d unit cells for n = 5/2. The unit
Note that there is already one pinned vortex at each of the
×
cell (a) is the lowest energy configuration for d/λ = 10. For
cornersofthesquare. ThedistancesAB=BC =BE =Dma
and A′B′ =B′C′ =B′E′ =D . muchlargervaluesofd/λ,theunitcellwiththelowestenergy
mb
will bethe onein (b).
form an equilateral triangle. Thus the nearest neighbor
distance A′B′ =B′C′ =B′E′ is between vortices m and n) compared to (a). But this
is done at the cost of bringing in interactions like those
d
D =sec(15o) (3) betweenvortexpairs(p,q)and(p,t)forevery“gain”ofa
mb
2 next-nearest neighbor interaction. The energies of these
twointeractionsford/λ=10arefoundtobequiteclose.
The angle φ in (a) is 24.49o, and the angle θ in (b) is
These energies are
15o. The interactionenergiescorrespondingto these two
distances for d/λ=10 are Umn K0(rmn/λ) 3.2 10−4,
∝ ≃ ×
√7 1
UAB ∝K0 3−λ d!≃2.2×10−3 Upq ∝K0(rpq/λ)≃1.9×10−4,
where r and r are the distances between vortices m
mn pq
sec(15o) andn,andpandq inFig.3,respectively. Itis clearfrom
UA′B′ K0 d 3.0 10−3. this comparisonthat the unit cell (a)wouldbe preferred
∝ 2λ ≃ ×
(cid:18) (cid:19) for d/λ=10.
From comparing these two energies it is clear that (a) The ratio of interaction energies of the next-nearest
is the global minimum, whereas configuration(b) can at and the nearest neighbors is 0.07 for this lattice when
best be a local minimum. d/λ = 10. This energy difference is appreciable here, so
Comingbacktothen=5/2case,wenowhavetobuild that we can expect that the unit cell we arrived at is
up the lattice with equalnumber of two types of squares the correct one. Note that any net force that might be
- one with two interstitial vortices and the other with presentononeoftheinterstitialvorticesduetotheasym-
one interstitial vortex. Note that here we have neglected metryinthestructurecanbecompensatedbyextremely
structures that have three or more interstitial vortices smalldisplacementsfromthepositionsobtainedfromthe
inside a square unit cell because such structures would “maximization of the shortest distance” method.
drasticallybringdownthenearestneighbordistance. Let
usnowlookatthepossibleunitscellsofsize2d 2dthat
canbemadeoutofthesetwotypesofsquares. T×heseare B. The ground state for n=3
shown in Fig. 3. If one constructs the lattice with these
unit cells, configuration (b) offers the least number of When the filling fraction equals 3, we have to build
next-nearest neighbors, the number of nearest neighbors up the lattice using blocks of type (a) in Fig 2. Again,
being the same inall the cases. So one canexpect (b) to looking at unit cells of size 2d 2d or smaller, we have
×
be the ground state unit cell, at least for large values of the configurationsshowninFig. 4to consider. Hereitis
d/λ. However,unit cell(a)is preferredif d/λis not very easy to see that unit cell (b) is preferredoverthe others.
large. This can be understood in the following way: the Thisisbecauseitistheconfigurationthatmaximizesthe
advantagethat(b)hasover(a)isthatithasonlyhalfthe minimum distance between any two vortices in different
number of next-nearest neighbors (interactions like that squares,thedistancesbetweenvorticeswithinonesquare
4
A’ A
q
p θ θ
C
C’
s ϕ
θ θ ϕ
ϕ
B’ D’ θ B
D
(a) (b)
E’ E
(a) (b)
FIG.5. Two possible ways of arranging three vorticesin a
squareso as to maximize the shortest distance. The configu-
rationin(a)offersthebestarrangementiflookedinisolation,
(c) (d) but the configuration in (b) wins out when one has to con-
struct a lattice of the unit cells. The angles θ = 15o and
FIG. 4. Possible 2d×2d unit cells for n = 3. Configura- ϕ=60o.
tion(b)canbeeasilyseentoofferthemaximumnext-nearest
neighbor distance, the nearest neighbor distances being the
same and equalin numberin each case. Considerthe figureshowninFig. 6. The solutionweare
looking for can be obtained by solving the equations
being the same in all the configurations. Again compar-
AB =BC =CA=D , (4)
s
ing the nearest interaction and the next-nearest one, we
have
CB′ =AC′′ =AP =D . (5)
Un K0(rsp/λ) 2.2 10−3, s
∝ ≃ ×
Note that here we have assumed a unit cell size d d.
×
Onsolvingtheseequations,weobtaintheunitcellshown
Unn ∝K0(rpq/λ)≃1.9×10−4. in Fig. 5(b). In the figure, the angle θ equals 15o and
the angle ϕ is 60o. In this lattice, the shortest distance
There is appreciable difference between these two val-
is D = sec(15o)d/2 and the next shortest distance is
ues, and hence, the ground state we have arrived at is s
D =3dsec(15o)/(4√2).
reasonable. When the filling lies between 2 and 3, one ns
This simple solution may not be the correctone if the
cansafelyassumethatthe groundstatestructurecanbe
latticespacingisnotlargeenough. Forexample,ifd/λ=
built up using squares of type (a) in Fig. 2, and squares
10, as we have been assuming when comparing energies,
that have one vortex at the center. In fact we make use
then the nearest and next-nearest neighbor interaction
of this in our simulations to arrive at the ground states,
energies turn out to be really close. Hence we can not
as described in section V.
rule out the possibility of the lattice arranging in such a
waythattheshortestdistanceisreducedsoastodecrease
the number of nearest or next-nearest neighbors. The
C. The ground state for n=4
relevant energies for d/λ=10 are
Here we have to place three interstitial vortices in one
D
square. Thisisnontrivialsinceevenifweensurethatthe K0 s 3.1 10−3,
λ ≃ ×
shortestdistanceismaximizedinonesquare,twovortices (cid:18) (cid:19)
in nearby squares may be closer to each other than the K0 Dns 2.1 10−3.
shortestdistancewithinasquarewhenwecreatethelat- λ ≃ ×
(cid:18) (cid:19)
tice. Note that we did not come across this problem in
the n=5/2 or n=3 fillings. To illustrate this problem, So the ground state obtained above is guaranteed to be
we show inFig. 5(a)a single cell ofa squarelattice with the correct one only for much larger values of d/λ. The
three interstitial vortices, which may very well be a lo- ground states we have obtained for n = 5/2 and n = 3
cal minimum configuration. But if we try to build the match with the images from experiments8 and also the
lattice using this cell, we can not do it without bringing results of simulations15. But for n = 4, the structures
the vortices in nearby squares closer than the minimum found in experiments and simulations are different from
distance in an individual cell. So what we need to do is theoneshowninFigs5(b)and6. Thisisexpected,since
to look for a pattern that will include vortices in differ- in the experiment the value of the ratio d/λ was close
ent squares while doing the minimization of the shortest to ten. Our simulation with d/λ= 10 gives the ground-
distance. We can solve this problem trigonometrically. statestructureshowinFig. 7(a),whichissimilartothe
5
the longer side of the rectangle. We shall consider here
A’’ only the cases where the filling is greater than one. In
the absence of the square symmetry, it is obvious that
B’’ the ground state structure will depend not only on the
C’’
penetrationdepthλ,butalsoontheratiol/b. Inthefol-
lowing analysis, we shall always assume that λ is much
R
P smaller than b, the shorter side of the basic rectangular
A A’
pinning cell. When the filling is 2, for values of l/b less
than √3, the ground state is one where each interstitial
B B’ vortexisatthecenteroftherectangle,sincethisensures
C C’
that the shortest distance is maximized (see Fig. 8(b)).
Butwhentheaspectratioexceeds√3andtheinterstitial
Q
C R
FIG.6. Thedistances onehastoequatewhenmaximizing
the shortest distance for n = 4, taking into consideration A
distances between vortices in neighboring squares. Note that
theunit cell size here is one square unit. P
D F
oneobtainedintheexperiment8. Inoursimulationswith B
verylargevaluesofd/λ,wegetstructuressimilartothat E Q
in Fig. 6. The simulation result for d/λ = 50 is shown (a) (b)
in Fig. 7 (b), which matches well with the predicted
structure. It is to be noted that the simulation result FIG. 8. The unit cells of the ground state, obtained
was obtained by starting the system near the expected by maximizing the shortest distance, when the filling is
ground state. So the claim is that it offers at least a 2 for a rectangular array of pinning sites. (a) The unit
localminimumoftheinteractionenergy. Thesimulations cell when the aspect ratio is greater that √3. The inter-
were carried out for different system sizes form 2d 2d stitial vortex is displaced horizontally from the center of
to10d 10dtoruleoutanydependenceofthe result×son the rectangle by distance Ds. In the figure, the distances
× AB = AC = AD = BF = BE. (b) The ground state unit
the boundary condition.
cell when theaspect ratio is less that √3, where theintersti-
tial vortex is at the centerof therectangle.
40 200
vortices are placed at the centers of the rectangles, the
30 150
distance between two interstitial vortices in neighboring
20 100 cells would be shorter than that between an interstitial
vortexandthe closestpinnedvortex. This wouldleadto
10 50 a displacement of the interstitial vortices sideways from
the center, along the bisector of the shorter sides of the
0 0 rectangle,tomaximizetheshortestdistance. Theresult-
0 10 20 30 40 0 50 100 150 200
ing structure is shown in Fig. 8 (a). The displacement
(a) (b)
of the vortex from the center is given by
FIG. 7. The ground state for n = 4, obtained by simu- l+√4l2 9b2
lated annealing when the ratio of the penetration depth to Ds = − − . (6)
6
the pin-lattice spacing is (a) 10 and (b) 50. The dark dots
denotethepinningsitesandtheaxislabelsareinunitsofthe It is worth noting that since the vortex in the center
penetration depth. wouldbemovingtowardstwoofthepinnedvortices,and
away from only one interstitial vortex per unit cell, the
displacement will approach the value given above only
when the ratio l/b is appreciably large and in the limit
IV. GROUND STATES FOR RECTANGULAR of small penetration depth compared to the sides of the
AND TRIANGULAR PIN ARRAYS rectangle. Forexamplewe havefoundinoursimulations
that even when the ratio l/b is 2, the ground state for
One can extend this type ofanalysis to pinning arrays b/λ=15is onewhere the interstitialvortexis veryclose
with other symmetries for finding the least energystruc- to the center, whereas for the l/b=3 and b/λ=15, the
tures for simple filling fractions. Let us firstconsiderthe ground-state structure is quite close to the one obtained
case of a rectangular array of pinning sites with a pin- frommaximizingtheshortestdistance. Someofoursim-
ning unit cell of dimensions l b, where we take l to be ulation results for filling equal to two are shown in Fig.
×
6
12 (a,b). Also, if the ratio l/b becomes too large, the l/b = 5/4. The unit cells that provide the largest min-
analysis will have to include more than two of the inter- imum distance are shown in Fig. 10. Note that when
stitial vortices,since now the solution like that shown in l/b = 4, the vortices are arranged parallel to the longer
Fig. 8 (a) can lead to two vortices being closer in the
next nearest cells or ones even further apart. C D
In trying to arrive at the lowest energy structures for A
fillings 5/2 and 3, it is important to determine how one B
E
can accommodate two vortices in a rectangularcell with
F
the shortest distance being maximized. There are two
possible minima that one has to consider: one where the G
(a)
vortices are arranged along the line dividing the shorter
sides, and one where they are arranged along the line
U V
dividing the longer sides, as shown in Figs 9 (a) and 9
(b), respectively. The shortest distance in each case is P
given by,
Q
√4l2+3b2 l W S
Ds1 = − , (7)
3
√4b2+3l2 b R
Ds2 = − . (8)
3 T
(b)
If one considers the distances within the cell, configura-
FIG. 10. The unit cells for n = 5/2 for a rectangu-
B E B E
lar array of pinning sites. (a) The unit cell when the as-
pect ratio is 4. The distances AB = AC = AE = BD.
A’
(b) The unit cell when the aspect ratio is 5/4. Here
A
PV =PU =PQ=QW =QS =Ds2 and QR=RS =RT.
D
When this structure is repeated periodically, the “image” of
D’ the vortex at P would be at the same distance QR from R.
This would ensure thestability of thevortex at R.
C C
(a) (b) side and in the other case, parallel to the shorter side.
Also when l/b = 5/4, the vortex in the single-vortex
FIG. 9. The two possible ways by which the shortest dis- rectangular cell is not located at the center, but slightly
tance can be maximized when there are two vortices in a displacedsidewaysalongthe bisectorofthe shortersides
rectangular cell. (a) In this case the vortices are placed on
to facilitate the maximizationof the second shortestdis-
the line that divides the shorter sides of the rectangle. The
tance involved. One should againkeep in mind that this
distances AB = AC = AD = DE = Ds1. (b) Here the vor- sort of analysis will go wrong if the aspect ratio is too
ticesareplacedparalleltotheshortersidesandthedistances
large, since then the distances between vortices in next-
A′B = A′E = A′D′ = CD′ = Ds2. For l/b less than 2, the nearest or further neighbor cells will become important.
configurationin(b)leadstoalargershortestdistance(within
In Fig. 12 (c) we have shown the ground state struc-
thecell)thantheonein(a). Ifl/bisgreaterthan2,thenthe
ture obtained from simulations for aspect ratio 4, with
distance Ds2 becomes greater than b and the vortices spill
b/λ = 5 and in Fig. 12 (d), the ground state for aspect
overinto thenext cell.
ratio 5/4, with b/λ= 24. In Fig. 12 (d), one can notice
tion (b) in Fig. 9 gives the lowest energy. But for large the slightshiftofthevortexfromthe centerinrectangu-
values of the ratio l/b, this configuration is disfavored lar units having a single interstitial vortex. The shift is
since it allows the vortices in one rectangle to get close notaslargeas expectedfromthe aboveanalysis. This is
to those in a neighboring one. Also, for l/b > 2, the due to the large but finite value of b/λ.
interstitial vortices “spill over” into the next cell, since For filling equal to 3, the lowest energy structures ob-
the distance Ds2 becomes greater than b. So, one has tained by considering 2×2 unit cells for two values of
to work out the structures for fillings 5/2 and 3 case by the aspect ratio, 2 and 5/4, are shown in Fig. 11. Here
case. For fillings between 2 and 3, one has to choose too, for large values of the aspect ratio, the structure is
appropriatenumber of two-vortexrectanglesofthe right composed of rectangular cells where the interstitial vor-
kind and single-vortex ones and arrange them so as to tices are aligned parallel to the longer side (Fig. 11(a)),
maximize the shortest distance appearing in the struc- whereas when the aspect ratio is smaller, the structure
ture. We have looked at the 2 2 unit cells possible for is made up of an alternating arrangementof rectangular
filling 5/2 for two values of the×aspect ratio, l/b=4 and cellsofbothtypes(Fig. 11(b)). Thesimulatedannealing
resultsforsimilarvaluesoftheaspectratio(Fig. 12(e,f))
7
C D
B
A
45 45
E F
30 30
15 15
(a) 0 0
0 45 90 135 0 20 40 60
S T X (a) (b)
20 96
P
V 15 72
Q W 10 48
Y 5 24
R U
0 0
0 20 40 60 80 0 30 60 90 120
(c) (d)
40 64
(b) 30 48
20 32
FIG. 11. Unit cells for n = 3 for a rectangular array of 10 16
pins. (a) The unit cell when the aspect ratio is 2. The 0 0
0 20 40 60 80 0 20 40 60 80
distances AB = AC = AE = BD = BF. (b) The unit
(e) (f)
cell for the same filling but for aspect ratio 5/4. Here the
distances PQ = PS = PT = QR = QU = Ds2 and
FIG. 12. The ground state structures obtained by simu-
TV =UV =VW =WY =WX =Ds1.
lated annealingwhenthepinningarrayhasrectangular sym-
metry. The parameters for different plots are, (a) n = 2,
yield the structures obtained from the above analysis. l/b = 3 and b/λ = 15, (b) n = 2, l/b = 4/3 and b/λ = 15,
Groundstatestructuresobtainedbysimulatedanneal- (c) n=2.5, l/b=4 and b/λ=5, (d) n=2.5, l/b=5/4 and
ing for a rectangular pin array with l/b=2 and integral b/λ = 24, (e) n = 3, l/b = 2 and b/λ = 10 and (f) n = 3,
values of n are reported in Ref. 17. In that study, the l/b=5/4andb/λ=16. Thedotsdenotethepinningcenters
andthecirclesrepresentthevortices. Notethattheboxsizes
intervortex interaction was assumed to depend logarith-
are not to scale.
mically on the intervortex distance. The ground state
structure found in Ref. 17 for n = 2 is similar to that
shown in Fig. 12(b), but the structure found there for
n=3isquitedifferentfromthoseofFig. 12(e,f). Thisis
anotherexampleoftheimportanceofthedetailednature
of the intervortex interaction in determining the struc-
ture of the ground state.
For a triangulararrayof pinning sites, it is easy to see
B R S
that when the filling is greater than 1 and less than 3,
the interstitial vorticeswill be placedon the centroidsof
the triangles in the limit where one can safely apply the Q
method ofmaximization of the shortestdistance. So the A
ground states when n is between 1 and 3 will be made P
C
up of parallelogram cells of the form shown in Fig. 13. D
These unit cells match wellwith the results of molecular (a) T (b) U
dynamics simulations18,15.
FIG.13. Basicbuildingblocksforgeneratinggroundstates
fortriangularpinningarrayswhenthefillingisbetween1and
V. EQUILIBRIUM MAGNETIZATION OF THE 3. (a) When there is a single interstitial vortex in a parallel-
GROUND STATES ogram, the shortest distance can be maximized by placing it
atoneofthetwocentroidsofthetrianglesinvolved. Thedis-
In this section, we describe a calculation of the zero- tancesAB=AD=AC =a/√3,whereaisthelengthofthe
temperature,equilibriummagnetizationofathin-filmsu- sideofapinningcell. (b)Whentwovorticesaretobeplaced
perconductorinthepresenceofasquarearrayofpinning in a single pin cell, they haveto be at the two centroids.
sites. The regionin the B H plain we are interestedin
−
is that just above Hc1, when the flux tubes start enter-
ing the sample. The idea is to find the free energy F of
8
the groundstate as a function ofthe magnetic induction 40 40
B, and then to obtain the applied magnetic field H by
30 30
taking a derivative of the free energy with respect to B.
20 20
Since we are considering the zero-temperature case, the
free energy is just the internal energy of the flux lattice. 10 10
Since we are looking for a nearly continuous variationof
0 0
the internal energy for taking the derivative, we need to 0 10 20 30 40 0 10 20 30 40
(a) (b)
locate the groundstates for filling fractionsseparatedby
40 40
smallintervals. Thiswouldbedifficulttodoanalytically,
30 30
since the unit cells for some filling fractions can be arbi-
trarilylarge. Also,as nbecomes large,the simple proce- 20 20
dureofmaximizationoftheshortestdistanceisnotgoing
10 10
to yield the correct ground-state structures. So we have
resorted to simulations to determine the ground states. 00 10 20 30 40 00 10 20 30 40
Inparticular,wehaveusedthesimulatedannealingtech- (c) (d)
nique to locate the global minima (or at least low-lying
local minima close in energy to the global ones) of the FIG.14. The ground states for (a) n = 9/4, (b) n = 5/2,
part of the internal energy associated with intervortex (c) n = 11/4, (d) n = 3 obtained by simulated annealing,
interactions. as discussed in the text. The unit-cell size is 4d 4d. The
×
The Helmholtz free energy per unit volume of the su- dark dots denote the pinningsites and thecircles denotethe
perconductor at zero temperature in the presence of the vortices. Theaxislabelsareinunitsofthepenetrationdepth.
pinning sites is
nǫ nǫ where N is the number of basic pinning squares in the
l p
F (n)= +E (9)
s d2 n− d2 systemandrij isthe separationbetweenvorticesi andj
in the ground state for the filling fraction n.
wherethefirsttermisthelineenergy,the secondtermis Thesizeofthesystemswesimulatedvariesfrom2d 2d
the interaction energy and the third term is the pinning to 8d 8d. In all cases, we used periodic boundary c×on-
energy. Here,ǫlisthelineenergyperunitlength,ǫpisthe dition×s to minimize surface effects. So the minimum
pinning energy per unit length and En is the interaction difference between two consecutive filling fractions was
energy per unit volume for filling fraction n. We note ∆n = 1/64. The ratio d/λ was taken to be 10, as
herethatthepinningenergyincreaseslinearlywithntill in our previous analysis. In order to save computation
nbecomesoneandthenremainsconstant,sincemultiple time, the vortices were allowed to stay only at the pin-
occupation of a pinning center is not allowed. Further, ning sites when the filling was less than one. For fill-
for simplicity, we express the pinning energy as ings between one and two, every pinning site was occu-
pied by a vortex which was never moved and the extra
ǫ =mǫ (10)
p l ones were allowed to move near the centers of the ba-
sic pinning squares. When the filling was greater that
where m is a positive number whose magnitude depends
two and less than three, the vortex configurations were
onthenatureofpinning. Theinteractionpartofthefree
constructed using basic units of squares containing one
energy, E , is the computational input. Once we know
n vortex at its center and squares containing two vortices
the free energy, we can compute the applied magnetic
placed such that the shortest distance within a square
field H using the relation
is maximized (as in Fig. 2(a)). These units were then
∂F H moved around and twisted while cooling to arrive at the
s
= . (11) minimum energy states. This procedure helped us to
∂B 4π
tracklow-lyingminima faster thanif we allowedvortices
Usingthe standardexpression1 forǫ andtakingthe log- to move freely. Once the basic structure was thus ob-
l
arithmoftheGinzburg-Landaunumbertobeequalto2, tained, the vorticeswere allowedto move freely during a
we get the following expression for the applied field as a second cooling schedule starting from a lower tempera-
function of the filling fraction: turetoobtainthe lowest-energystructure. InFig. 14we
have shown some of the ground state structures we have
H = Φ0 [1 m 1 Θ(n 1) ]+ ∂En′ . (12) obtained this way for fillings between 2 and 3. For fill-
2πλ2 − { − − } ∂n ings 5/2 and 3, we find that the structures match those
obtainedinexperiments8,aswellasinouranalysisusing
Here E′ is given by the expression
n maximization of the shortest distance. The structures
for n=9/4 and n=11/4 may not be the actual ground
En′ = 2πΦλ02N K0(rij/λ), (13) states, either due to the smallness of the unit cell of our
i>j simulation or due to the presence of many nearly degen-
X
9
erate local minima. From the energy versus filling fraction data, one can
findtheappliedfieldusingEq.(11)andthencomputethe
0.8 magnetization M using the relation
B =H +4πM (14)
0.6
With Substrate In Fig. 16 we have plotted B versus H in the entire
E range of filling for which simulations were carried out,
n
0.4 from n = 0 to n = 3. Figs 17 and 18 show magnified
versions of this plot in the regions between filling frac-
Homogeneous
tions 0 and 1, and between 2 and 3, respectively. Note
0.2
that we have not explicitly included the pinning energy
term in our analysis. This term would just add a con-
stant contributionto H for fillings up to 1. The features
0
0 1 2 3 of the curve from n = 1 to n = 2 are the same as those
n in the interval between n = 0 to n = 1. This is due to
the fact that the ground state structures are similar in
FIG. 15. The energy of the vortex lattice in the presence the two regions (see section III). The B H plot shows
−
of a square array of pinning sites (upper dashed curve). The flat regions at values of B corresponding to fillings 1/8,
energy of the triangular vortex lattice in the absence of pin- 1/5, 1/4, 1/2, 3/4, 4/5, 7/8, 1, 9/8, 6/5, 5/4, 3/2, 7/4,
ning is also shown for comparison (lower solid curve). The 9/5,2inthefillingfractionrangebetween0and2. Also,
system size is 8d 8d. The abscissa is the filling fraction n in the rangeof n between2 and3, there areroughly two
×
andtheordinateisthetotalenergyperunitthicknessinunits plateaus, appearing near n=2.3 and n=2.6.
of Φ20/(8π2λ2).
1
In Fig. 15 we have plotted the ground state energies
obtained from the simulation for different fillings. The
simulation unit cell was 8d 8d and the energies were 0.75
×
computed for fillings 1/64 to 3. The upper curve shows
φ
the results obtained in the presence of the pinning sites B
B/ 0.5
andthelowercurveistheenergyofthetriangularlattice
for the same density of vortices. Note that we have not
includedthepinningenergyintheplot. Thiswouldbring 0.25
down the upper curve below the curve for the pin-free
case.
15.7562 15.7569 15.7577
3
H/B
φ
2.5
FIG.17. ExpandedviewoftheplotofFig.16intheregion
2 between n = B/B = 0 and n = 1. One can see plateaus
φ
Bφ appearing at n= 1/4, n =1/2 and n = 3/4. The dark dots
B/ 1.5 arethedatapoints. Thelightdottedlineisdrawnasaguide
to theeye.
1
0.5 The observed values of the filling fractions between
n = 0 and n = 1 at which the plateaus occur indicate
that these values of n correspond to fillings where the
15.75 15.85 15.95 16.05
introductionofanewvortexintothesystemleadstothe
H/B
φ appearance of a shorter distance than those existing in
the lattice or in its dual (i.e. the lattice obtained by re-
FIG.16. The dependence of the magnetic induction B on placing particles by holes, andvice versa)23. This makes
the applied field H in the entire region of our simulation. A sense because the introduction of this shorter distance
number of plateaus can be seen at points corresponding to brings in a larger energy scale, leading to a discontinu-
simple rational filling fractions mentioned in the text. Only ous change in the derivative of the energy with respect
theverticallinescorrespondtothedatapointsobtainedfrom to the filling n. In the simulations, we have not scanned
our calculation; the dotted lines are guides to the eye. Both very small intervals of n. Also, for some of the fillings,
B and H havebeen scaled by thematching field Bφ. thegroundstatemaynothavebeenobtainedinoursim-
ulated annealing calculation. For these reasons, we can
10
