50 40 30 20 10 0 50 40 30 20 10 0 Magnetic field ( T) (a) (b) (c) (d) FIG. 4: Simulations of systems consisting of (a) 6, (b) 27, (c) 122, and (d) 500 vortices which have the interaction en- ergy calculated from a two-gap GL theory for the type-1.5 condition (see text). A scale bar corresponds to 10 µ m. FIG. 3: Scanning SQUID microscope images of vortices at 20 µ T obtained from different cool downs. The left image is the same as Fig. 2(b). 3 by 46 and 116 estimated from SB/φ$_{0}$ , respectively. We note that the vortex array varies considerably in the locations of clusters measured by different cool downs, while keeping the same general features, as seen in Fig 3. Since the pinning sites are fixed, this implies that the vortex clustering is not due to the presence of pinning centers but is rather to be attributed to the pe- culiar vortex-vortex interaction combining repulsion and attraction. In a Nb film with the strong pinning, disordered vortex patterns are observed at low vortex densities, as shown in Fig. 1(d). We now compare qualitatively these vortex patterns. Figure 1(e) shows the 2D multi-fits 23 of the scanning SQUID image in Fig. 1(b) to numerical cal- culations from a London model $^{20,21}$. Fits were made by parameterizing the coordinates of each position of vor- tices with z and λ determined above. From the fits it turns out that the number of vortices is 23, which is in- deed consistent with 23 estimated by SB/φ$_{0}$ , where S is the areas of a sample. The distribution of the intervortex distance can be derived from the fits, as shown in Fig. 1(f). This distribution is a bimodal and it clearly splits into two: one ranging from ∼ 2 µ m to ∼ 7 µ m and another from ∼ 15 µ m to ∼ 30 µ m. In principle, the former corre- sponds to the intragroup distribution and the latter cor- responds to that concerning the intergroup. On the other hand, the vortex distribution in a Nb film can be fitted by the Gaussian function, as shown in Fig. 1(g). The cen- ter of the peak in the Gaussian function is at 14.5 µ m, indicating that the average intervortex distance is 14.5 µ m. This agrees with the intervortex distance ( φ$_{0}$/B ) 1 / 2 of ∼ 14 µ m estimated by assuming the presence of the Abrikosov lattice. This fact indicates that the vortex pattern in Fig. 1(b) is not a disordered Abrikosov lattice but vortex cluster array. On the other hand, the vor- tex pattern in Fig. 1(d) does correspond to the strongly disordered Abrikosov lattice. Figures 4(a)-4(d) show results of the simulation men- tioned above. In type-1.5 superconductors vortices begin to form clusters consisting of a few vortices at very low vortex densities, [Fig. 4(a)] which does not coincide with the experiment displayed in Fig. 1(a). Perhaps, the exis- tence of some pinning centers may be responsible for this fact. This will be considered in the future work. The number of small clusters increases with applied fields, [Fig. 4(b)] which apparently corresponds to Fig. 2(g). Then small clusters start to merge with each other,[Fig. 4(c)] to grow into bigger clusters with a further increase of an applied field [Fig. 4(d)]. This simulation result co- incides quite well with experimental data [Figs. 2(h) and 2(i)]. The simulations show that in a vortex cluster the dis- tances between vortices and their nearest neighbors are almost independent of the applied field ( ∼ 1 µ m). There- fore, these distances are not determined by the expres- sion ( φ$_{0}$/B ) 1 / 2 but rather by the competition between short-range repulsive and long-range attractive interac- tions between vortices. In Fig. 5 the average intervortex distance in clusters obtained from scanning SQUID and decoration images 5 is shown as a function of applied field, compared with the dependence ( φ$_{0}$/B ) 1 / $^{2}$. The average distances indeed do not follow any dependence close to ( φ$_{0}$/B ) 1 / 2 but follow instead a very weak linear field de- pendence. The intervortex distance in clusters, indepen- dent of the applied field, can be identified as a unique property of a system of vortices with short-range repul- sion and long-range attraction, i.e., vortices induced in a type-1.5 superconductor. This fact indicates that vor- tex clusters definitely emerge even in crystals grown by a different crystal maker, implying that in general vor- tices tend to form vortex clusters in high quality single crystals of MgB$_{2}$. 5 Figures 2(a)-2(c) show the variation of the vortex structure with applied fields. After vortices form local clusters in Fig. 2(a), some islands of magnetic flux and the semi-Meissner regions come out of the grouping with increasing applied field. [Figs. 2(b) and 2(c)] The posi- tions of the vortices can be determined by the 2D multi- fits [Figs. 2(d)-2(f)] of Figs 2(a)-2(c), as shown in Figs. 2(g)-2(i). The position maps show that vortices line up locally in the early stage [Fig. 2(g)] and then small vor- tex clusters merge with each other, [Fig. 2(h)] to grow into bigger clusters with applied fields. [Fig. 2(i)] The number of vortices changes to 49 in Fig. 2(h) and then 114 in Fig. 2(i). These are reasonably well approximated A difference between ( φ$_{0}$/B ) 1 / 2 and the intervortex distance in clusters decreases as an applied field increase