Thermally active nanoparticle clusters enslaved by engineered domain wall traps

The stable assembly of fluctuating nanoparticle clusters on a surface represents a technological challenge of widespread interest for both fundamental and applied research. Here we demonstrate a technique to stably confine in two dimensions clusters of interacting nanoparticles via size-tunable, virtual magnetic traps. We use cylindrical Bloch walls arranged to form a triangular lattice of ferromagnetic domains within an epitaxially grown ferrite garnet film. At each domain, the magnetic stray field generates an effective harmonic potential with a field tunable stiffness. The experiments are combined with theory to show that the magnetic confinement is effectively harmonic and pairwise interactions are of dipolar nature, leading to central, strictly repulsive forces. For clusters of magnetic nanoparticles, the stationary collective states arise from the competition between repulsion, confinement and the tendency to fill the central potential well. Using a numerical simulation model as a quantitative map between the experiments and theory we explore the field-induced crystallization process for larger clusters and unveil the existence of three different dynamical regimes. The present method provides a model platform for investigations of the collective phenomena emerging when strongly confined nanoparticle clusters are forced to move in an idealized, harmonic-like potential.

S table localization of nanoparticle clusters on a surface may lead to several technological developments in different research fields, ranging from drug delivery 1,2 , to microfluidics 3 , optics 4 , and photonics 5 . Understanding the adsorption and mobility of nanoparticles on a substrate is also crucial in different surface-based technologies such as friction 6,7 and heterogeneous catalysis 8 . For example, localized impurities as inorganic fullerene-like nanoparticles showed excellent lubricating properties on layered materials, with the enhancement of sliding friction and exfoliationmaterial transfer 9 . In addition, investigating the dynamics of confined nanoparticle clusters is also important from a theoretical point of view. These clusters represent a simple yet nontrivial model system for studying multi-body effects in condensed matter physics, and understanding their mesoscopic dynamics may shed light on other systems on different length-and timescales. Examples span from the arrangements of electrons in a parabolic potential, where they form a Wigner crystal 10 , to the dynamics of ions in a lateral optical confinement 11 .
In colloidal science, the formation of clusters of microscale particles via isotropic, attractive interactions and their complex transition pathways have been investigated in two 12 and in three dimensions 13,14 . More recently, the effect of active perturbations on the assembly pathways has been reported with levitated granular particles assembled by acoustic forces. 15 . These works demonstrated that current experimental techniques allow manipulation and control of cluster at the colloidal length scale and above. However, applying similar approaches to investigate clusters of magnetic nanoparticles remains a challenging task. At smaller length scales, large field gradients are required to overcome thermal fluctuations and provide a stable trapping site for nanoparticles. Some successful approaches have been recently proposed based on the use of plasmonic landscapes 16,17 or hardwall confinements 18,19 . Most of such techniques require the use of lithographic patterns composed of fixed microstructures which may influence the dynamics due to steric interactions and cannot be easily changed by external control.
Here we demonstrate the stable trapping and control of thermally active clusters of magnetic nanoparticles in solution by using extended circular traps made of magnetic Bloch walls. These domain walls generate strong and tunable magnetic gradients which induce the assembly of nanoparticles into fluctuating clusters in two dimensions (2D). The research on controlled motion of domain walls in magnetic thin films is currently pushing the limit of magnetic data storage technology and is also providing applications in logic devices 20,21 , spintronics 22 , nanowires 23,24 , and ultracold atoms 25 . We use these nanoscale entities to trap and control soft magnetic nanoparticles, as an alternative approach to optical tweezers 26 or dielectrophoretic traps 27 .

Results
Experimental system. We realize cylindrical ferromagnetic domains, or magnetic "bubbles", by using a single crystal, uniaxial ferrite garnet film (FGF), see Fig. 1a-c. The FGF is grown via dipping liquid phase epitaxy on a 〈111〉 oriented single crystal gadolinium gallium garnet (Gd 3 Ga 5 O 12 ) 28 , and has the composition Y 2.5 Bi 0.5 Fe 5−q Ga q O 12 (q = 0.5−1). The grown film displays a labyrinth pattern of stripe domains with alternating perpendicular magnetization vector. The domains are separated bỹ 20 nm thick Bloch walls (BWs) where the magnetization vector rotates by 180°in the (x, z) plane. To transform the stripe pattern into a triangular lattice of magnetic bubbles, we anneal the lattice by keeping it for $ 15 min to a high frequency (0.5 kHz) magnetic field of amplitude B z = 3 mT. After switching off the field, the FGF displays a triangular lattice of cylindrical ferromagnetic domains with uniform magnetization, lattice constant a = 11.8 μm and diameter D = 8.8 μm. The size of the magnetic bubbles can be tuned by an external field perpendicular to the film, B ext ¼ B zẑ . When the field is parallel (antiparallel) to the bubble magnetization, it increases (decreases) the radius of the cylindrical domains. From the analysis of the experimental data, Fig. 1d, we extract a saturation magnetization of B s = 21.3 mT (critical field B c = 11.4 mT 29 ). More technical details can be found in "Methods". We also note that recently a strong edge stray field has been generated by an array of regular nanorod assembled in a triangular lattice 30 . However, the presence of the oppositely magnetized surrounding film impedes the formation of localized vortices in our system. Nanoparticle trapping. Above the magnetic lattice we deposit a water dispersion of paramagnetic nanoparticles with diameter d = 360 nm (Microparticles GmbH), and doped with~40% by weight of iron oxide. Before placing the particles, we use soft lithography to cover the FGF with a 1 μm thin layer of a polymer film to avoid sticking due to magnetic attraction, see "Methods" for more details. Due to density mismatch, the nanoparticles sediment in water reaching the FGF surface. As shown in Fig. 1c, see also Video S1 in the Supporting Information, the particles are effectively 2D confined within the cylindrical domains, and can only exchange position by moving within the (x, y) plane, but not along the perpendicular direction.
The BW trapping at the center of the magnetic domains could be explained by calculating the magnetostatic potential generated by the stray field of the film, B stray . The energy of one paramagnetic particle at an elevation z above such lattice is given by U 1 (r) = −υχB 2 (r)/(2μ 0 ), see "Methods". Here, B = B ext + B stray is the total magnetic field, υ the particle volume, χ~2 the effective magnetic volume susceptibility 31 , and μ 0 = 4π × 10 −7 H/ m is the magnetic permeability of the medium (water). To determine the energy landscape U 1 (x, y), we calculate numerically B stray by summing up a two-dimensional triangular lattice of cylindrical ferromagnetic domains, see "Methods" for more details. Figure 1e shows the result for zero external field (B z = 0) at the particle elevation of z = 1.3 μm. The energy landscape displays a triangular lattice of radially symmetric potential wells that are centered at the locations of the magnetic bubbles. Further, from the calculations we find that when the applied field is parallel to the bubble magnetization (B z > 0), these potential wells become deeper, while these minima disappear when the field is antiparallel to the bubble magnetization (B z < 0); we will come back to this effect later.
Single-particle fluctuations. We proceed to determine the effective shape of the confining potential from the single-particle fluctuations. As shown in Fig. 2a, a trapped nanoparticle performs a confined, angle-independent diffusive motion relative to the center of the magnetic bubble at the origin. The external field B z increases the bubble area (red circle), and makes the potential stiffer, which induces a stronger confinement (red trajectory). This effect is confirmed in Fig. 2b, where from the particle trajectory we extract the effective confining potential U 1 (r), applying the standard Boltzmann distribution 32 . Even if the real trapping potential displays a complex shape (see Fig. 1e), we find that the nanoparticle explores only the central well staying away from the boundaries where the Bloch walls are located. In such region, the potential well can be well approximated with a harmonic function, ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-25931-7   where k e is the effective spring constant. We note that exactly such radial dependence comes out naturally from our theoretical model when considering the close proximity of a magnetic bubble, see "Methods". Further, Fig. 2b shows that increasing the applied field makes the potential effectively stiffer and thus k e grows with B z . Here, our model predicts a linearly growing function k e (B z ) (see "Methods"), which is in perfect agreement with the experiment. The corresponding experimentally measured values of the spring constant, extracted from different potential distributions, are shown in Fig. 2c. As expected, the potential stiffness varies linearly with the field amplitude, and we measure a maximum value of k e = 0.034 pN/μm for B z = 4.7 mT, which corresponds to a bubble of diameter D = 9.5 μm.
The particle dynamics can be further characterized by the mean-square displacement (MSD) given by, hΔr 2 iðtÞ ¼ hðrðt þ t 0 Þ À rðt 0 ÞÞ 2 i, where r(t) = (x(t), y(t)) is the particle coordinate measured relative to the center of the bubble, t is the lag time, and 〈. . . 〉 a statistical average, which we performed over~15 independent experiments. The longtime behavior of the MSD can be described by a power law, 〈Δr 2 〉(t)~t α with an exponent α which is used to distinguish between normal (α = 1) and anomalous (α ≠ 1) diffusion. In a harmonic trap, the particle dynamics is described by an overdamped Langevin equation, ζ _ r ¼ Àk e r þ ξðtÞ, which follows directly from our model, see "Methods". Here, ζ is the particle friction coefficient, and ξ is a random force that describes thermal fluctuations. This equation is used to arrive at the 2D time-dependent MSD 33 , which describes a crossover at t ≃ τ = ζ/2k e from free diffusion at small times (t ≪ τ), 04 μm 2 /s to confined Brownian motion, 〈Δr 2 〉 (t) = 2k B T/k e at long times (t ≫ τ). The long-time limit of this equation provides an independent way to measure the effective spring constants.
When using Eq. (2) to plot all the MSD curves, Fig. 2d, we detect slight differences in the values of spring constants k e compared to those determined from the potentials in Fig. 2c. This is not unexpected because the nanoparticle subjected to thermal noise explores the weakly anharmonic nature of the genuine confining potential and small deviations relative to the harmonic model are captured unequally well by the different approaches. As shown in "Methods", the account of weak anharmonicity of our magnetic trap keeps valid the harmonic model, Eq. (1), with an effectively suppressed spring constant, k e − Δk e , where the small correction Δk e ∝ k B T > 0. Although accounting for such slight shift in k e allows us to achieve a better quantitative agreement for MSDs shown in Fig. 2d, it remains otherwise unessential for our study and therefore we neglect Δk e further.
Pair interaction and collective states. The circular magnetic traps can be used to enclose pairs of interacting nanoparticles and measure the corresponding pairwise interaction potential U 2 from the distribution statistics. As predicted by our theoretical model (see "Methods"), such interactions arise from the magnetic moments induced in paramagnetic nanoparticles by both the stray field of the FGF B stray % 3e Àκz B sẑ and external field B ext ¼ B zẑ , where κ ∝ a −1 is a damping exponent. Note that the total field B = B stray + B ext and hence the induced dipoles remain always perpendicular to the confining plane, for any external field B z , as shown in the schematic in Fig. 3a, including the case B z = 0 mT. This naturally suggests that the pairwise dipolar interactions are central and strictly repulsive, with a field-dependent strength γ = γ(B z ). We show the validity of this hypothesis first for the case B z = 0 mT, see Fig. 3a. In the presence of the external potential, the distribution of particle positions is dictated by the balance between the external and the repulsive dipolar forces. For a pair of particles with the displacements r 1 and r 2 from the trap center and interparticle distance r, we measure their total potential U = U 1 (r 1 ) + U 1 (r 2 ) + U 2 (r). A good analytic approximation for the total potential is given by U(r) ≈ 2U 1 (r/2) + U 2 (r), which is compared with U(r) evaluated directly from the distribution of interparticle distances; more details are given in "Methods". We note that the experimental data, approximate formula and results of simulations are in good quantitative agreement, see Fig. 3a. However, as expected U(r) computed from the particle distribution exhibits a slightly better correspondence to the experimental data than the prediction from the analytic formula with the same fitting parameters. Therefore, we stick to the more accurate representation when setting the mapping between the model and the experiments. Application of the external field, B z , induces a visible narrowing of the potential well: notice not only the right part of the graphs in Fig. 3b, but also steeper repulsion tails with the growth in B z , as can be observed from the left part of the graphs. Fitting the simulation data to the experiment for U(r) fixes the field-dependent strength γ(B z ) of the pairwise repulsive interactions, U 2 (r).
When increasing the number of nanoparticles in the magnetic trap, the repulsive colloids are forced to coexist and they compete for filling the central minimum. Since the confining potential and interparticle interaction potential are now completely specified by Eqs. (1) and (3) with the parameters extracted from the experiment, we construct a Brownian dynamics simulation model. Combining these deterministic forces with irregular forces caused by thermal fluctuations, we arrive at the Langevin equations for N interacting particles (see "Methods"), used to rationalize the experiments and to characterize the emerging collective states. The evolution of particle i with the position is determined by the restoring, repulsion and thermal noise forces, as given by the three terms in the right-hand side of Eq. (4), respectively. We will now report collective states for an increasing number N = 1, …, 6 of particles confined to the trap, as shown in Fig. 4a for B z = 0 mT. In Fig. 4b we show the stationary spatial density distribution ρ(x, y) at different fields B z ranging from 0 to 4.7 mT. The corresponding one-dimensional displacement distributions P(δx) obtained by cutting ρ(x, y) along one direction (y = 0) are presented in Fig. 4c, as drawn from experiments (left column) and simulations (right). Both experiments and simulations show a very good agreement in describing the particle behavior within the trap. Under zero field (first column Fig. 4b) the density field remains highly localized only for a single particle, N = 1. Already for two particles, N = 2, the density distribution significantly flattens across the trapping region and each position within the magnetic trap has almost the same probability to be visited. With the further growth in N, this feature becomes more pronounced: the density distribution broadens, while staying nearly uniform in the most part of the trap. Note that although the particles repel, tending to keeping apart and filling larger regions, the repulsion strength remains relatively weak to break the quasi-uniform structure of ρ(x, y). b Two-dimensional density distribution of the particle positions projected onto the (x, y) plane. High (low) probabilities correspond to yellow (blue) regions. Scale bar at the bottom is 2 μm. c Probability distributions P(δx) of the particle displacement along the x-axis from experiments (left column) and numerical simulation (right column).
An increase in B z strengthens the repulsion between the particles and qualitatively changes the density distribution. While one particle still shows the similar Gaussian distribution but increasingly sharper at large B z , for N = 2, 3 and 4 the density distributions ρ(x, y) become significantly nonmonotonic. Indeed, for relatively stronger repulsion, the restoring trap force is no longer able to keep particles distributed uniformly about the central potential well. Instead, they are radially displaced from the trap center and localize at a certain distance from it. Thus, the density field ρ(x, y) shows an angle-independent, ring-like structure, which corresponds to the bimodal distributions P(δx) with a local minimum at the center, see Fig. 4c. This picture, however, does not preserve for the further growth in N. Starting from N = 5, the local minimum at the center flips to a global maximum, while other maxima do not disappear and the distribution slightly broadens. At N = 5 and 6, the clusters contain enough particles to stabilize one particle at the trap center by forming a stable outer ring out of the other particles, which is impossible for smaller N. The structure of density plots reflects the shell-like ordering of particles, which is determined by the balance of confining and repulsive forces. Under these forces, the number of local maxima can be understood from the number of rings that can be formed by the given number of particles 34 .
Cluster crystallization and trap escape. Apart from contrasting the experimental data for a few trapped particles, our numerical mapping model, Eq. (4), allows exploring the assembly and dynamics for large field amplitudes, which are currently unreachable by our experimental setup. These results, however, could be readily tested with other ferromagnetic thin films able to support larger field modulations 35 . Increasing the applied field, leads to a stronger localization of nanoparticles in the trap and simultaneously strengthens the repulsion between them. In a related context, melting of confined few-body systems has been investigated in experiments 36,37 and numerical simulations [38][39][40] for dipolar particles under hard-wall confinements. However, in contrast to these works, our system is characterized by a nearly harmonic confining potential that can be controlled together with the pairwise interaction. Figure 5 shows the field-induced crystallization process for a cluster composed of N = 29 nanoparticles. We first analyze the structural transition from disorder to order in terms of the timeaveraged bond-orientational parameter Ψ 6 , see "Methods". Generally, its values range from 0 to 1, indicating complete disorder and perfect order, respectively. As shown in Fig. 5a, we can identify three different regimes. For low fields, B z < 60 mT, the particles are thermally delocalized and can explore the whole bubble area. Increasing B z , the parameter Ψ 6 raises rapidly towards a maximum value of Ψ 6 ≈ 0.7, which is kept stable over a range of B z from 60 to 150 mT. In this regime of intermediate fields, the nanoparticles start to localize along three concentric shells, with frequent intershell jumping which, however, neither destroys the shell ordering nor suppresses the radial mobility, see insets in Fig. 5a. Finally, for high fields, B z > 150 mT, we observe a decrease in Ψ 6 similar to a re-entrant behavior. Now the pairwise interactions are so strong that they impede the intershell jumping, significantly reducing the radial mobility and forcing the particles to stay within the formed concentric shells. In this situation, the angular dynamics within the shell intensifies, while the radial motion freezes, a re-entrant effect that in many cases appears similar to the case of hard-wall confinement of particles 36 . To further quantify the crystallization process we compute another parameter u rel , the relative interparticle distance fluctuations 41 (see "Methods"). As shown in Fig. 5b, the relative interparticle distance fluctuations u rel decrease with the increase of the external field, with an amplitude of order similar to that reported for microscopic particles 40 . Further, we observe a clear change in slope of u rel close to the second transition at high field where only intrashell motion occurs.
While the experiments reported until now demonstrate the trapping within the circular Bloch walls, our system allows us also to change the particle location across the lattice and thus to exchange particles between magnetic bubbles. Thus, the trapped nanoparticles are enslaved within the magnetic bubble, and can leave them only by an external command, such as an applied field. This feature is demonstrated in Fig. 6a-c by switching the direction of the applied field B z → − B z thus, the field becomes antiparallel to the bubble magnetization. In this situation, the nanoparticles are forced to jump into the interstitial region (first and second image), while a subsequent field inversion (−B z → B z ) forces the nanoparticles to return back to one of the nearest bubbles (second and third image), see also Video S2 in the Supporting Information. As shown in Fig. 6a, the central bubble is initially filled by N = 5 particles which are then ejected toward six interstitial places, where they meet other nanoparticles coming from neighboring domains. Thus, each interstitial region receives a fraction of the nanoparticles from the three nearest neighbor bubbles, and after the second field reversal only N = 3 have returned back to the original bubble. The exchange process can be understood by analyzing the magnetostatic potential in the particle plane (x, y), see Fig. 6c. When the field is antiparallel to the bubble magnetization, it induces a reversal of the energy landscape and the bubbles correspond not to minima but maxima of the potential, while six potential wells of triangular shape are nucleated around each bubble. This induces a particle displacement at a speed of v~20 μm/s, which is proportional to the local magnetic gradient. As shown in the third image of Fig. 6a, some particles may stay within the interstitial region even after the second field reversal. They are located in an unstable equilibrium, and can easily jump back to a bubble due to thermal fluctuations. Reversing the field polarity may be used to remove nanoparticles from each bubble, which later can be refilled with another field reverse. Thus our virtual magnetic trap may be also used as a rudimentary form of logic memory based on magnetic nanoparticles, which can be stored within the platform by localizing them, and later "erasing" their location by an external command. In particular, a recent work 42 demonstrated the possibility of combining lithographic confinement, electrostatic levitation and external actuation to store and retrieve logic information from levitated nanoparticles. Since our FGF could be easily coupled with a lithographic structure, similar operations could be parallelized on the whole bubble lattice using an external magnetic field.

Discussion
We demonstrate the controlled magnetic trap to assemble and manipulate clusters of magnetic nanoparticles in solution. We combine experiments with theory to show that the equilibrium dynamics of the clusters arises from the balance between the confining harmonic potential and repulsive dipolar interactions induced by the substrate and the applied field. The external field serves as the only control knob that may be used to tune both the stiffness of confinement and the strength of interparticle interactions. We further construct a numerical simulation model based on the experimental parameters as a mapping between the experiment and theory, which we exploit to extrapolate the experiment for larger ranges of fields and particle numbers. We investigate the field-induced crystallization process and reveal three different dynamics regimes. Finally, we experimentally demonstrate delocalized transport of nanoparticles across the lattice and the exchange process between the magnetic domains by inverting the applied magnetic field. We note also that trapping of microscale colloids has been achieved both by lithographic wedges 43,44 and high-frequency electric field 45 . However, our work allows us to trap nanoscale magnetic particles via a magnetic landscape with tunable stiffness, keeping them confined to a surface in stable two-dimensional clusters in absence of any topographic relief.
Magnetic nanoparticles are the focus of intense scientific research for their special physical properties which make them widely used in biomedicine 46 , microfluidics 47,48 , magnetic imaging 49,50 or nanomaterial-based catalysts 51 , among others. However, the doping with iron oxide makes these nano-size objects light adsorbing, and thus rather difficult to be reliably trapped via optical tweezers. Although recent progress in this direction has been demonstrated at the single-particle level 52,53 , extensions to stable entrapment of ensembles/clusters of nanoparticles in more than one spatial direction remains challenging. With our study, we are able to overcome these restrictions because our setup provides a versatile platform to localize, stably confine and tune the behavior of nanoparticles on a spatially extended surface.
The magnetic bubble lattice can be easily configured by an out of plane magnetic field and the corresponding bubble domains manipulated via a small gradient when the field is still on. Further, given the fast response of the Bloch wall to magnetic field, with a propagation speed of the order of~1km s −1 54 , these features make this type of approach a promising candidate for magnetic-based nanoscale trapping. While our magnetic domains are relatively large, much smaller bubbles could be also synthesized in garnet film 55 that could eventually lead to further system miniaturization. Moreover, the proposed technique can also be extended to biological systems such as magnetotactic bacteria 56 and different dispersing media such as viscoelastic fluids 57 . In the latter case, plans to use our trapped nanoparticles as strongly thermal microrheological probes are currently under research.

Experimental system
Preparation of the ferrite garnet film (FGF). We coat the FGF with a thin polymer film (1 μm) made of a positive photoresist AZ-1512 (Microchem Newton, MA). The coating increases the effective particle elevation from the FGF, which has a twofold effect. First, it makes the potential effectively parabolic. Second, it reduces the otherwise strong attraction by the Bloch walls, preventing the particle sticking to the FGF. The photoresist is applied via spin coating (Spinner Ws-650Sz, Laurell) and subsequent UV photo-crosslinking (Mask Aligner MJB4, SUSS Microtec). We use the following procedure. First, we clean the FGF in an ultrasonic bath filled with acetone (Merck) for $ 15 min and then with isopropanol (Sigma) for the same amount of time. After that, the FGF is dried under a stream of N 2 . Few drops of the photoresist are deposited above the cleaned FGF, and then dispersed using a spin coater working at 3000 rpm for 30 s. After that, the polymer is baked by placing FGF above a plate heated at 95 ∘ C for 1 min, and then photo-crosslinked via exposure to 5 s of UV irradiation at a power of P = 30 mW/cm 2 . A final post-bake process is applied by placing the FGF above a hotplate heated at 115 ∘ C for 50 s. The latter process induces further hardening of the photoresist above the FGF.
Preparation of colloidal solution. We dilute a drop of the stock solution of the nanoparticles in highly deionized water (MilliQ, Millipore). The particles become electrostatically stabilized in water by the negative charges acquired from the dissociation of the surface carboxylic groups (COOH). We tune the amount of H 2 O to reach a density of~10 7 beads/ml and add few drops above the magnetic film. After a few minutes, the particles sediment to the FGF surface due to density mismatch and they remain confined above it without sticking. We wait for 15 min to allow the sedimentation of all particles, place above a coverslip (no.1, Thermo Scientific Menzel) and use an immersion oil (Immerso 1111-806, Zeiss) between the coverslip and the microscope objective (100 × 1.3NA, Nikon).
Magnetic field and particle tracking. The external magnetic field perpendicular to the FGF is applied by using a custom-made coil placed below the magnetic film (see Fig. 1a). The coil is composed of 700 turns of a 0.5 mm thick copper wire and is connected to a DC power supplier (TTi El 302), which allows to apply a spatially uniform, static magnetic field up to B z ≃ 20 mT. We use a teslameter (FM 205, Projekt Elektronik GmbH) to calibrate the field amplitude and determine the homogeneity of the field distribution around the sample plane. We find that for the amplitudes used (B z ≤ 6 mT) the field is spatially uniform above the sample area of 1 cm 2 . The particle positions are determined via digital video microscopy 58 . We use an upright optical microscope (Eclipse Ni, Nikon) equipped with a 100 × 1.3 NA oil immersion objective and a CCD camera (Basler Scout scA640-74fc, Basler) to record experimental movies at 75 frames per second.

Theoretical model
Derivation of the model General framework. The experimental measurements are interpreted within a model that describes overdamped motion of an ensemble of i = 1, …, N nanoparticles with positions r i = (x i , y i ) in a potential U and subject to thermal fluctuations, where ζ is the friction coefficient, ξ is the Gaussian white noise with zero mean and diagonal covariance of strength 2ζk B T with k B T being the thermal energy. The potential U(r i , r j ) = U 1 (r i ) + ∑ j>i U 2 (r i , r j ), comprising the single-particle interaction with the external field, U 1 , and pairwise interactions with all other particles, U 2 . A spherical paramagnetic particle of diameter d in a magnetic field H behaves as an induced magnetic moment m = υχH, where χ is the effective magnetic susceptibility and υ = πd 3 /6 is the volume of particle. Therefore, considering U 1 (r i ) = − (m i ⋅ B i )/2 with B i = B(r i ) as the energy of the dipole-field interaction and U 2 as the energy of dipole-dipole interactions, we have where r ij = r i − r j ,r ij ¼ r ij =r ij and r ij = |r ij |. Note that because the solvent is nonmagnetic, the magnetic induction B = μ 0 H with μ 0 = 4π × 10 −7 H m −1 the magnetic permeability of free space.
Exact stray field above the FGF. The stray field generated at the surface of the magnetic bubble lattice, B stray , can be calculated exactly by summing up the field from a triangular lattice with period a of magnetic bubbles with lattice vectors p = na − + ma + with a ± :¼ ð ffiffi ffi 3 p a=2; ± a=2Þ, and n, m integers. Each bubble is considered as a cylindrical uniformly magnetized ferromagnetic domain, generating a stray field above its surface written in cylindrical (r, z) coordinates as 59 : b ¼rb r ðr; z; tÞ þẑb z ðr; z; tÞ with Here, B s is the saturation magnetization of the FGF, r ¼ ðx À p x Þx þ ðy À p y Þŷ, Q n is the Legendre function of the second kind and order n, Π(n, m) gives the complete elliptic integral of the third kind, K ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2rD=½z 2 þ ðr þ D=2Þ 2 q ; n ± ¼ 2r=ðr ± ffiffiffiffiffiffiffiffiffiffiffiffiffiffi r 2 þ z 2 p Þ, the bubble diameter depends on the external field B ext and can be expressed as: q and The field generated by such array of bubbles is given by B b = ∑ n,m b nm with the indexes n and m over the entire triangular lattice. The overall substrate field is obtained as B stray ≔ B b − B f , where B f is the contribution due to the oppositely magnetized film calculated for a cylindrical domain covering the entire sample area. This approach is used to calculate the landscapes in Figs. 1e and 6c.
Reduced description. The total magnetic field above the substrate B = B ext + B stray , where the external field B ext ¼ B zẑ and an approximate solution for the stray field that captures triangular symmetry of the magnetic bubble lattice (see supplementary information in Ref. 60 ) can be represented as B stray (x, y, z) κy=2Þ. Here, z is the particle elevation above the substrate, and κ ¼ 4π=ð ffiffi ffi 3 p aÞ with the lattice constant a. Expanding the solution about the center of a magnetic bubble, κr ≪ 1, we arrive at the reduced analytic expession with r 2 = x 2 + y 2 . Employing this result in Eq. (6), we obtain a harmonic approximation for the confining potential, U 1 (r) = k e r 2 /2 (Eq. (1)), with a fielddependent spring constant of the form, k e (B z ) = k 0 (1 + c 1 B z ), where k 0 and c 1 are fitting constants. To evaluate the pairwise interaction potential, Eq. (7), we stick to the limit κr → 0 in Eq. (9), which reveals central strictly repulsive interactions between the particles, U 2 (r) = γ/r 3 (Eq. (3)), with the field-depenent repulsion strength γðB z Þ ¼ γ 0 ð1 þ c 2 B z Þ 2 and fitting constants γ 0 and c 2 . Accounting for these results in Eq. (5), we obtain the governing equation, see Eq. (4), serving as a mapping between the theory and experiment.
The originally anharmonic nature of confining potential may lead to deviations in effective spring constants measured from the particle distribution P 1 (r) and from the MSD. For the latter, our measurements indeed show slightly smaller values of k e , a known effect that is well captured by small anharmonic corrections to the potential, see Appendix B 61 ; for limited statistics, they may not be well detectable from measuring P 1 (r) because of its sharp exponential form. Within our model, they are equivalent to a coordinate-dependent spring constant k 0 e ðrÞ ¼ k e ð1 À c 4 r 2 =4Þ. Evaluating its average, we remain within the harmonic approximation for U 1 with a slightly modified spring constant hk 0 e i ¼ R k 0 e ðrÞP 1 ðrÞdr ¼ k e À Δk e ; Δk e ¼ c 4 k B T=2 > 0, in quantitative accord with our MSD data. Because the constant correction Δk e is otherwise irrelevant for our study, we do not take it into account.
Pairwise interactions. The parameters of pairwise repulsion interaction are evaluated from the case of two particles, for which the total potential U(r 1 , r 2 ) = U 1 (r 1 ) + U 1 (r 2 ) + U 2 (r). Here, r = |r 1 − r 2 | is the interparticle distance. The corresponding two-particle stationary probability distribution P 2 ðr 1 ; r 2 Þ / exp½ÀβUðr 1 ; r 2 Þ can be projected to P 2 ðrÞ ¼ R P 2 ðr 1 ; r 2 Þδðjr 1 À r 2 j À rÞdr 1 dr 2 / exp½ÀβUðrÞ. Therefore, we evaluate U(r) as Àln ½P 2 ðrÞ=P max by measuring P 2 (r), where P max ¼ max r P 2 ðrÞ; with this normalization U min ¼ min r UðrÞ ¼ 0. We match the corresponding data from the experiment and simulations, see Fig. 3, by keeping the trap parameters as described above and tuning the strength γðB z Þ ¼ γ 0 ð1 þ c 2 B z Þ 2 of repulsive interactions, U 2 (r) = γ/r 3 . We find the best fit for βγ 0 ≈ 2.4 μm 3 for B z = 0 mT and c 2 = 0.07 mT −1 , which accounts for the field dependence.
For a system with an elongated, quasi-one-dimensional trap, for which the most probable configuration corresponds to r 1 = − r 2 = (r/2, 0) with the trap center at origin, and a central pairwise potential U 2 (r) a reliable analytic approximation is available 62 , U(r) ≈ 2U 1 (r/2) + U 2 (r). For our two-dimensional trap, this estimate is still valid, because due to axial symmetry, cf. Fig. 4b, the most probable configuration can be similarly represented as r 1 = − r 2 = r/2, leading to the approximation P 2 ðrÞ / exp½Àβð2U 1 ðr=2Þ þ U 2 ðrÞÞ. The relative weights of other possible configurations seem weak to introduce significant errors, and the approximate formula U(r) ≈ 2U 1 (r/2) + U 2 (r) remains a robust reference also in two dimensions, see Fig. 3a.
Analysis of collective states. To study collective states, cluster crystallization and trap escape, we intergrate Eq. (4) numerically. The parameters of the trap and repulsion, including their field dependence, are taken as determined in the singleand two-particle setups. We have additionally ensured that every particle interacts with all other particles. This way, we have access to particle positions within the numerical model. Experimentally, particle positions are extracted by particle tracking. These data are used to evaluate and compare such quantities as meansquare displacement and particle distributions.
To quantify the degree of local hexagonal ordering, from particle positions we also calculate the bond-orientational order parameter, defined as where N b is the number of neighbors of particle k, and θ kj is the angle between a fixed axis and the bond joining particles k and j. Further we average it at each time step over all particles. Following a previous work on optically trapped colloidal cluster 63 , to minimize artefacts from the curved frontier we calculate ψ 6 (t) = 〈ψ 6,k 〉, by considering only the particles not adjacent to the outer boundary. Finally, we perform a running-time average and obtain the measure Ψ 6 ¼ ψ 6 ðtÞ for an entire simulation, as shown in Fig. 5. Another order parameter that we use to quantify the crystallization process is the relative interparticle distance fluctuations 41 , u rel ¼ 2=ðNðN À 1ÞÞ ∑ N 1 ≤ i<j ½hr 2 ij i=hr ij i 2 À 1 1=2 , where r ij is the distance between particles i and j, N is the total number of particles, and we have also performed the running-time averaging.

Data availability
The data that support the findings in this study are available within the article and its Supplementary information. Further details are available from the corresponding authors upon reasonable request.