Composite vortex ordering in superconducting films with arrays of blind holes

The pinning properties of a superconducting thin film with a square array of blind holes are studied using the nonlinear Ginzburg–Landau theory. Although blind holes provide a weaker pinning potential than holes (also called antidots), several novel vortex structures are predicted for different size and thickness of the blind holes. Orientational dimer and trimer vortex states as well as concentric vortex shells can nucleate in the blind holes. In addition, we predict the stabilization of giant vortices that may be located both in the pinning centers and/or at the interstitial sites, as well as the combination of giant vortices with sets of individual vortices. For large blind holes, local vortex shell structures inside the blind holes may transfer their symmetry to interstitial vortices as well. The subtle interplay of shell formation and traditional Abrikosov vortex lattices inside the blind holes is also studied for different numbers of trapped vortices.


Introduction
The behavior of superconducting vortices in the presence of periodic arrays of antidots reveals a wide range of interesting commensurability and matching effects when the number of vortices is an integer or a simple rational multiple of the number of pinning sites [1,2]. The resulting vortex lattice structures have a unit cell that is a multiple of the underlying antidot lattice. These stable vortex configurations give rise to specific jumps in the magnetization M(H ) and the enhancement of the critical current j c (H ) (see [3] and references therein). In this paper, we consider the more general case-superconducting films with arrays of blind holes in the presence of a uniform applied magnetic field. These not fully perforated holes can trap several single quanta vortices, whereas all vortices trapped by an antidot fuse in a piercing magnetic flux corresponding to a single multiquanta vortex. Furthermore, the interaction of a flux line with an antidot is different from that with a blind hole. In both cases, the normal/superconducting boundary condition which restricts the supercurrents to flow parallel to the walls of the hole leads to an attractive force due to an image antivortex inside the hole [4]. For the antidots, the interaction force acts along the total length of the flux line, whereas for blind holes a weaker force, roughly proportional to the depth of the hole, is expected. In addition, once a flux line is trapped by an antidot, the notion of the vortex core is lost. In contrast to this, vortices captured by blind holes remain as single-quantum units with a well-defined core and are accommodated inside the hole while minimizing their mutual repulsion [5].
Calculations within the linearized Gizburg-Landau (GL) theory show that the critical field of a superconducting film with a blind hole is sensitive to the bottom layer thickness, but the number of vortices which nucleate inside the hole is not [5]. Recent ac susceptibility measurements [6] also show a weaker vortex-pinning potential and smaller saturation number of the blind holes as compared with the antidot case. However, our previous calculation for a superconducting disk with a single blind hole, based on nonlinear GL theory [7], showed that both the critical field and the number of vortices inside the blind hole strongly depend on the radius and bottom thickness of the blind hole and can be equal to the values found for a superconducting ring.
Different vortex states and transitions between them in thin superconductors with periodic pinning arrays in which multiple vortices can be trapped at individual sites were studied in [8] using molecular dynamic calculations. A variety of vortex states, including collective dimer, trimer and composite states with orientational ordering, was found as a function of the pinning strength and the size of the pinning sites. Although the general behavior of the vortex lattice was accurately described in these calculations, the approximations made there were only valid in a certain range of parameters. Namely, those vortices are crudely considered as classical point particles and the pinning was simply introduced through a model potential, usually a simple attractive potential. The situation becomes more complicated if vortices are treated as extended objects: they can overlap and even form a giant vortex [9]- [11]. In addition, due to the finite size of the vortex core and its elastic properties there exists a local repulsive component in the generally attractive pinning force, which depends on the size of the pinning site [12], which has never been considered in molecular dynamics simulations [2,8]. This modulated pinning force may also influence different dynamic vortex phases [13]- [15].
Our superconducting film with arrays of blind holes can be closely related to the system of colloidal particles on an array of optical traps [16]- [20]. For example, two-dimensional colloidal particles exposed to periodic external electromagnetic fields exhibit a variety of molecular 3 crystalline phases, e.g. orientation dimer and trimer states [17], which is also obtained in our numerical calculations. Moreover, the melting process of these colloidal crystals [18] is very similar to the melting of vortex lattices on pinning arrays [15].
Note also that vortex configurations in superconducting films with an array of blind holes can be thought of as an ideal analogue of macroscopic Wigner crystals in a tunable pinning potential landscape, such as the ones recently studied experimentally in [21]. The transitions shown there between the square-pinned and the partially pinned states occur by locally decreasing the electrostatic pinning potential (achieved by nano-engineered electrodes). Contrary to the case of antidots, the pinning strength of blind holes can be tuned by varying the bottom thickness of the blind holes for arbitrary blind hole radius. In the case of antidots, the radius of the holes should be decreased considerably in order to obtain the partially pinned vortex phases [22]. In addition, the tunability of the pinning potential of the blind holes is most beneficial for several modern vortex ratchet realizations [23]- [26].
Recently, Berciu et al [27] have shown that the spin of electrons, having great potential in storing and processing of information in nano-structured devices, can easily be manipulated by the stray fields of superconducting vortices. They have shown that in a hybrid system consisting of a superconducting film and a diluted magnetic semiconductor (DMS), the inhomogeneous magnetic field of a superconducting vortex creates sufficient Zeeman potential landscape in the DMS to localize the charge carriers. In other words, the control of the vortex dynamics in the superconductor translates into controlled manipulation of the spin/charge textures in the DMS. It is of particular interest to induce interaction of adjacent spin textures by, e.g. vortices approaching each other, and one of the suggested mechanisms for that behavior was the case of two vortices in a blind hole under an applied dc current [27]. This strengthens the need for an exact knowledge of vortex configurations in the blind holes as a function of relevant parameters. The imaging of the vortex structures in large blind holes was performed some time ago [5], but no systematic theoretical study has been done up to now.
In the present paper, we study for the first time vortex structures in superconducting thin films with arrays of blind holes in the presence of a uniform applied magnetic field using the full nonlinear GL theory. The paper is organized as follows. The theoretical formulation of the approach is given in section 2. Phase diagrams for the considered parameter space are shown in section 3, with focuses on the first-order phase transitions between states with different blind hole occupation number (i.e. number of captured vortices) and second-order structural phase transitions for the fixed vorticity inside the blind hole. Vortex shell structures inside the blind holes and the process of shell filling for different number of pinned vortices are studied in section 4. Section 5 summarizes our findings.

Theoretical formalism
We consider a superconducting film (of thickness d λ, ξ ) with a regular square array of blind holes (of radius R, period W and thickness d b < d) in the presence of a perpendicular uniform magnetic field (see figure 1). For the given system we solved two nonlinear GL equations averaged over the sample thickness [9], which can be written in dimensionless units in the following form [7]: where is the density of superconducting current and κ * = λ 2 /dξ is the effective GL parameter. The last term in equation (1) describes the effect of the sample thickness variation on the superconducting condensate. The coordinate-dependent thickness takes the values d(x, y) = d b inside the blind hole and d(x, y) = d outside the blind holes. Equations (1), (2) are solved using the iterative procedure from [9] after they were discretized using the link-variable approach [28]. The magnitude of the applied magnetic field H = n 0 /S = H n is determined by the number n of flux quanta 0 = hc/2e = 2.07 × 10 −7 G cm 2 piercing through the simulation region area S. Simulations are done for a 2 × 2 unit cell with the following periodic boundary conditions for the vector potential A and the order parameter around the square simulation region [29]: where b i , i = x, y are lattice vectors, and η i is the gauge potential. The dimensionless Gibbs free energy is calculated as where integration is performed over the primitive cell volume V and A 0 is the vector potential of the applied uniform magnetic field. In this paper, we restrict ourselves to superconducting samples significantly thinner than the magnetic field penetration depth (i.e. κ * = λ 2 /dξ 1).

Vortex phase diagrams at low applied fields
Contrary to the situation for an antidot lattice, where flux is confined in the hole due to the surrounding supercurrents, vortices trapped inside the blind hole can move more freely within the imposed boundaries, interacting with each other repulsively. In what follows, we will consider the case of integer matching fields, i.e. an integer number of flux quanta per blind hole. Figure 2 shows the equilibrium vortex phase diagram for the second matching field as a function of the blind hole radius R and period W , where the solid lines denote first-order transitions between states with different numbers of vortices inside the blind hole n 0 and the dashed ones depict second-order configurational transitions for unchanged vorticity inside the blind hole. The ground-and metastable states are determined in our calculation by comparing the energy of all stable vortex states found when starting from different randomly generated initial conditions. As in the case of an antidot lattice [30], the occupation number (n 0 ) of each blind hole depends not only on the hole radius R, but also on the proximity of the neighboring blind holes in the lattice. For small radii of the blind holes, one vortex is pinned by the blind hole and the second one is located in the interstitial site (inset 1). The pinned vortices are confined and appear more compressed in comparison with interstitial vortices. In the n 0 = 2 state, we first have a giant vortex trapped in the blind holes (inset 2) and with increasing the radius and the period these giant vortices split into individual vortices (note that in bulk superconductors and nonstructured films the energy of a giant vortex state is always larger than the energy of a multivortex state). The transition between these two states is indicated by the black dashed curve. The vortices inside a blind hole repel each other and move to the edges of the pinning sites. The interaction between vortices in the neighboring blind holes give rise to a parallel ordering of the dimers along the 45 • angle (inset 3). By further increasing the blind hole radius, vortices adjust 6 themselves in such a way that the adjacent dimers are rotated 90 • with respect to each other (inset 4). For comparison, the open circles in figure 2 show the transition between the n 0 = 1 and the n 0 = 2 state for the superconducting film (thickness d = 0.1ξ ) with an array of antidots instead of blind holes.
In order to see the influence of the film thickness-to-blind hole thickness ratio on the vortex structure, we calculated the phase diagram also for a thicker film (d = 0.5ξ ), as shown by the thicker curves in figure 2. In this case, the threshold blind hole radius for the n 0 = 2 state decreases compared with the case of a thinner film due to the enhanced screening of the applied field. Due to the latter effect, the giant vortex state (region (2) in figure 2) becomes energetically more favorable and is stable over a larger portion of the phase diagram. For thicker films we did not find the orthogonal dimer state which was found for thinner films (inset 4 of figure 2), indicating clearly the weaker interaction between the vortices in neighboring pinning centers in the case of thicker films. Note also that the effective GL parameter κ * = λ 2 /ξ d decreases (hence vortex-vortex repulsion weakens) with increasing film thickness d.
The equilibrium vortex lattices for H = H 3 are shown in figure 3 as a function of the radius R and period W of the blind holes. In the n 0 = 1 state, one vortex is pinned by the blind hole and for large W interstitial vortex dimers in adjacent cells are rotated over 90 • , preserving the twofold symmetry (inset 2 of figure 3). Note that the vortex lattice unit cell is four times that of the blind hole lattice. This state is identical to the vortex state observed experimentally [1] in a superconducting film with a square array of antidots in the case of the third matching field. With decreasing the period W , the two vortices are strongly caged at the interstitial sites by the blind holes, resulting in the formation of giant vortices (inset 1). For increasing radius of the blind holes, each hole captures two vortices forming a giant vortex in the blind hole. With further increasing the blind hole radius and the period these giant vortices split into multivortices with different possible orientations (insets 4 and 5), as was found for H = H 2 , where no interstitial vortices were present (cf insets 3 and 4 of figure 2). In the n 0 = 3 state, pinned vortices form first a giant vortex (inset 6) and further split into individual vortices (insets 7 and 8). Vortices show an orientational phase transition with changing lattice parameters. Although some of the vortex configurations (insets 2 and 8 in figure 3) for H = H 3 were obtained previously using molecular dynamic calculations [8], we found here new configurations with different orientation of vortex dimers (insets 4 and 5 in figure 3) and trimers (inset 7 in figure 3) and the novel coexistence of multivortex states with giant vortex states (insets 1, 3 and 6 in figure 3) that were not predicted in [8]. Because of the large spacing between the blind holes in the experimentally considered sample, the ordered vortex lattices were not studied [5]. Figure 4 shows the vortex phase diagram for the fourth matching field (H = H 4 ), as a function of W and R. For small radii, each blind hole pins only one vortex and the remaining vortices reside between the holes forming a slightly deformed hexagonal vortex lattice [1] (together with the pinned vortices; see inset 2). With decreasing the spacing between the blind holes interstitial vortices are compressed into a giant vortex (inset 1). In the lower part of the phase diagram for the n 0 = 2 state, we found giant vortices both inside the blind hole and at the interstitial sites (inset 3). With increasing the period W first the weakly confined interstitial giant vortices split into individual vortices (inset 4) and subsequently the pinned ones as well (inset 5). In the latter case, we found compositely ordered vortex structures in the blind holes and in the interstitials, with aligned vortex dimers in the blind holes and orthogonal dimers between them. Note that the earlier molecular dynamic calculations did not show a completely ordered lattice [8]. In the n 0 = 4 state, the vortices in the blind hole form a square lattice with the same orientation as the pinning lattice (inset 8). Contrary to the cases of low matching fields (figures 2 and 3), for H = H 4 , we did not find orientational phase transitions between the pinned vortices. Vortex-antivortex states were not found at the interstitials, in contrast to what was found previously for an antidot lattice [30], as sufficient tight interstitial confinement could not be reached. Figure 5 shows typical ground state vortex configurations in a superconducting film with an array of blind holes at fractional matching fields. For smaller size of the blind holes each blind hole pins a single vortex and additional vortices are located at interstitial sites (figure 5  a,b,d)). The ordering of the interstitial vortices is different from the one obtained at the integer matching fields (see figures 2-4), which is due to the different number of vortices in the adjacent cells. Similar vortex structures have been found in samples with an antidot lattice [30]. Because of the extra interstitial vortices at every other site, the confinement does not act equally on all pinned vortices, which leads to the change in the ordering ( figure 5(c,e,f)). For example, in the case of the n 0 = 2 state at H = H 5/2 (figure 5(c)) pinned dimers are rotated by 90 • from one another and surround the interstitial vortices. Note that similar orientation of dimers was found as a stable state for the third matching field in molecular dynamic simulations [8].

Large blind holes: vortex shell structures
As we have shown in the previous section, blind holes create a confinement potential for the pinned vortices. Vortex configurations in circularly confined systems, e.g. in superconducting disks, have already been studied both theoretically [11], [31]- [34] and experimentally [35]. These studies show that vortex structures are governed by rich and subtle physics arising due to the competition between vortex-vortex interaction and confinement, thus strongly dependent on the number of vortices involved and on the radius of the superconducting disks [33]. It was shown that vortices form well-defined 'concentric' shell structures, filling these shells according to specific rules (so-called 'magic numbers'). Moreover, the formation of the shells was catalogued into a 'periodic table'. Similar shell structures were also found in different systems, such as vortices in superfluids [36], confined clusters of charged particles [37] electrons on liquid helium and dust in plasma [38,39]. Multiple occupation by single quantum vortices in large blind holes was imaged experimentally in [5] using Bitter decoration technique. It was shown that the density of vortices is considerably higher inside the blind hole than around it even if the bottom thickness is very close to the film thickness, and that vortices mostly concentrated close to the perimeter of the blind hole.
In this section, we study vortex configurations inside the blind holes at different applied magnetic fields (i.e. different occupation number n 0 ) for larger radius and period of the blind holes. To determine the stable vortex patterns, we performed field-cooling simulations each time starting from random vortex configurations for a given applied magnetic field. First, one should note that in the case of closely spaced blind holes their square arrangement interferes with the circular confinement and vortices mostly reside in the four corners of the unit cell, where the distance between the neighboring blind holes is the largest (see inset 8 of figure 4). The influence of the anisotropy of the effective vortex confinement to the structure of logarithmically interacting particles (such as superconducting vortices in extreme type-II limit) have been studied in [34,40]. It was shown that this anisotropy has a strong influence on the shell structure of vortices, e.g. shells can even be converted into straight lines. To see the effect of this complex geometry on the vortex distribution, we considered a superconducting film with an array of blind holes of radius R = 6ξ and period W = 14ξ (i.e. W − 2R = 2ξ minimal spacing between the holes). Figure 6 shows stable vortex patterns in the system for different number of pinned vortices n 0 . The first influence of the asymmetry seen in the pinning potential is the irregularities at the pentagonal structure at n 0 = 5 (b). Moreover, the (1,4) state (c) at this field becomes stable (although not a ground state), which was found earlier only in the case of vortices in superfluid helium [36]. The complex symmetry of the pinning also affects the states with larger n 0 . For example, the lowest energy states for n 0 = 6 (d), n 0 = 8 (g) and n 0 = 9 (i) have a clear square symmetry, which is the same fourfold symmetry as in the n 0 = 4 state (a). This kind of square symmetric (1,8) state was observed experimentally [35] (even more frequently than the circular (1,8) state), which is a noticeable exception from the general rule where vortices should form concentric rings. The other ground (f,j) and metastable (e,h,k) states are also clearly influenced by the square symmetry of the pinning lattice. Pure circular confinement of the blind holes is restored if we add extra vortices into the interstitial sites (see figures 6(k,l)). Figure 7 shows the vortex structure in a sample with the same radius of blind holes (R = 6ξ ) as in figure 6 but now for W = 24ξ . The effect of the interstitial vortices is clearly seen from this figure: a symmetric fourfold pattern (a) changes its orientation due to the interaction with interstitial vortices. All other states now exhibit circular symmetry. The evolution of the ground state vortex configurations is as follows: up to n 0 = 5, the vortices form regular polygons (figures 7(a,b)). Differently from the case of superconducting disks [11,35], where a maximum of five vortices fitted in the first shell, now the first eight vortices are added one by one to form the single shell (figures 7(c,d,f)). The second shell appears in the lowest energy state only at n 0 = 9 in the form of one vortex in the center and eight vortices in the second shell ( figure 7(h)), i.e. the (1,8) state. This configuration remains stable until the n 0 = 12 state is reached ( figure 7(l)), i.e. the excess vortices are added to the outer shell. Imaging experiments using the Bitter decoration technique [5] showed that the states (8), (1,7), (9) and (1,8) are stable inside the blind holes and vortices are mostly concentrated along the perimeter of the blind holes, which agrees with the present numerical calculations. We constructed a periodic table of ground state vortex configurations and compared the latter with the results obtained for individual disks. We use the standard notation for confined geometries [37] for the vortex shell state, i.e. (n, m) for the state with n vortices in the inner shell and m vortices in the outer ring. We used the following criterion to determine the number of vortices in the shells: the vortex belongs to a given shell if its distance from the center of the blind hole differs by less than 10% from the average shell radius.
In brief, the evolution of vortex shells in superconducting disks is as follows: for vorticity (i.e. number of vortices in the disk) L 5 all the vortices are arranged in a single shell and the  [33] and experiment from [35]), in the system of charged particles [37] and in superfluid He ( [36]). n 0 Blind hole Disk theory Disk exp Charged particles Vortices in He  (3,9,20) . . . (6,9,17) (1,5,11,15) 33 (3,10,20) . . . (1,6,9,17) (1,6,11,15) 34 (4,10,20) . . . (1,7,9,17) (1, 6,12,15) formation of a second shell starts from L = 6. Similarly, the formation of a third shell starts at L = 17, and of a fourth shell at L = 33. These results are contrasted in table 1 with the ground state vortex configurations inside the blind hole. The vortex structures inside the blind holes up to n 0 = 5 agree both with previous theoretical calculations in related systems and with experimental results. However, the state (6), where all six vortices are located in one shell, was found only as a metastable state in the case of superconducting disks [11,33,35], whereas the (1,5) state was found to be the ground state. Similarly, the predictions for superconducting disks and parabolically confined charged particles [39] found (1,8) as a metastable state, while in our case it becomes a ground state as in the case of vortices in He. The theoretically predicted and experimentally observed ground state configuration (2,7) for superconducting disks [11] was not found to be stable inside the blind holes either in our calculations or in the experiment [5]. The formation of the second shell starts for larger number of vortices, but also the number of vortices in the first shell grows faster than in the case of a superconducting disk. The inner shell occupancy changes from 1 to 4 with increasing n 0 from 11 to 14 (i.e. states (1,10), (2,10), (3,10) and (4,10)). The third shell in our case appears at n 0 = 19, which is larger than the theoretical predictions and the experiment. The rule for filling of the first shell (counting from the center) is similar to what was found for charged particles in a parabolic potential [37]: (i) the first shell never exceeds five vortices, second 13 and third 22; (ii) when all the shells are filled up to their maximum allowed number of vortices, a new shell consisting of only one vortex appears in the center. As we have shown in figures 6 and 7 we found that more than one vortex configuration can be stable for a given number of pinned vortices n 0 , as in the case of individual superconducting disks. In the latter system the ground state vortex configurations for a certain vorticity were recently found to depend on the dimensions of the sample: depending on the radius of the disk, the cross-overs between states with the same vorticity L are found [33]. The ground state transition between different states with the same L can also take place with changing applied magnetic field for a fixed radius of the disk [33].
Therefore, we analyzed the region of stability of the ground state configurations and performed a systematic study of all possible vortex states. As an example, we present in table 2 the number of vortices in shells inside the blind hole of radius R = 6ξ (period W = 16ξ ) corresponding to the global minima for different magnetic fields. Note that the periodic boundary conditions in equations (4), (5) restrict us to the matching fields, i.e. the magnetic field should provide flux quantization over the simulation area and the mismatch between H/H 1 and n 0 equals the number of interstitial vortices. For this particular case, the ground state transition takes place between the states (8) → (1,7) and (1,11)→(2,10) with increasing magnetic field, i.e. the number of vortices in the outer shell decreases with increasing H . The latter is a consequence of the contribution of the interstitial vortices to the effective confinement which increases with increasing H . This effect was also related to the change in the steepness of the potential energy profile in individual disks [33]. For the parameters of our system we did not find the (2,7) state as a ground state configuration, i.e. the (1,8)→(2,7) transition does not occur. However, for a larger size of the blind holes, a variety of different vortex states becomes stable (see e.g. figure 8).
For a small number of pinned vortices (n 0 4), we obtained only one stable vortex configuration. For larger n 0 , however, at least two or more stable vortex patterns are found (see figures 6(b,c), (d,e), (g,h) and (j,k), and figures 7(d,e), (f,g) and (h,i)). Note that states with the same n 0 are plotted in figures 6 and 7 in ascending order of the free energy for a given n 0 . As another example, the four vortex configurations that we found for n 0 = 19 are shown in figure 8 in ascending order of free energy. The ring structure of vortices is quite evident in these figures. The ground state configuration ( figure 8(a)) consists of one vortex in the center, surrounded by an apparent ring of five vortices, which is surrounded by another ring of 13 vortices, i.e. (1,5,13). The first excited state ( figure 8(b)) has an unusual central shell with two 'rings' of three vortices just slightly different in radius (the (3,3,13) state or the (6,13) state according to our earlier convention of number of vortices in a shell). In each of these rings, the three vortices are symmetrically placed on a perfect circle. This structure is surrounded by a ring of 13 vortices. This kind of vortex structure was obtained earlier for the system of vortices in superfluid helium [36]. The next stable vortex configuration for n 0 = 19 is the one with the two shell structure with five vortices in the inner shell and 14 vortices in the outer shell ( figure 8(c)). The three shell (1,6,12) configuration ( figure 8(d)) has the highest energy for the n 0 = 19 vortex state.
For larger superconducting disks [32] the structure of the inner shells is very close to a perfect triangular lattice, i.e. all the vortices are sixfold coordinated, whereas outer shells consist of lower-coordinated vortices. The latter leads to a competition between two types of ordering: ordering into a triangular-lattice structure in the center and ordering into a shell structure at the edge of the vortex cluster. In order to see the transition between these two vortex arrangements in our system, we calculated vortex configurations for larger number of pinned vortices, and examples of some of them are shown in figure 9. For small n 0 , we have a clear shell structure (figure 9(a)) and with increasing n 0 vortices in the center experience less presence of the boundary and start forming triangular structure, which becomes apparent already in figure 9(b). For larger n 0 , the vortices are clearly arranged in an Abrikosov lattice in the center of the bind hole, which is surrounded by a perfect circular shell at the edge of the blind hole (see figure 9(d)).
This kind of perfect triangular lattice structure for the inner part and circular ordering of the outer vortices was also found in a large system of charged particles in a parabolic confinement [37]. The density of particles in such a system is nearly constant inside and decreases slightly towards the outer shells. However, in the case of a hard-wall confinement, the density increases with radius and particles are arranged into well-pronounced shells even for large systems [37]. Particles in the hard-wall confinement potential tend to occupy first the edge positions at the wall due to their mutual repulsion. When a critical density at the edge is achieved, the remaining electrons form shell structures inside the cluster. In this sense, our system is more similar to the system of hard-wall confined charged particles, because the vortices in our case have higher density close to the perimeter of the blind holes (see figure 9). However, the density of vortices in the different shells does not increase monotonically (see figure 10) as in the case of hard-wall confined particles [37]. The inter-vortex distance in general increases slightly with radius (i.e. density decreases) except for the outer shell where vortex density is always highest.
As we mentioned above, the interstitial vortices localized between the blind holes help to restore the circular confinement of the pinned vortices (see figure 7). Moreover, at the same time the pinned vortices may impose their symmetry on the interstitial ones. In order to see this effect, we conducted simulations for larger inter-blind hole distances. Figure 11 shows the vortex structures in a superconducting film of thickness d = 0.1ξ with an array of blind holes of radius R = 10ξ and period W = 40ξ at the matching fields H = H 25 (a) and H = H 100 (b). The ring structure inside the blind hole is quite evident from figure 11(a). Due to higher density of vortices and their repulsive interaction, the vortex rings inside the blind hole transfers their shape on the interstitial vortices. As a consequence, shells of vortices are formed outside   the blind holes, as indicated by white circles in the figure. Note also that each blind hole is surrounded by a region which is practically free of vortices. The number of shells at the interstitial regions increases with increasing applied field due to the stronger vortex-vortex interaction (see figure 11(b)). However, the effect of pinned vortices on the interstitial vortex configuration fades away in the central region between the blind holes where a transition from the vortex ring structure to a triangular vortex lattice takes place.
Besides the direct visualization of vortices in superconducting samples with blind hole arrays [5], the effect of multiply trapped particles on the interstitially located ones can be studied experimentally in colloidal systems [18]. An alternative system for this study could be the clusters of macroscopic charged particles, recently developed by Saint Jean et al [21]. The simplicity of the approach, a very controllable system, and the large contrast between the particles makes this approach very promising for such investigations.

Conclusions
A rich variety of stable vortex configurations is studied in a thin superconducting film with a square array of blind holes. While the maximum number of flux quanta trapped by a blind hole is smaller than the corresponding number for antidots, we found an exciting plethora of ordered vortex states as a function of the size and period of the blind holes. These states include collective dimer, trimer, and composite states in which the vortex structures inside and outside the pinning sites exhibit an orientational ordering with respect to each other. Besides these vortex states, a combination of giant vortices in the blind holes and at the interstitials, as well as a combination of giant vortices with multivortices were obtained. The interaction of the pinned vortices inside a blind hole with the interstitial vortices gives rise to an orientational ordering with respect to the vortex structures in the adjacent pinning sites and the transitions between these ordered vortex states were obtained as a function of the blind hole parameters. Ordered vortex configurations are also found for rational (n + 1 2 ) matching fields.
For larger radius of the blind holes, the vortices are arranged into shells with an average inter-vortex distance which depends both on the radius for a particular structure and on the total number of vortices. The maximum number of vortices which fill different shells is always larger than the one in the case of superconducting disks, which we relate to a weaker confinement imposed on vortices by the blind holes. The formation of these concentric shells of vortices was studied for a broad range of blind hole occupation number n 0 . For large values of n 0 , the vortices are arranged in an Abrikosov lattice in the center of the blind hole, which is surrounded by a perfect circular outer shell of vortices. These vortex shells may even impose their symmetry on the interstitial vortices and cause a transition from Abrikosov-like vortex lattice to vortex shell structures at the interstitial sites.