Description of the simulations :
Bound state configurations of $N_v$ vortices found in high-precision large-scale numerical minimzation of the two-component Ginzburg-Landau free energy in the type-1.5 regime. The variational problem is defined using a finite element formulation provided by Freefem++ framework. A nonlinear conjugate gradient algorithm is used to find the minima of the energy. More details about the numerical methods employed to study can be found in a Supplementary Material Letter .

In the following, we present evolution from initial configurations being either a giant vortex, or a diluted collection of vortices.

The displayed quantities are: From left to right, on the first line are plotted the magnetic field, and condensate densities $|\psi_1|^2$ and $|\psi_2|^2$. On the second line, the supercurrents moduli $|J_1|$ and $|J_2|$ and $~\mathrm{Im}(\psi_1^*\psi_2)$ which is nonzero when there appears the phase difference between both components.
1. Superconductors with two active bands:
We first consider the case where both components are below their critical temperature (i.e. they are both active). The first band is the stronger while the second one is weaker.

Movie 1: The ground state of an $N_v=9$ flux quanta configuration in type-1.5, $U(1)\times U(1)$ superconductor (i.e. $\eta=0$).The parameters of the potential are a1b1 $(\alpha_1,\beta_1)=(-1.00,1.00)$ and $(\alpha_2,\beta_2)=(-0.60,1.00)$ , while the electric charge is $e=1.48$ (corresponds to Figure 2 of the paper).

Movie 2: The ground state of a $N_v=12$ flux quanta configuration using the same parameter set as in Movie 1. Here the global 8-folded discrete symmetry of the cluster has been broken toward a less symmetric configuration (a 3-folded discrete symmetry) in favor of a regular vortex lattice in the cluster. There is again a competition between type-I-like (normal circular cluster with a boundary current) and type-II-like tendencies (vortex lattice). This corresponds to Figure 3 of the paper.

Movie 3: Magnetic ground state of an $N_v=9$ vortex configuration with a stronger electric charge coupling $e=1.55$ (parameters of the potential are the same as in Movie 1). This shows the behavior of the system with respect to a charge increase. Increasing the electric charge decreases penetration length and thus pushes the system towards type-I regime i.e. with more circular boundary (corresponds to Figure 4 of the paper).

Movie 4: Magnetic ground state of an $N_v=9$ vortex configuration with the parameter set given by Movie 1, but with $e=1.59$ and added Josephson coupling $\eta=0.1$. Although the Josephson term introduced an energy penalty for phase difference, it has little effect on the vortex cluster boundary where the high magnetic pressure still generates strong phase difference gradients. This corresponds to Figure 3 of the paper.

Movie 5: Magnetic ground state of an $N_v=12$ vortex configuration. Parameters of the potential are $(\alpha_1,\beta_1)=(-1.00,1.00)$ and $(\alpha_2,\beta_2)=(-0.0625,0.25)$, while the electric charge is $e=1.30$ and added Josephson coupling $\eta=0.5$. Although the Josephson term introduced an energy penalty for phase difference, it has little effect on the vortex cluster boundary where the high magnetic pressure still generates strong phase difference gradients.

Movie 6: Magnetic ground state of an $N_v=18$ vortex configuration. Parameters of the potential are $(\alpha_1,\beta_1)=(-1.00,1.00)$ and $(\alpha_2,\beta_2)=(-0.0625,0.25)$, while the electric charge is $e=1.30$ and added Josephson coupling $\eta=0.5$. Although the Josephson term introduced an energy penalty for phase difference, it has little effect on the vortex cluster boundary where the high magnetic pressure still generates strong phase difference gradients.

2. Superconductor with strong interband Josephson coupling:
Here we consider the case where one band is below its critical temperature (i.e. active) while the second band has superfluid density only because interband proximity effect (i.e. passive).

Movie 7: A bound state of an $N_v=25$ vortex configuration in case when superconductivity in the second band is due to interband proximity effect and the Josephson coupling is strong $\eta=7$. $(\alpha_1,\beta_1)=(-1.00,1.00)$ and $(\alpha_2,\beta_2)=(3.0,0.5)$ and $e=1.30$. This corresponds to Figure 6 of the paper.
Simulations show the increased importance of non-pairwise interactions which lead to formation of vortex chains rather than a compact cluster.

Movie 8: A bound state of an $N_v=25$ vortex configuration in case when superconductivity in the second band is due to interband proximity effect and the Josephson coupling is strong $\eta=7$. $(\alpha_1,\beta_1)=(-1.00,1.00)$ and $(\alpha_2,\beta_2)=(3.0,0.5)$ and $e=1.30$. This corresponds to Figure 7 of the paper. We start with a diluted initial configuration where vortices are located with the attractive range of their two-body potential. However instead of contracting, initially the vortex cluster expands, thus showing that the evolution is dominated by nonpairwise interaction. Once the group breaks into several subgroups of vortices and thus gets rid of some of the multibody forces, the subclusters start to shrink. The intervortex distance within each subcluster in the final configuration is smaller than that in the initial configuration despite the expansion in the initial stage of the evolution.