# Vortex states in nanoscale superconducting squares: the influence of quantum confinement

###### Abstract

Bogoliubov-de Gennes theory is used to investigate the effect of the size of a superconducting square on the vortex states in the quantum confinement regime. When the superconducting coherence length is comparable to the Fermi wavelength, the shape resonances of the superconducting order parameter have strong influence on the vortex configuration. Several unconventional vortex states, including asymmetric ones, giant multi-vortex combinations, and states comprising giant antivortex, were found as ground states and their stability was found to be very sensitive on the value of , the size of the sample , and the magnetic flux . By increasing the temperature and/or enlarging the size of the sample, quantum confinement is suppressed and the conventional mesoscopic vortex states as predicted by the Ginzburg-Laudau (GL) theory are recovered. However, contrary to the GL results we found that the states containing symmetry-induced vortex-antivortex pairs are stable over the whole temperature range. It turns out that the inhomogeneous order parameter induced by quantum confinement favors vortex-antivortex molecules, as well as giant vortices with a rich structure in the vortex core - unattainable in the GL domain.

###### pacs:

74.78.Na, 74.20.De, 74.25.Dw## I Introduction

Vortex states in mesoscopic superconductors have been extensively studied in the past two decades, both theoretically and experimentally.meso1 (); meso2 (); meso3 (); meso4 (); meso5 (); meso6 (); meso7 (); meso8 (); meso9 (); meso20 (); meso21 (); meso10 (); meso11 (); meso31 (); gv2 (); meso32 (); meso33 (); meso34 (); meso35 (); meso36 (); meso37 (); meso38 (); meso39 (); meso40 (); meso41 (); meso42 (); meso43 () Two main interactions have been found to govern vortex behavior in a mesoscopic system. The first one is the vortex-vortex interaction, which causes vortices to form Abrikosov triangular lattices in bulk type-II superconductors. The second one is the interaction between vortices and sample boundaries, which makes vortex configurations strongly dependent on the size and geometry of mesoscopic samples-whose dimensions are of the order of the penetration depth or the coherence length . For example, in square mesoscopic samples, vortex configurations try to best match the symmetry. When there is only one vortex in the sample ( state where is the winding number or vorticity), the vortex always sits in the center of the sample in order to balance the boundary effect from all sides. For the state, two vortices sit on the diagonal such that the vortex-vortex separation is maximal in order to minimize the intervortex interaction. A giant vortex with can be induced when the boundary confinement pushes two single vortices together, as predicted theoreticallymeso2 () and observed experimentally.gv1 (); gv2 (); gv3 () For state, because of its incompatability with the four-fold symmetry, the theory predicts that the ground state corresponds to an anti-vortex sitting at the center surrounded by four vortices.anti (); anti2 () In short, the symmetry of the sample largely determines the vortex configurations when the size of the superconductor is reduced.

However, the properties of nanoscale superconductors, whose sizes are of the order of the Fermi wavelength , are very different from those of mesoscopic superconductors. This is because the distance between electronic levels becomes comparable to the superconducting energy gap due to quantum confinement Anderson (). As a consequence, the number of Cooper pairs is suppressed which leads to the quantum-size effect (QSE),qse1 (); qse2 (); qse3 () quantum-size cascades,cascades () the shell effectshell () and inhomogeneous spatial distribution of the order parameter.qse2 () The latter is the most important for the present work because it is expected to strongly influence the vortex states in nanoscale superconductors. A similar behavior was shown for an isolated vortex core, where oscillations of there order parameter on the order of the Fermi wavelength were predicted. meso44 ()

Inhomogeneous superconductivity has been studied in various systems in the last decades and shows more complex behavior than homogeneous ones. It is known that vortices tend to migrate and get pinned in areas where superconductivity is suppressed.inrmp () The reason is that it is more favorable energetically for a vortex to suppress the superconducting order parameter in a region where it has already been suppressed, although sometimes vortices can be pinned where the gap is large.pin () Some three-dimensional (3D) samples can also be treated as inhomogeneous systems.3d1 (); 3d2 (); 3Dtip (); 3d3 (); 3d4 (); 3d5 (); 3d6 () For example for a 3D tip geometry, an asymmetric vortex state can be the ground state because the thick region prevents the vortex from penetrating it.3Dtip () In multi-layered superconductors, vortices enter first and reside favorably in the weak layers. Then, vortices will penetrate into the strong layers only after weak layers become saturated and various vortex clusters and asymmetric vortex states are induced.meso6 () Also, the fabrication of anti-dots in superconductors results in a spatially varying superconducting energy gap with a barrier at the interfaces. In these systems, the combination of the giant vortex, multi-vortex and anti-vortex states can be found as ground state, which depends strongly on the detailed geometry of the antidots. ando1 (); ando2 (); ando3 ()

For conventional superconductors, ( is the Fermi wave vector and is the BCS coherence length), systems of size comparable to will not be large enough to host a vortex (being much smaller than the coherence length). However, materials with small coherence lengths, e.g. high-T cuprate superconductors, will have and therefore in such systems it is possible to obtain vortex states in the quantum confinement regime. Another such system is a graphene flake deposited on top of a superconductor. Because of the proximity effect, Cooper pairs will diffuse in graphenegr1 (); gr2 (); gr3 (); gr4 (). In graphene the scattering length is large, therefore such a system is in the clean limit. More importantly, near the Dirac point, the Fermi wavelength is very large and can be easily manipulated by doping. In other words, can be tuned, which will allow for different vortex patterns to be realized in the graphene flake in the quantum confinement regime, but for more accessible sample sizes (even above 100nm). A similar configuration was also recently proposed by Knopnin et. al in Ref.[gr5, ]. Yet another system where effects of quantum confinement on vortex matter can be probed systematically are the optically trapped cold gases xxx (), which are nowadays extremely controllable.

For studying such nanoscale systems, microscopic Bogoliubov-de Gennes (BdG) theory is required. Previous works used the BdG method to focus on isolated single vortex linesdos2 (); dos3 (); dos4 (); meso44 (), giant vorticesdos5 (); dos6 () and to describe the local density of states modifications due to vortex-vortex and vortex-boundary interactionsdos62 (); dos7 (); dos8 (); dos9 () but are in the mesoscopic limit as opposed to the nanoscale limit considered here. Although Refs.[meso43, ] and [jp, ] studied the groud state vortex states in a mesoscopic-nanoscopic crossover region by solving the BdG equations self-consistently, quantum confinement effects do not play any role. Recently, we investigatedprl () the vortex states in nano-scale superconducting squares. We found unconventional vortex states in the quantum limit due to shape-induced resonances in the inhomogeneous Cooper-pair condensate. Vortex-antivortex structures, asymmetric vortex states and vortex clusters were found as ground states over a wide range of parameters. They are distinct from previous results obtained in mesoscopic superconductors using the GL theory. However, there are still several aspects that remained unclear. For example, how does the size of the sample affect the vortex states in nanoscale superconductors? Under which conditions, does one recover the conventional GL results? Why are the antivortex states more stable in the nanoscale limit while giant vortex states are unfavorable? How do the vortex states change if temperature is increased?

In order to answer these questions, in this paper, we study vortex states in nanoscale superconducting squares systematically. Vortex states for different sample sizes, parameters and temperatures are investigated and the stability of the symmetry-induced vortex/anti-vortex molecules is discussed. More unconventional states, very different from the ones obtained within GL theory, are found. This study is therefore fully complementary to what is known for vortex matter in superconductors.

The paper is organized as follows. In Sec. II, we introduce the theoretical approach, i.e. the BdG approach for a square geometry. In Sec. III, we present the inhomogeneous superconducting state in the absence of the magnetic field in order to better understand the QSE in nano-scale superconductors. In Sec. IV, the ground states and metastable states are studied at zero temperature and on the sample size dependence of the vortex states is discussed. In Sec. V, the finite temperature ground states are studied. In Sec. VI, we discuss the generation of vortex/anti-vortex molecules and study the structure of the vortex core. Finally, we summarize our findings in Sec. VII.

## Ii Theoretical Formalism

In the presence of a magnetic field, the Bogoliubov-de Gennes (BdG) equations read

(1) |

where () are electron-(hole-)like quasi-particle eigen-wave functions and are the quasi-particle eigen-energies. The single-electron Hamiltonian reads with being the Fermi energy and the vector potential (we consider a gauge such that ). Furthermore, we take the electron band-mass to be isotropic (i.e., ) and a circular Fermi surface. The pair potential is determined self-consistently from the eigen-wave functions and eigen-energies:

(2) |

where is the coupling constant, is the cutoff energy, and is the Fermi distribution function at temperature .

We consider now a two-dimensional superconducting square with the size . The confinement imposes Dirichlet boundary conditions (i.e. ) such that the order parameter vanishes at the surface. In an extreme type-II superconductor (and/or very thin sample), it is reasonable to neglect the contribution of the supercurrent to the total magnetic field.

In this case, the free energy FreeEg1 (); FreeEg2 () of the system is given by

(3) | |||||

where the spatial dependence of and is implicit.

In order to solve Eqs. (1) and (2) numerically, we expanded the wave functions and as

(4) |

where the basis set

(5) |

is the basis eigen-states of the Hamiltonian in the absence of the magnetic field. The corresponding eigen-energies of such states are . Through comparison with results obtained by using the finite difference method, we found that the results are converged when we include the states with energies as large as times the cutoff energy , i.e. where is the flux quantum and the flux through our sample.

Sample | I | II | III |

4 | 9 | 25 | |

1.245 | 1.85 | 3.14 | |

2.04 | 3.09 | 5.07 |

In our numerical investigations, we restrict ourselves to the three materials (or samples) with parameters given in Table 1. For convenience, we measure the distance in units of the bulk coherence length at zero temperature and the energy in units of . Here, where is the Fermi velocity and is the bulk value of the order parameter at zero temperature. Note that and are not independent when up to a constant.

To find the different vortex configurations, which include all stable states, we search for the self-consistent solutions in the following two steps. (1) Global scanning: Starting from any reasonable vortex state (usually, we start from the Meissner state at ), we slowly sweep up/down the applied flux with regular intervals and recalculate the superconducting states each time, until a new state is found. Then, we repeat the sweeping process from the new state until no new vortex configurations appears. (2) Special initial states: Starting from usual states obtained in GL theorymeso10 (), we sweep up/down the applied flux. If a new state appears, we repeat the step (1). In such a way, we are able to trace back and forth all found vortex states in the whole region of their stability and make sure that the usual GL states are always considered.

The minimum tolerance in the change of the order parameter between two steps in the self-consistent iteration is

(6) |

where and is the order parameter at the and iteration. We use the absolute tolerance since the relative tolerance can be abnormally high in the vortex core where . This is quite strict when we compared to but is necessary in order to ensure the precision in finding the true ground states in the BdG calculation. We show an example the evaluation of the vortex configuration for sample II and at flux . Based on the GL theory, the ground state of for such a system should be giant vortex state or multi vortex state (for larger squares) where two vortices sit on the diagonal of the sample. We start the calculation with such the multi vortex state as initial state. As seen from the Fig. 1, the vortices merge into a vortex anti-vortex molecule and the result converges quickly. The error between each steps reach as low as . However, the error increases gradually when we continue the self-consistent procedure. After the second-order phase transition, the new state with two vortices sitting parallel to one of the sides has lower energy. Finally, the symmetry of the state does not change and the error is always around which comes from non-physical factors, i.e. numerical accuracy.

In the calculation, we found some situations where results do not converge and this usually comes from the change in the number of the quasi-states contained in the Debye window. To avoid this, we set the smearing energy for the quasi-states. Then, the self-consistent condition reads

(7) |

where is the Fermi distribution function. The choice of is empirical. It should be enough small in order not to affect the results. Meanwhile, it should be enough large to make the result converge through the iteration. Our experience show that is suitable for our current work.

## Iii Spontaneous Inhomogeneous Superconductivity induced by Quantum confinement

First, we study the spatial distribution of the order parameter as a function of sample size in the absence of the magnetic field since the vortex configurations will be strongly affected by that distribution once the magnetic field is applied. In Fig. 2 contour plots of the order parameter and the corresponding diagonal profiles are shown for sample I-III with different sizes. As expected, all the order parameters are fourfold symmetric and show Friedel-like oscillations in space which result in four well defined peaks at each corner. Hence, the superconductivity is inhomogeneous. For example, Figs. 2 (a)-(c) show the results for sample I with size , respectively. For , there are three oscillations in the order parameter along the diagonal and the resonant central peak prevents vortex from sitting here. However, the profile of the order parameter can be changed dramatically when the size of the sample changes. For [shown in Figs. 2(b)], we found that the central peak disappears and a relatively flatter area generates in the center. When is increased to , the flat area enlarges and the Friedel-like oscillations can be neglected at center when we compared it with oscillations near boundary. It indicates that the oscillations of the order parameter result from the quantum confinement effect (or boundary effect).

In order to study the quantum size effect on the order parameter, we show the amplitude of the order parameter in the center of the sample and the spatially averaged value over the whole sample for the sample I with sizes in Fig. 3. As seen, changes dramatically with increasing and converge to when . At the same time, the increases gradually with increasing. In principle, both of the parameters will converge as where the boundary effect can be totally neglected. Since and show strong quantum size effect between and , we limit ourselves in following sections to study the samples for these particular size .

The profile of is also strongly affected by . Figs. 2(d) and (e) show results for sample II with size and and (f) is for sample III with size . Comparing to the sample I with the same size, we found that the wave number and the amplitude of the oscillations along diagonal are larger. It indicates that the superconducting order parameter shows more inhomogeneous behavior with larger .

## Iv Vortex configurations at zero temperature

In this section, we consider the zero temperature case for which the system is always in the quantum limit (). First, we study the sample I with size and show in Fig. 4 the free energy of the stable vortex states for flux . Different curves (colors) represent states with different winding number and the states among them which reached the ground state are marked by dots. Vertical lines and shadows show the flux range for each state as the ground state. The top dashed line stands for the free energy of the system in the normal state when the coupling constant is set to zero. When compared with the GL theorymeso10 (), one of the differences is that the free energy of the normal state depends on the magnetic field while it is a constant in GL theory. The reason is that the energy levels of the confined electrons are different for different magnetic fields. In our case, the change of the energy is relatively small when compared with the energy gap (energy difference between the normal state and the superconducting state) especially in weak fields. Although the shown energy curves look conventional, there are significant differences with the GL theory.meso10 ()

By sweeping the magnetic field up and down, we can get the full energy spectrum and the corresponding vortex states. For a certain magnetic field, it is common to have more than one converged solution. The lowest energy state is the ground state while the states with higher energy are referred to as metastable states. In Fig. 5, we show the contour plot of the absolute value of the order parameter of the corresponding ground states for various winding numbers.

As can be seen from Fig. 4, the system favors states with winding numbers and because they have relatively large ground state flux range (excluding the Meissner state). From Fig. 5, we observe that these states have fourfold symmetry which is compatible with the sample geometry. One interesting feature of this system is the richness of metastable states. These states appear for all winding numbers except and . The number of metastable states reaches a peak for and equals . From the free energy curves, we notice that the energy difference between the ground state and the metastable states can sometimes be very small. This makes the ground state difficult to find in simulations unless we sweep the field up and down many times. We also note that most metastable states are focused at lower magnetic fields from the corresponding ground state flux range. The reason is that vortices get easily stuck at the boundary due to the pronounced oscillations of the order parameter. For the same reason, their stability range is narrower due to asymmetry.

The number of metastable states decreases for higher . In this case, the shorter distances between vortices cause strong interaction between them. This limits the choice of stable positions for vortices and therefore the number of metastable states are lower. Due to this reason, metastable states are less favorable for smaller samples because of easy saturation with vortices. For example, no metastable states were found in sample I with size .

### iv.1 Ground states

Next we discuss the ground states configurations for different vorticities in more detail. For the state, shown in Fig. 5, the vortex sits at the center of the sample which is compatible with the conventional picture. Although such result is observed for different parameters of the sample, the state with diagonal location of the vortex can also be found in some cases.prl () It is clear that the order parameter around the vortex core is suppressed and the profile shows the competition between symmetry and . It means that the vortex has long range (longer than ) interaction with other objects such as the other vortices and/or boundaries.

By increasing the flux to , the ground state shifts from to and two vortices sit along the diagonal. This again coincides with the result from GL theory, and results from the competition between the confinement imposed on vortices by the Meissner currents and the vortex-vortex repulsion. In this case, these effects can be clearly seen from the profile of the order parameter, especially from the suppressed area around the vortices. The vortex-vortex interaction suppresses the order parameter mostly in the area between them. This can not be found in GL theory because the order parameter is always smooth and changes slowly in space. The state with shown in Fig. 5 becomes the ground state when the fieldapplied flux is between and . It resembles the multi-vortex state obtained within GL theory where the three vortices are at the apices of a equilateral triangle. However, the perpendicular bisector of the triangle always coincides with one of the diagonals of the square sample in GL theory while it is parallel to one of the edges in our case. This is because the grid-like pattern in the inhomogeneity of the order parameter imposes preferential positions for the vortices inside the square. The state has a similar feature but the configuration is compatible with the GL result.

From these states, we conclude that when the GL vortex configuration, which minimizes the vortex-vortex and vortex-boundary interaction, matches the oscillation pattern due to quantum confinement, then the state has a wider flux stability range.

Two ground states, and , are found for in the flux ranges and , respectively. Both of them have a pentagonal vortex configuration. This is because the particular shape resonance at the considered field causes the order parameter to be peaked at the center. Therefore, it costs energy for vortices to sit in the center of the sample. For the same reason, the ground states and do not have vortices in the center. Moreover, when , vortices start to be compressed in the sample. If they do not form giant vortices, they will be very close to each other and form string-like structures [see state].

States with , , and shown in the panels of Fig. 5 have a common feature, as all of them keep the fourfold symmetry. contains a antivortex at the center while and have a single vortex and a giant vortex with at the center, respectively. The state with is the only ground state which contains an antivortex. The antivortex is closely surrounded by four vortices and forms the core structure for . The outer shell is formed by the remaining four individual vortices sitting at four corners. States with and contain vortex dimers, i.e. two vortices close to each other at each corner. The fourfold symmetry makes both former states have a larger ground state flux range. also keeps the symmetry but the energy is sometimes even higher than the state , which has only symmetry. The reason is that the giant vortex costs extra energy.

In order to visualize the changes in the ground states when key parameters change, we plot the phase diagram for samples I-III for and in Fig. 6. Different shadings of blocks in Fig. 6 indicate different vortex types. The plain white background represents conventional multi-vortex states as found within the GL theory, while the blue background with square grid represents giant vortices, also compatible with the result obtained from GL theory. Asymmetric vortex states attained only by BdG theory are represented by yellow background with horizontal grid pattern. States containing parallel vortex chains, represented by orange background with vertical grid pattern, and part of the vortex-antivortex molecules represented by pink(grey) background, are new compared to GL theory.

We can conclude from Fig. 6 that the quantum size effect is important not only for the transition field and the stability range of the different vortex states, but also for the vortex configurations. The reason is that the oscillation patterns of the order parameter are very sensitive to and may cause totally different behaviors even for two samples of identical size.

As can be seen from Fig. 6(a) and (b), the vortex phase transition fields vary greatly with sample size except for the phase boundary between and state. All the phase boundaries oscillate with . We find that some samples favor vortex states with even winding number while other disfavor them. For example, sample I with favors state whereas the one with disfavors to the point of non-existance. When becomes large, the transition fields start to converge and the flux stability range of each state will be roughly one flux quantum.

When compared to samples I and II, sample III (with large ) shows a more conventional picture. Moreover, the phase transition field increases only slightly with increasing . It means that for large , the quantum size effect on the transition field, at least when winding number is small, can be neglected. Nevertheless, a plethora of different vortex configurations is found. For example, an asymmetric state is found in all three samples. Moreover, sample I always shows asymmetric states when . One other interesting phenomenon is that, for , sample II can host all five types of vortex states with increasing.

From the phase diagram in Fig. 6, we notice that nanoscale superconductors favor anti-vortices and disfavor giant vortices. For example, the giant vortex state appears for in sample II with . Based on GL theory, only smaller samples will exhibit a giant vortex configuration. However, when , the two individual vortices form an asymmetric vortex state. On the other hand, we find that the probability of forming anti-vortices is much higher than in GL theory. In GL case, antivortex states usually appear for when four vortices are at the four corners and surround a centered antivortex. Usually the distance between vortices and the antivortex is small (less than ). In the BdG calculation, at least two more antivortex states can be found. One is still the state and the antivortex is still at the center but four vortices are at the edges instead of the corners. This configuration was first presented by us in Ref. [prl, ] where we showed that the size of the vortex/anti-vortex(V-aV) molecule is larger than the one obtained with GL theory. The other antivortex state appears for . In this case, the four vortices still sit at four corners but the centered antivortex carries two flux quanta, e.g. it is a giant antivortex! Due to the strong vortex-antivortex interaction, the size of such V-aV molecule is small.

### iv.2 Metastable states

Next we will briefly discuss the metastable states of this system. To do this, we start from sample I with . Metastable states are important in the BdG formalism because the energy difference between the ground and the metastable states can be very small. This suggests that these states could be easily found in experiments. Alternatively, some metastable states can become ground states as the parameters are changed.

All six found metastable states for and their free energy curves are shown in Fig. 7. The state is similar to the ground state , but rotated over , hence their free energies are very close to each other. Actually, the difference in the orientation of the vortex pattern always results in a small difference in energy. State (f) is not obtained within the GL theory. In our case, due to the shape-resonant inhomogeneity of the order parameter, the rotation of the vortex pattern to the ground-state configuration is prevented by the spatial oscillations of the order parameter.

The metastable state is only stable at higher field and its free energy is very close to the ground state . Therefore we zoomed on the energy difference in the inset of Fig. 7. From the figure, one sees that the energy of the state is lower than the ground state, , when the applied flux is larger than . In fact, the state can exist even up to . From the vortex configuration shown in Fig. 7, we find that this state is a giant vortex. Such a state has been predicted by the GL theory because the magnetic field pushes the two vortices towards each other and makes them merge into a giant vortex. Usually, the phase transition between the multi vortex state and the giant vortex state is continuous (second-order). However, the barrier induced by the inhomogeneity of the order parameter leads to a first order phase transition in our case. One more difference between the BdG giant vortex and the one in GL theory is its core structure. Due to the shape resonances, the contour plot of the core shows a diagonal cross shape while the giant vortex core in the GL case is always circular. Furthermore, the giant vortex state in our results has two allotropes: see state , compared to the state . The state exists up to higher field while only exists in lower field. Hence, the size of the giant vortex seen in is larger than the one seen in . Another difference between them is the orientation of the core. has diagonal cross shape while the state has edge cross shape.

Other three metastable states , and are observed only in lower field. They have in common the fact that at least one vortex is stuck at the boundary since the Meissner current pushes the vortex outward at low fields. It is obvious that they have lower energy when the applied field is lower. The energy of state is always lower than the one of state because of the longer distance between the vortices. As is usual in mesoscopic superconductors, vortices in these states avoid to be located at the very corners of the sample, due to strong local superconductivity there.

For the metastable states, the results are summarized in Fig. 8. The state is a V-aV state and exists only in higher field while the giant vortex only exists in lower field. Note again that the energy of the giant vortex states is much higher than the other metastable states with three single vortices. has lowest energy for around , when there is one vortex located at the boundary. States and follow when the field decreases and there are two and three vortices stuck at the boundaries, respectively. States and are disfavored and have higher energy due to the close distance between vortices. Note again that no vortex sits at the corners in this states.

States with show a wide ground state flux range and different metastable states, which is the largest variety of all vortex states. From Fig. 9, we find that the metastable states concentrate around the applied flux .

For low fields, we conclude again that vortices are close to the surface and these states are always the lowest energy state for a given state at low fields. States and have symmetry and all of the four vortices are trapped close to the boundary. Please note that in the state four vortices sit at the corners. This kind of state is rare in the BdG results because corners give the highest potential energy contribution for vortices. From the free energy curve of the state , we can find that the slope of the energy curve is opposite to the other states [such as ] in this field range. This indicates that vortices are repelled by the Meissner current in order to balance the inward force the vortices experience from the corner. When the field is too low, a vortex is expelled from the sample and the state jumps to a state. The vortex configuration has been found experimentally in conventional mesoscopic superconductorsmeso9 () but was a result of the presence of pinning sites. This state can not be obtained in plane squares within the GL theory . Another state to notice is whose energy curve does not cross any other curve. When the superconductor is in this state and the field is swept down, this state will be the first to jump to the state. This is understandable from the vortex configuration of because the vortex at the corner is easily expelled when field is lowered.

At high fields only one metastable state exists . It can be seen as the state obtained after a degrees rotation of the ground state-. This is a consequence of the fact that the inhomogeneous pattern of the order parameter changes with field. At such a high field, the corner vortex position in state becomes unstable, which forces the vortices to sit at the edges, similar to the case. At the same time, the strong field pushes vortices closer together so that the distance between vortices in is shorter than the one in the ground state .

## V Vortex states at finite temperature

So far, all our calculations were done at zero temperature, , where is the bulk critical temperature at zero flux . In what follows, we investigate the effect of temperature on the vortex configuration. First, we show all the vortex states for the flux range , for sample I with size at where the system is NOT in the quantum limit since . The corresponding free energy curve as a function of flux is presented in Fig. 10. Contrary to the results for , which were shown in Fig. 4, the figure looks more conventional (similar to the results obtained by GL theory in Ref. meso10, ) and there is only one stable state for each winding number . Moreover, only giant vortex states are found in this case for . For the size of the square sample considered here, , GL theorymeso10 () predicts that multi-vortex states should exist. Here we find instead that multi-vortex states are absent since increases as temperature increases.

In order to see how temperature affects the coherence and the profile of the superconducting order parameter, we show the order parameter for and in Fig. 11. The diagonal profile of at shows the strongest Friedel-like oscillations. As temperature increases to , the profile is similar to the one obtained at zero temperature, but with less oscillations at the vortex core. Both cases are in the quantum limit and the order parameter shows rapid variation in the core. When temperature reaches , we find that both the average and the oscillations of the absolute value of the order parameter are suppressed which indicates that the vortex states become more conventional. Finally, the order parameter is smooth at and the GL results are approached.

As seen from the Fig. 11, the coherence length, which represents the vortex core radius, increases with increasing temperature. As defined by Kramer and Peschkramer (), we calculate the coherence length as

(8) |

where is the distance to the vortex core. We plot in Fig. 12 as a function of temperature, . As discussed in Ref. [dos3, ], can be described by when is close to ( in our case). In the intermediate temperature regime, there is a substantial suppression of the coherence length because of the bound states. At low temperature, the shrinkage of the coherence length stops and saturates when the system is in the quantum limit. Note that at is around three times larger than the one at zero temperature. This explains why only giant vortex states can be found at such temperatures.

Fig. 13 shows the order parameter for sample I at for in panel (a) and in panel (b), respectively. Both are giant vortex states and the symmetry grid pattern is strongly suppressed. As can be seen from the figure, the vortex cores show perfect circular symmetry, which is in agreement with the results from GL theory. Of course, the size of the vortex core shown in panel (b) is larger than the one shown in panel (a) because its vorticity is larger.

Finally, we end this section with the (T-) phase diagram for lower fields for sample I with . This is shown in Fig. 14. The thick black curve indicates the phase boundary between the superconducting and the normal state. When the system is in the quantum limit, for these parameters, only unconventional vortex states, such as asymmetric and states and edge-parallel states, are found as ground states. When temperature increases, the vortex states become conventional and the symmetry of the states is always preserved. Note that the asymmetric state goes through a continuous phase transition to the symmetric state, which means the vortex moves gradually as the temperature changes. However, for higher winding numbers, the system usually goes through a first order phase transition. For example, the phase transition between the parallel vortex state and the giant vortex state of is of first order. This is different from the GL result, where vortices merge into a giant vortex through a continuous phase transitionbenxu (). For the state, we note that the ground state flux range for the four-fold symmetric V-AV state is larger than the asymmetric one due to the compatibility of its symmetry with the geometry of the sample.

Concluding this section, higher temperature: 1) makes vortex states look more conventional (closer to the GL results); 2) smoothens the order parameter; 3) suppresses the influence of the oscillation of the order parameter and 4) increases the superconducting coherence length . As a consequence, the number of metastable states is also lowered. The effect of temperature is very different (more complex) from the effect obtained by simply changing the effective size of the sample as is usually done within the GL theory.

## Vi Giant anti-vortex and the structure of the vortex core

In this section, we discuss the appearance and stability of anti-vortex states in the BdG theory in order to explain the existence of the giant antivortex. Actually, such a state was already found in Ref. anti2 () through the linear GL method by introducing artificial pinning.

From the phase diagram shown in Fig. 6, we found that anti-vortex states are surprisingly stable within BdG theory. This is due to the fact that the grid pattern oscillation of the order parameter gives an additional contribution to the symmetry of the vortex states and therefore, in a square sample, the symmetry is enhanced. The other reason to form an anti-vortex is that the oscillations induced by the order parameter are seen in the vortex core where the order parameter is already suppressed. These oscillations can easily lead to a shift in the phase of the order parameter by and, thus, result in the formation of vortex-antivortex molecules.

In order to explain this, we first discuss briefly the vortex profile in the GL theory. Fig. 15 shows a schematic diagram of the vortex states for different winding number in GL theory. The diagonal profiles of the order parameters vary smoothly in space and the vortex emerges where the order parameter vanishes. Note that the phase of the order parameter is adjusted such that along the diagonal the order parameter is real. Panel (a) from Fig. 15 shows the simplest case when only one vortex sits at the center. As can been seen from the figure, the order parameter changes sign, which indicates the phase shift of the order parameter. The profile is an odd function and near the vortex core. Panels (b) and (c) from Fig. 15 show the diagonal profiles for . Both profiles are even functions due to the phase shift between the opposite corners. The order parameter exhibits property. When there is only one root, as can be seen from panel (b), the vortex is a giant one. When there are two roots, as shown in panel (c), the configurations are multi-vortex states. Similarly, the profiles of the order parameter shown in panels (d) and (e) from Fig. 15 for show a spatial dependence. One root means that we have a giant vortex state whereas three roots represent a vortex-antivortex configuration. Note that, in order to generate the central anti-vortex, the order parameter has to oscillate around the center of the square.

Now let us move to nano-size superconductors where the BdG theory has to be used and the spatial oscillation of the order parameter cannot be neglected. As can be seen from Fig. 16 the oscillation plays an important role in generating vortices, especially when the value of the order parameter is comparable to the amplitude of the oscillation. For instance, panels (a) and (b) from Fig. 16 show the symmetric and the asymmetric vortex states. The reason for the appearance of the asymmetric vortex state is the fact that the order parameter has an odd number of oscillations across the diagonal. Thus, the vortex can not sit at the center. For the states, panel (c) shows a giant vortex state where the sign of the profile of the order parameter is always positive across the diagonal. However, due to the oscillations, the case shown in panel (d) of Fig. 16 can easily exist and shows a giant anti-vortex () at the center. Further, the configurations show a large diversity for a fixed winding number . When , the configuration can be , , , and so on. Panels (e) and (f) from Fig. 16 are for states. Apparently, they are similar to the GL case shown in Fig. 15, but the probability of the occurrence of the V-aV state is much larger than in GL case. The reason is that the result with one root is just a special case while the general case shows oscillations at and around the vortex core.

The V-aV molecules do not only exist for smaller winding number , but they can also appear for large in the BdG results. Fig. 17 shows an example for sample III with =7 at and with a winding number . We find that vortices concentrate in the central dark (blue) area where the order parameter is strongly suppressed. From the phase of the order parameter, which is shown in Fig. 17(b), the total winding number is found but it is difficult to distinguish each vortex. After a careful analysis, we plot schematic diagram of the vortex configuration in Fig. 17(c). The dark(blue) dots and open diamonds indicate vortices and anti-vortices, respectively. As can be seen, the lattices of vortices and antivortices are nested within each other. Since anti-vortices attract vortices, all the vortices ( vortices and antivortices) can be condensed in the cental area of the sample. This picture becomes more accurate when is large. The stronger the oscillations of the order parameter the more V-aV pairs are generated. However, the size of the V-aV pair can only be of the order of the Fermi wavelength. Thus, it will be very hard to detect them in experiments. This is why these states are mostly treated as a giant vortex in conventional superconductors. In other words, the suppressed central area of the order parameter, after coarse graining, will look like a giant vortex with .

## Vii Conclusion

To summarize, we investigated the vortex states in a nanoscale superconducting square for different sizes , parameters , and temperatures . First, we found that the inhomogeneous pattern of the order parameter in the absence of magnetic field strongly depends on and the size . This oscillation pattern will give an additional contribution to competing effects that determine the vortex configurations when the field is applied. Due to the inhomogeneous order parameter induced by the quantum topological confinement, samples with different and will favor different winding numbers .

We find unconventional vortex states such as asymmetric, edge-parallel and vortex-antivortex states as the ground state of our nanoscale system. These were never seen in the Ginzburg-Landau approach. The inhomogeneous pattern of the order parameter, especially the strong oscillation at the boundaries causes additional potential wells for vortices which in turn generates a lot of metastable vortex states. Furthermore, in the quantum limit, nano-size superconductors favor vortex-antivortex molecules while disfavoring giant vortex states.

We observe that vortex ground states and the phase transition fields are very sensitive to changes in the parameter , size and temperature . This is a direct consequence of the quantum size effect. However, this effect is suppressed when the size is large or when temperature is high. In this case most metastable states become unstable and the ground states become compatible with GL theory.

For high magnetic fields, vortex-antivortex pairs can be easily found when is large because the absolute value of the order parameter becomes smaller than the amplitude of its oscillations. However, detection of such states is beyond the current experimental abilities.

The peculiar vortex states uncovered in the present work should be observable in superconducting systems where is small. Such systems could be high- superconducting nano-grains for which the coherence length is small or cold-atom condensates with small , i.e. large Fermi wavelength. Of special interest could be hybrid systems made of superconducting substrates and graphene sheets for which the Fermi wavelength is highly tunable near the Dirac point. Future work could also address the fundamental vortex-vortex and vortex-antivortex interactions for systems with a small , for which the oscillations of the order parameter on the order of become important.

## Viii Acknowledgments

This work was supported by the Flemish Science Foundation (FWO-Vlaanderen) and Methusalem Funding of the Flemish government.

## References

- (1) A. K. Geim, I. V. Grigorieva, S. V. Dubonos, J. G. S. Lok, J. C. Maan, A. E. Filippov, and F. M. Peeters, Nature (London) 390, 259 (1997).
- (2) V. A. Schweigert, F. M. Peeters, and P. S. Deo, Phys. Rev. Lett. 81, 2783 (1998).
- (3) J. J. Palacios, Phys. Rev. Lett. 84, 1796 (2000).
- (4) B. J. Baelus, F. M. Peeters, and V. A. Schweigert, Phys. Rev. B 63, 144517 (2001).
- (5) M. V. Milošević, S. V. Yampolskii, and F. M. Peeters, Phys. Rev. B 66, 024515 (2002).
- (6) O. Bourgeois, S.E. Skipetrov, F. Ong, and J. Chaussy, Phys. Rev. Lett. 94, 057007 (2005).
- (7) L. F. Chibotaru, A. Ceulemans, V. Bruyndoncx, and V. V. Moshchalkov, Phys. Rev. Lett. 86, 1323 (2001).
- (8) V.R. Misko, V.M. Fomin, J.T. Devreese, and V.V. Moshchalkov, Phys. Rev. Lett. 90, 147003 (2003).
- (9) H. J. Zhao, V. R. Misko, F. M. Peeters, V. Oboznov, S. V. Dubonos, and I. V. Grigorieva, Phys. Rev. B 78, 104517 (2008).
- (10) J. Bonča and V. V. Kabanov, Phys. Rev. B 65, 012509 (2001).
- (11) T. Mertelj and V. V. Kabanov, Phys. Rev. B 67, 134527 (2003).
- (12) B. J. Baelus, and F. M. Peeters, Phys. Rev. B 65, 104515 (2002).
- (13) A. S. Mel’nikov, I. M. Nefedov, D. A. Ryzhov, I. A. Shereshevskii, V. M. Vinokur, and P. P. Vysheslavtsev, Phys. Rev. B 65, 140503 (2002).
- (14) K. Tanaka, I. Robel, and B. Jankó, PNAS 99, 5233 (2002).
- (15) D.S. Golubović, M.V. Milošević, F.M. Peeters, and V.V. Moshchalkov, Phys. Rev. B 71, 180502 (2005).
- (16) M.R. Connolly, M.V. Milošević, S.J. Bending, J.R. Clem, and T. Tamegai, EPL 85, 17008 (2009).
- (17) F.R. Ong, O. Bourgeois, S.E. Skipetrov, and J. Chaussy, Phys. Rev. B 74, 140503 (2006).
- (18) P.J. Pereira, V.V. Moshchalkov, and L.F. Chibotaru, Phys. Rev. E 86, 056709 (2012).
- (19) R. Geurts, M.V. Milošević, and F.M. Peeters, Phys. Rev. B 81, 214514 (2010).
- (20) A.S. Mel’nikov and V.M. Vinokur, Nature (London) 415, 60 (2002).
- (21) S.-Z. Lin, T. Nishio, L.N. Bulaevskii, M.J. Graf, and Y. Hasegawa, Phys. Rev. B 85, 134534 (2012).
- (22) T. Nishio, T. An, A. Nomura, K. Miyachi, T. Eguchi, H. Sakata, S. Lin, N. Hayashi, N. Nakai, M. Machida, and Y. Hasegawa, Phys. Rev. Lett. 101, 167001 (2008).
- (23) T. Cren, D. Fokin, F. Debontridder, V. Dubost, and D. Roditchev, Phys. Rev. Lett. 102, 127005 (2009).
- (24) J. Tóbik, V. Cambel, and G. Karapetrov, Phys. Rev. B 86, 134433 (2012).
- (25) P.J. Pereira, L.F. Chibotaru, and V.V. Moshchalkov, Phys. Rev. B 84, 144504 (2011).
- (26) H.T. Huy, M. Kato, and T. Ishida, Supercond. Sci. Technol. 26, 065001 (2013).
- (27) X.-H. Hu, A.-C. Ji, X.-G. Qiu, W.-M. Liu, Eu. Phys. J. B, 79, 473 (2011).
- (28) T. Cren, L. Serrier-Garcia, F. Debontridder, and D. Roditchev, Phys. Rev. Lett. 107, 097202 (2011).
- (29) A. Kanda, B. J. Baelus, F. M. Peeters, K. Kadowaki, and Y. Ootuka, Phys. Rev. Lett. 93, 257002 (2004).
- (30) L. F. Chibotaru, A. Ceulemans, V. Bruyndoncx, and V. V. Moshchalkov, Nature (London) 408, 833 (2000).
- (31) R. Geurts, M. V. Milošević, and F. M. Peeters, Phys. Rev. Lett. 97, 137002 (2006). ibid. Phys. Rev. B 75, 184511 (2007).
- (32) P. W. Anderson, J. Phys. Chem. Solids 11, 26 (1959).
- (33) A. A. Shanenko, M. D. Croitoru, and F. M. Peeters, Phys. Rev. B 75, 014519 (2007).
- (34) M. D. Croitoru, A. A. Shanenko, and F. M. Peeters, Phys. Rev. B 76, 024511 (2007).
- (35) M.D. Croitoru, A. Vagov, A.A. Shanenko, and V.M. Axt, Supercond. Sci. Technol. 25, 124001 (2012).
- (36) A. A. Shanenko, M. D. Croitoru, and F. M. Peeters, Phys. Rev. B 78, 024505 (2008).
- (37) S. Bose, A. M. Garcia-Garcia, M. M. Ugeda, J. D. Urbina, C. H. Michaelis, I. Brihuega, and K. Kern, Nat. Mater. 9, 550 (2010).
- (38) M. Machida and T. Koyama, Phys. Rev. Lett. 90, 077003 (2003).
- (39) G. Blatter, M. V. Feigelman, V. B. Geshkenbein, A. I. Larkin, and V. M. Vinokur, Rev. Mod. Phys. 66, 1125 (1994).
- (40) D. Valdez-Balderas and D. Stroud, Phys. Rev. B 76, 144506 (2007).
- (41) M. M. Doria, A. R. de C. Romaguera, and F. M. Peeters, Phys. Rev. B 75, 064505 (2007).
- (42) B. Xu, M. V. Milošević, and F. M. Peeters, Phys. Rev. B 77, 144509 (2008).
- (43) Y. Chen, M.M. Doria, and F.M. Peeters, Phys. Rev. B 77, 054511(2008).
- (44) B. Xu, M. V. Milošević, and F. M. Peeters, New J. Phys. 11, 013020 (2009).
- (45) B. Xu, M. V. Milošević, and F. M. Peeters, Phys. Rev. B 82, 214501 (2010).
- (46) C.-Y. Liu, G. R. Berdiyorov, and M.V. Milošević, Phys. Rev. B 83, 104524 (2011).
- (47) M. A. Engbarth, S. J. Bending, and M. V. Milošević, Phys. Rev. B 83, 224504 (2011).
- (48) G. R. Berdiyorov, B. J. Baelus, M. V. Milošević, and F. M. Peeters, Phys. Rev. B 68, 174521 (2003).
- (49) G. R. Berdiyorov, M. V. Milošević, and F. M. Peeters, Phys. Rev. B 74, 174512 (2006).
- (50) G. R. Berdiyorov, M. V. Milošević, and F. M. Peeters, New J. Phys. 11, 013025 (2009).
- (51) H. B. Heersche, P. Jarillo-Herrero, J. B. Oostinga, L. M. K. Vandersypen, and A. F. Morpurgo, Nature (London) 446, 56 (2007).
- (52) H. Tomori, A. Kanda, H. Goto, S. Takana, Y. Ootuka, and K. Tsukagoshi, Physica C: Superconductivity 470, 1492 (2010).
- (53) A. M. Black-Schaffer and S. Doniach, Phys. Rev. B 78, 024504 (2008).
- (54) L. Covaci and F.M. Peeters, Phys. Rev. B 84, 241401 (2011).
- (55) N. B. Kopnin, I. M. Khaymovich, and A. S. Mel’nikov Phys. Rev. Lett. 110, 027003 (2013).
- (56) S. Aubin, S. Myrskog, M. H. T. Extavour, L. J. LeBlanc, D. McKay, A. Stummer, and J. H. Thywissen, Nature Physics 2, 384 (2006).
- (57) L.-F. Zhang, L. Covaci, M.V. Milošević, G.R. Berdiyorov, and F.M. Peeters, Phys. Rev. Lett. 109, 107001 (2012).
- (58) I. Kosztin, S. Kos, M. Stone, and A. J. Leggett, Phys. Rev. B 58, 9365 (1998).
- (59) A. Spuntarelli, P. Pieri, and G. C. Strinati, Physics Reports 488, 111 (2010).
- (60) C. Caroli, P. G. de Gennes, and J. Matricon, Phys. Lett. 9, 307 (1964).
- (61) F. Gygi and M. Schlüter, Phys. Rev. B 43, 7609 (1991).
- (62) N. Hayashi, T. Isoshima, M. Ichioka, and K. Machida, Phys. Rev. Lett. 80, 2921 (1998).
- (63) S. M. M. Virtanen and M. M. Salomaa, Phys. Rev. B 60, 14581 (1999).
- (64) C. Berthod, Phys. Rev. B 71, 134513 (2005).
- (65) A.S. Melnikov and M.A. Silaev, JETP Lett.83, 578 (2006).
- (66) A. S. Melnikov and V. M. Vinokur, Phys. Rev. B 65, 224514 (2002).
- (67) A. S. Melnikov, D. A. Ryzhov, and M. A. Silaev, Phys. Rev. B 78, 064513 (2008).
- (68) A. S. Melnikov, D.A. Ryzhov, and M. A. Silaev, Phys. Rev. B 79, 134521 (2009).
- (69) H. Suematsu, T. Ishida, T. Koyama, M. Machida, and M. Kato, J. Phys. Soc. Jpn. 79, 124704 (2010).
- (70) Ben Xu, M. V. Milošević, Shi-Hsin Lin, F. M. Peeters, and B. Janko, Phys. Rev. Lett. 107, 057002 (2011).
- (71) L. Kramer and W. Pesch, Z. Phys. 269, 59 (1974).