Revealing emergent magnetic charge in an antiferromagnet with diamond quantum magnetometry

Whirling topological textures play a key role in exotic phases of magnetic materials and are promising for logic and memory applications. In antiferromagnets, these textures exhibit enhanced stability and faster dynamics with respect to their ferromagnetic counterparts, but they are also difficult to study due to their vanishing net magnetic moment. One technique that meets the demand of highly sensitive vectorial magnetic field sensing with negligible backaction is diamond quantum magnetometry. Here we show that an archetypal antiferromagnet—haematite—hosts a rich tapestry of monopolar, dipolar and quadrupolar emergent magnetic charge distributions. The direct read-out of the previously inaccessible vorticity of an antiferromagnetic spin texture provides the crucial connection to its magnetic charge through a duality relation. Our work defines a paradigmatic class of magnetic systems to explore two-dimensional monopolar physics, and highlights the transformative role that diamond quantum magnetometry could play in exploring emergent phenomena in quantum materials.


I. INTRODUCTION
Topologically protected states in magnetic materials are promising candidates for novel spintronics architectures, such as neuromorphic computing [1,2].In particular, topological textures in antiferromagnets (AFMs) could provide additional advantages over their ferromagnetic counterparts including enhanced stability, as well as faster and richer dynamics [3][4][5][6][7][8][9][10][11][12][13][14][15][16].The vanishing net moment brings forth its own challenges, particularly making the detection of AFM textures difficult.Synchrotronbased dichroic X-ray techniques, sensitive to the staggered magnetization, are at the imaging forefront and have for the first time revealed the existence of twodimensional topological AFM spin textures in hematite, α-Fe 2 O 3 [17,18].This powerful technique, however, comes with the caveat that it is insensitive to the sign of the staggered magnetization and thus the associated vorticity, i.e. the actual whirling of the spin textures, is not observed.
When viewed instead through the lens of the canted magnetization, we uncover weak but finite magnetic fields emanating from the twists of the staggered moments.The presence of such magnetic fields can equivalently be described by the magnetic analogue of Gauss's law [19], therefore pointing to the existence of a distribution of emergent magnetic charges in a topologically-rich AFM landscape.The state of the art vectorial weak-field sensing of diamond quantum magnetometry (DQM), employing a single nitrogen-vacancy (NV) colour centre as a point field sensor [20][21][22][23][24], therefore puts it in a unique position to study the above proposed conceptual idea of emergent magnetic charges in a completely new class of magnetic materials -canted AFMs.
In this letter, we demonstrate DQM imaging of topological textures in the archetypal AFM α-Fe 2 O 3 and we show that these AFM textures host a rich tapestry of magnetic charge distribution.In particular, the duality relation between staggered vorticity and magnetic charge allows us to associate the AFM Bloch meron with a spatially extended emergent magnetic monopole.Distinct from emergent magnetic monopoles in other realizations, such as spin ice [25], we observe here that the positively and negatively charged monopolar textures are topologically equivalent, while their topological antiparticle has a magnetic quadrupolar character.In addition to being a promising tool to investigate topological AFM textures, our results demonstrate the DQM potential for discovering emergent phenomena in model magnets.

II. DETAILS OF α-FE2O3 AND DQM
Hematite α-Fe 2 O 3 is an AFM oxide insulator, which hosts a wide family of topological spin textures [17,26,27].Figure 1(a) illustrates the trigonal corundum atomic structure of α-Fe 2 O 3 .It comprises a stack of anti-parallel FM sublattices along the c-axis, with magnetisation textures M 1 (yellow) and M 2 (green), as shown in Fig. 1(b).Spin re-orientation occurs at the Morin transition temperature, T M ∼ 200 K [17](see supplementary note S1); below and above T M , the FM sublattices lie predominantly out-of-plane and in-plane, respectively.The Néel vector l = M 1 − M 2 characterises the AFM order, while m = M 1 + M 2 is the net magnetisation, as shown in Fig. 1(c).Above T M , m has a predominantly in-plane orientation with an average magnitude m ∆ ∼ 2×10 3 A/m.This weak magnetisation is due to the slight in-plane canting of angle ∆ ∼ 1.1 mrad (see S1) between M 1 and M 2 , as a consequence of the intrinsic bulk Dzyalonshinskii-Moriya interaction (DMI) vector along the c-axis [28](see S5).Consequently, m lies in-plane and satisfies the relation m • l = 0. Since m ∆ is around three orders of magnitude weaker than | l|, this canting-induced weak magnetisation has no discernible effect on the AFM character of α-Fe 2 O 3 .Finally, m vanishes in situations where l turns out-of-plane below T M or due to the local formation of AFM spin textures.
We quantify the magnetic-field distribution from these spin textures via a custom-built cryogenic DQM (see S2). Figure 1(d) illustrates our diamond cantilever probe hosting a single NV centre, which is scanned at a constant height above the sample surface.The NV centre is an atomistic single spin defect with a paramagnetic ground-state manifold and state-selective optical transitions.This allows the Zeeman splitting between the ground states |±1 and |0 to be probed with a microwave frequency f mw sweep and optical excitation, i.e. via optically detected magnetic resonance (ODMR).Consequently, in the weak-field approximation [29] with negligible strain, we infer the local B-field projected onto the NV-axis (B NV ) from the energy difference between |0 and |+1 given by ∆E + = h(|f + −D|−γB bias ) = hγB NV , where h is Planck's constant, f + is the resonant frequency corresponding to the transition, D ≈ 2.87 GHz and γ = 28 MHz/mT.A small bias field B bias ≈ 0.5 mT is applied along the NV-axis to enable extraction of the field orientation.Figures 1(e) and (f) illustrate the typical variation in ODMR frequency across a linescan over the α-Fe 2 O 3 surface at T = 4 K and 300 K, respectively.The colour plot displays the amplitude of the ODMR signal, while the white curves demarcate f + used to extract B NV .An ODMR raster scan across the sample surface provides a B NV image.We then transform to the laboratory coordinates (B x,y,z ) via the Fourier reconstruction technique [30], where z coincides with the c-axis of the crystal.Figures 1(g) and (h) present example images of B z collected at 4 K and 300 K, respectively.The qualitative differences observed in these images already reveal distinct magnetic phases for temperatures below and above T M [17].The B z image below T M is generally weak, consistent with the absence of net magnetisation.By contrast, the B z image above T M shows strong and distinct features, consistent with the expectations of nonzero net magnetisation in this phase.

EMERGENT PROPERTIES IN α-FE2O3
To gain a physical interpretation of the underlying magnetisation distribution from the measured B z images, we begin from a thin-film approximation given by (see S4.1): where * indicates convolution, α i (i = xy, z) are effective point spread functions [31], t is the film thickness, d is the height above the film surface, and m xy and m z are the in-plane and out-of-plane components of m. m xy and m z contribute to B z through the divergence and Laplacian, respectively.The α i account for the B-field decay above the surface, acting as blurring kernels with their size of order d.Hence, the spatial resolution of our DQM imaging technique is set by the effective NV-sample distance, d NV (see S3). Due to the symmetry of the DMI in α-Fe 2 O 3 (see S5), m z = 0 and m xy = 0, rendering the second term in Eq. 1 identically zero.Therefore, B z images are the divergence of the canted magnetisation ∇ • m xy , convolved with α xy .Moreover, for a ẑ-oriented DMI, the net magnetisation is also given by m = ∆(ẑ × l), where ∆ is the DMI-set canting angle.This yields an additional expression, namely ∇• m xy = ∆[ẑ •( ∇× l )] (see S5).The striking consequence of this expression is that B z images also offer a projected measure of the staggered vorticity, i.e. the curl of the Néel vector, V = ∇ × l.

CHARACTERISTIC FIELD SIGNATURES AND VORTICITY READOUT
With the gained insight on the relationship between the B z images and m xy , we show next the B z images obtained in Fig. 1 are produced by AFM antiphase domain walls (ADWs), merons, antimerons and bimerons -consistent with their recent observation in α-Fe 2 O 3 [17].Below T M , we model the B z images with a linear AFM domain wall ansatz [17], characterised by width w and phase ξ a (see S5.1).The phase ξ a controls the spatial variation of l, resulting in an AFM Néel (a-Néel) or an AFM Bloch (a-Bloch) ADW profile for (ξ a = 0, π) and (ξ a = π 2 , 3π 2 ), respectively.For a linear ADW profile centred at x = 0 along the x-axis, ∇ • m xy = m ∆ π w sin πx w sin(ξ a ) for |x| ≤ w 2 , and zero elsewhere [17] (see S5.1).In other words, we expect ADWs to display a sinusoidal profile in ∇ • m xy and hence in B z , with zero crossing at the centre and amplitude and sign modulated by sin (ξ a ).In particular, an a-Néel ADW will not yield a B z signal as ∇ • m xy = 0, whereas an a-Bloch counterpart will show maximal signal.Based on Eq. 1, these characteristics are reflected in the simulated B z image of an ADW in Fig. 2(a) (see S5.4) assuming a constant phase ξ a = π.The measured B z image in Fig. 1(g), as well as a close-up image in Fig. 2(b), capture the signature zero crossing of an ADW.DQM also reveals variations of B z along the wall boundary capturing the spatially varying nature of the phase ξ a , as expected in α-Fe 2 O 3 [17].In effect, the full phase dependence of B z allows us to identify unambiguously the zero-signal sections along the varying ξ a wall boundary as a-Néel ADWs. Figure 2(c) presents a reconstruction of the m xy distribution of the multi-chiral ADW (illustrated by arrows), obtained by fitting its B z field to the experimental data in Fig. 2(b) through systematic regularization strategies (see S7).
In contrast, DQM at 300 K captures relatively larger spatial features of strong B z signal, as showcased in Fig. 1(h).This above-T M regime is where we expect finite net magnetisation forming a tapestry of localised whirling topological structures, such as multi-chiral merons, antimerons and bimerons [17], as well as topologically-trivial in-plane domain walls.Each of these AFM textures produces a distinctive B z signal, allowing us to develop a  systematic procedure to disambiguate them, see S5.Here we focus on two-dimensional (2D) topological textures and model isolated merons and antimerons based on a linear ansatz [17,32], described by the phase ξ a and the topological winding number N (see S5.3).The magnetisation divergence for such an object in polar coordinates , where f (r) is a radial function dependent on the (anti)meron phase (see S5.3).A meron (N = +1) produces a radially symmetric B z distribution about its core with magnitude and polarity controlled by sin (ξ a ).Analogous to ADWs, a-Néel merons (ξ a = 0, π) are divergence-free and exhibit 2 ) show maximal B z amplitude.In contrast, for an antimeron (N = −1), the B z distribution is two-fold symmetric and ξ a controls a trivial azimuthal offset.Simulated B z images in Fig. 2(d, g, j) of the a-Bloch meron of both polarities and the antimeron reinforce these observations (see S5.4).Thus, DQM unambiguously reveals the topological winding number N and the staggered vorticity V for each spin texture.The final attribute of such a spin texture, the sign of its topological charge, is however rendered indistinguishable by DQM due to the vanishing canted moment at its core.Figure 2 (panels e, h, k) presents the measured B z images of an anti-clockwise a-Bloch meron (N = +1, ξ a = π/2), a clockwise a-Bloch meron (N = +1, ξ a = 3π/2), and an antimeron (N = −1), respectively, in good agreement with their simulated counterparts.Further, in Fig. 2 (f, i, l) we once again reconstruct the corresponding m xy distributions and use them to simulate the measured B z image for each texture discussed above (see S7).Given the typical density of spin textures evident in Fig. 1(h), the reconstructed m xy approach better captures the finer details of the measured B z images in the absence of true isolation.We note that the initial density of (anti)merons can be reduced via meron-antimeron annihilation mediated by an external in-plane magnetic field [17] (see S8). Finally, a meron and an antimeron in close proximity can also form a stable bimeron.Figure 2(m) displays the corresponding simulated B z image of an isolated bimeron, while Fig. 2(n) shows the measured B z image of one such occurrence.As before, Fig. 2(o) displays the simulated B z image from its reconstructed m xy .It is interesting to note that while (anti)merons are always topologically protected, this is not necessarily true for meron-antimeron pairs.Labelling a meronantimeron pair as topologically protected would require the knowledge of the sign of the topological charge of its constituents [28,34].

EMERGENT MAGNETIC CHARGE
Reconstructing the m xy distribution from B z images is useful in understanding the underlying spin textures, but the fact that DQM actually provides a direct measure of ∇ • m xy creates a unique opportunity to consider a magnetic analog of the electric Gauss's law.Namely, the non-zero divergence of the magnetisation manifests the existence of an areal 2D magnetic charge density via σ m = −t( ∇ • m xy ).Here, t is the film thickness and m xy is independent of t -valid in the thin film limit [17].Therefore, AFM textures in α-Fe 2 O 3 have associated emergent magnetic charge distributions, which locally act as sources or sinks of the magnetic field.In fact, we can define a formal duality relation that connects the 2D magnetic charge density σ m to the staggered vorticity V via which scales with sin(ξ a ), highlighting the influence of the texture phase.Crucially, retrieving the full 2D emergent charge density σ m only requires the Fourier deconvolution of the measured B z images from the transfer function α xy (see S4.2).We can also perform a downward (upward) continuation [35] (see S4.2) of the planar B distribution captured in Fig. 2 (b, e, h, k, n), down to (away from) the sample surface.This allows a threedimensional visualisation of B (= H in vacuo) in the volume above the emergent magnetic charge distribution.Figure 3 panels (a) to (d) illustrate the representative field lines of B above a (anti-)clockwise a-Bloch meron, an ADW and an antimeron, respectively.For the two Bloch merons, B is consistent with the profile expected for spatially extended sources and sinks of magnetic field emanating from a monopolar distribution, which we retrieved via Eq. 2 in panels (e) and (f).This implies that a-Bloch merons host a class of emergent monopoles piggybacking on the whirling AFM textures.Interestingly, such monopolar magnetic distributions are not observed in ferromagnetic materials, as the presence of long-range demagnetizing fields favours divergence-free Bloch textures (e.g. in FM vortices).In our case, we are freed from this constraint due the presence of a weak demagnetizing contribution relative to the much stronger AFM exchange in α-Fe 2 O 3 .In contrast to merons, the ADW in Fig. 3(c) and the antimeron in Fig. 3(d) are associated with σ m distributions that exhibit dipolar and quadrupolar characters, respectively, as illustrated in panels (g) and (h).Finally, we emphasise that the observation of emergent monopoles is fully consistent with the simulations of AFM topological textures in Fig. 2 and does not violate Maxwell's equation, ∇ • B = 0.
On the basis of the above analysis, it is tempting to attribute to isolated a-Bloch merons a non-zero net monopolar charge quantified by Q m ≡ S σ m ds within an area S. We pick a circular integration area S of radius r centred on a given spin texture, as shown in Fig. 4(a and b) as an example in the case of antimerons.The particular (1/r) dependence of ∇ • m xy for two-dimensional magnetic charges hosted by spin textures then yields (see S5.9) Figure 4(c) presents the radial dependence of Q m for four measured merons (light blue dashed curves) and antimerons (light red dashed curves), while the dark blue (red) dashed curve is the average Q m radial dependence for merons (antimerons).For an isolated linear meron model, Q m scales linearly with r, and the measured Q m radial dependence is in agreement with this predicted behaviour.Q m itself is not a topological invariant, as a smooth transformation of an a-Bloch meron to an a-Néel meron would tune Q m from non-zero to zero value.For an isolated antimeron, the two-fold rotational symmetry should ensure Q m = 0 for all r values (see S5.9) and at short distances the experimentally obtained Q m value of an antimeron agrees well with this prediction.Beyond typical measurement-induced variations, deviation from the strict Q m = 0 condition arises when an antimeron is influenced by neighbouring spin textures.This reduces the two-fold symmetry and reveals finite Q m .Figures 4(a) and 4(b) capture this reduced symmetry and the consequential bias towards negative and positive Q m , respectively.Finally, as r goes to infinity, the integration area starts would overlap with surrounding spin textures, which leads to a further deviation from the Q m = 0 condition for an isolated antimeron model.
The assumption that we have a collection of isolated spin textures somewhat oversimplifies the reality.While the AFM topological textures are mesoscopically discernible via DQM and X-ray techniques and thus appear localised, they are in fact the constituents of the complex multi-textural ensemble that interacts via a 2D magnetic charge canvas.The magnetic charge per constituent is therefore not just dictated by their nature as merons and antimerons, but modified through their interaction with other constituents in the ensemble.For example, an isolated meron and antimeron pair forming an AFM bimeron, as shown in Fig. 2(n), would have a non-zero Q m , whose sign is determined predominantly by its meron (see S5.9).However, this clearly cannot be the case for a bimeron embedded in a uniform magnetisation field in the far field (see S5.10), since the divergence theorem ensures that Q m = 0.This goes to show that the interaction among AFM textures produces additional magnetic charge density away from the cores and highlights the interacting nature of this multi-textural ensemble.

DISCUSSION AND OUTLOOK
Our ability to identify the duality between topological AFM textures and magnetic charges is due to the direct readout of staggered vorticity enabled by DQM, as the NV centre senses both amplitude and orientation of the weak magnetic field.As such, we go beyond detecting antimerons and merons, to further distinguish between clockwise and anti-clockwise a-Bloch components, which otherwise appear indistinguishable in X-ray dichroic images.Our imaging approach can be easily extended to a wider family of topological textures, including a-Néel merons and bimerons, skyrmions, and distorted AFM textures that are otherwise divergence-free (see S5.6-5.8),relevant for topological AFM circuitry [7,26,36].
Although hematite provides especially favourable conditions for DQM imaging due to DMI-induced spin canting, it is by no means unique among AFMs in possessing a weak net magnetisation or quasi-isotropic spins in 2D [37][38][39][40].In fully compensated antiferromagnets without bulk DMI, staggered spin textures can also generate local net magnetisation, either statically or dynamically [41][42][43].Moreover, the DQM technique can be useful in detecting preferential vorticity in ultra-thin films induced by interfacial interactions --a key requirement for applications of quantum materials in topological spintronics [26,27].
In light of our discovery of the duality between magnetic charges and topological AFM textures, it is tempting to compare hematite with other systems that harbour emergent magnetic monopoles, such as the pyrochlore spin ice [25].Although intriguing, monopoles in spin ice are intrinsically distinct, as they have an underlying gauge charge, which is topological and quantised.On the contrary, the emergent magnetic charges in hametite are two-dimensional, not quantised and are only topological in the sense that they dress topological AFM textures underpinning them.In fact, we have demonstrated that hematite supports a rich tapestry of interacting magnetic charge distributions which could open up new and complementary ways to detect, manipulate and functionalise -via their magnetic charge -AFM topological textures.For example, our capability to classify different AFM spin textures could be combined with conventional spin manipulation techniques such as various spin torques, allowing for natural read-out and write-in schemes.Realizing this in a complex manifold of states endowed with highly non-linear interactions could be attractive for unconventional computing [1,2].Finally, the intriguing physical insights revealed in α-Fe 2 O 3 are a testament to the relevance and potential of DQM as a versatile table-top platform to explore emergent phenomena in antiferromagnets and other quantum materials.

C2. NV-to-Sample Distance
The NV-to-sample distance d N V is determined by measuring the magnetic field emitted across the edge of an out-of-plane magnetic material with the sensor of interest.Our calibration sample is a CoFeB magnetic strip (details in [44]) with a magnetization that remains saturated at remanence.This allows for accurate d N V value extraction following the proposal in [45].An example ODMR linescan across the edge of the magnetic strip, obtained with the diamond sensor used for imaging in Figs.1-3 in the main text, is given in Fig. S7 (a).A distribution of extracted d N V values from multiple ODMR linescans, such as the one in (a), is given in Figure S7 (b), resulting in an average d N V value of 70 ± 6 nm which sets the effective spatial resolution of our imaging.This is sufficient to study the AFM topological textures and piggybacked charges in the study, whose length scales are larger.
Appendix D: Magnetic Field Analysis

D1. Magnetic Field Components Above Sample Surface
Assuming a density of magnetic dipoles m(x, y) of thickness t, the magnetic field B(x, y) generated at a distance d above the sample surface (source-free region) can be described in Fourier space via the dipolar tensor [30,46,47]: In the absence of a canted m z component in α-Fe 2 O 3 , we have: In equation D3 we use the convolution theorem, and define ) and F and F −1 are the Fourier transform and its inverse.In real space, one can intuitively interpret equation D3 as the convolution of a point spread function, α xy [31], with the divergence of the in-plane canted magnetization, ∇ • m xy .In fact, as the distance d above the sample surface increases, spatial features of ∇ • m xy get increasingly 'blurred out', decreasing the effective spatial resolution of the NV sensor.A general discussion involving all components of the magnetization can be found in [31].

D2. Downward (Upward) Continuation and Charge Retrieval
Given the constrains set out in section D D1 and following equation D1, one can write the relation [35,48]: where B = F ( B).A positive or negative ∆d in equation D4 corresponds to an upward or downward continuation of the B-field distribution captured at d. Notice that when ∆d is positive, the exponential term e −k∆d acts to smooth out high frequency components and physically translates to losing spatial resolution.In the negative case, the exponential term tends to amplify high frequency noise present in the measurements.For this reason, we apply a Hanning window with a cut-off frequency at 2π/d N V , which is a common method to suppress the undesired amplification of high-frequency noise [23,48,49].By applying the downward continuation method down to the sample surface, it is therefore possible to retrieve the σ m = −t ∇ • m xy given by: where αxy = F (α xy ).This is a deconvolution process with the function α xy .
this, we will now describe the spatial distributions of the canted moment in various AFM quasiparticles considering their respective linear ansatz.Herein, we assume a canted moment m ∆ ∼ 2 × 10 3 A/m, and a canting angle ∆ ∼ 1.1 mrad.This is consistent with values reported for bulk α-Fe 2 O 3 (Table 4.2 of [50]) and comparable to the magnetic response of our samples, as discussed in our previous studies [17,33].In the small angle approximation, the average magnitude of the sublattice magnetization M s = m ∆ /(2∆) ∼ 9.1 × 10 6 A/m.

E1. Magnetic Field and Charge of AFM Antiphase Domain Walls
Assuming a linear ansatz for the antiphase domain wall profile centered at x = 0 along the x-axis, with phase ξ a and width w, we have: where, Hence,

E2. Magnetic Field and Charge of AFM In-plane Domain Walls
Uniform IP AFM domains do not generate any stray fields since their magnetisation would be non-divergent (see main text Eq. 1).However, in-plane domain walls above the Morin transition, which are expected in α-Fe 2 O 3 to span 60 • and 120 • , can generate stray fields.Their stray field signature depends on the domain wall width w, winding direction l 0 , internal phase of the domain wall -given by the angular span η 0 , and phase offset φ 0 : where, Ψ = Hence, IP domain walls generate different stray fields depending on their internal phase as illustrated in Fig. S9.

E3. Magnetic Field and Charge of AFM Merons and Antimerons
Assuming a linear ansatz for an (anti)meron centered at the origin, with phase ξ a , topological winding number N , and radial size R M , we have: where, θ = Hence, The last term is referred to as a radial function f (r) in the main text.

E4. Magnetic Field Simulations of Various AFM Entities
We have simulated the magnetic fields generated d = 70 nm above m grids containing antiphase domain walls in Fig. S8, 60 • and 120 • in-plane domain walls in Fig. S9, and merons and antimerons in Fig. S10.

E5. Distinguishing Topological and non-topological AFM textures
Above the Morin transition, topological textures are always intertwined with IP domain walls -in fact, each topological texture is a 'pinch-point' of IP domain walls [17,18].It is therefore crucial to be able distinguish them apart.Whilst IP domain walls can generate stray fields (fig.S9), their strength, distribution, and evolution are markedly different from counterparts generated by merons and antimerons (Fig. S10).In the following we outline the salient differences and how one can distinguish them: • Locally, topological merons and antimerons are two-dimensional entities which produce confined stray fields, characteristic of monopolar and quadrupolar distributions.Alternatively, IP domain walls are deconfined quasione-dimensional trivial texturesThis is illustrated in Fig. S11 using the linear meron ansatz.
• The maximum stray field magnitude generated by merons and domain walls is (i) inversely proportional to their characteristic length scale (e.g., meron core radius; IP domain wall width, W ) and (ii) directly proportional to  the total angle they span in spin space ( π for meron; η 0 for IP domain wall angle) and is given by: Since merons trap out-of-plane cores, the anisotropy energy penalty causes these cores to be small [17].On the other hand, IP domain walls are much wider as they do not have out-of-plane regions and are only limited by the very weak basal anisotropy in α-Fe 2 O 3 [51].This causes R M < W , at temperatures well above the Morin transition, which is the case at 300 K.Moreover, since IP domain walls are expected to have η 0 ≤ π.Consequently, the magnitude of the stray field generated by a-Bloch merons at their core is significantly more intense, by a factor of ∼ 2 compared to trivial textures.This allows unambiguous identification of topological meron cores.
• Although both topological and trivial textures are present above the Morin transition, application of in-plane magnetic fields results in the large-scale annihilation of merons and antimerons only.This triggers the selective elimination of strong and confined monopolar stray fields, leaving behind trivial textures, as illustrated in SI section H and our previous work [17,18].This experimentally evinces the different evolution of topological and non-topological entities.
• As for antimerons, the quadrupolar stray field distribution with a zero-field cross pattern (Fig. S10) while distinct from IP domain walls (Fig. S9), can be challenging to disambiguate in isolation.For this reason, we restrict antimeron identification and fitting to regions in the immediate proximity of strong monopolar field distributions generated by merons.This is justified, given that merons and antimerons occur with equal probability after executing a Kibble-Zurek transition [17,18] and that, based on our X-ray work, meron-antimeron pairs are extremely common.

E6. Invisibility of Divergence-free Textures
For a divergence-free texture (i.e.ik x mx + ik y my = 0), the z-component of the field in Eq.D2 vanishes, rendering it invisible to magnetic field detection.This is further illustrated in Figs.S12 and S13, where ADWs and merons turn invisible when ξ a = 0, π.In reality, textures are never perfectly isolated, and the distortion produced by the proximity of other textures will make even a-Néel meron visible through the resulting magnetic field.This is illustrated in figure S14, which display the magnetic field of an a-Néel meron undergoing progressive geometrical distortion, and will be further discussed in Fig. S19 in the context of a-Néel bimerons.

E7. Magnetic Field and Charge of AFM Bimerons
An AFM bimeron in a uniform m xy background can be described analytically as a skyrmion model rotated by 90 • along an in-plane axis.Assuming a phase ξ a and radial size R B , we have:  where, θ Hence, Equation E14 allows us to model the σ m distribution of bimerons, see Fig. S15.Here, we observe that individual monopolar and quadrupolar σ m distributions contributed by the meron and antimeron components, albeit significantly distorted in comparison to their isolated counterparts, can be discerned in the σ m distribution of an isolated AFM bimeron.This is slightly more evident in the case of a-Bloch bimerons (ξ a = π/2, 3π/2).For a-Néel bimerons   the σ m distribution only reflects a strongly distorted quadrupolar character as the meron component is relatively divergence-free.More importantly, we observe a non-zero σ m distribution, irrespective of the ξ a value, implying that bimerons of all characters generate distinct and detectable B z signatures, see Fig. S16.This enables DQM to play a versatile role in the future discovery and study of isolated AFM bimerons in engineered AFM systems with a defined in-plane anisotropy and inhomogeneous DMI.Interestingly, while the integrated charge Q m over a sufficiently large area centred at the core of the bimeron vanishes for all ξ a , the asymmetric σ m distribution suggests that it may be possible to spatially manipulate AFM bimerons via external magnetic fields.

E8. Magnetic Field and Charge of AFM Skyrmions
We can adapt the treatment from the previous section to construct the magnetic distribution in AFM skyrmions.Assuming a phase ξ a and radial size R S , we have: where, θ Hence, In actuality, Néel skyrmions undergo deformation/distortion [52], which could result in detectable magnetic field signatures, akin to what is expected for distorted a-Néel merons, as illustrated in Fig. S14.

E9. Integrated Magnetic Charge of AFM Merons and Antimerons
The integrated charge Q m of an isolated (anti)meron with phase ξ a and winding number N as a function of radial distance r from the core can be analytically written as: For isolated symmetric textures, Q m is either zero for antimerons (N = −1) and a-Néel merons (N = +1, ξ a = 0, π), or accumulates radially (N = +1, ξ a = 0, π), as illustrated in Fig. S19(d).
Moreover, the integrated charge Q m of an isolated bimeron with phase ξ a can be written as: While isolated bimerons with ξ a = 0, π accumulate some charge at short distances near their cores, see Fig. S19(d), the total integrated charge, outside R B , for all isolated bimerons embedded in a uniformly aligned magnetic background is identically is zero.are linearly dependent on the first one.We chose to work with B z since it intuitively reveals a low-pass filtered ∇ • m xy (eqn.D3).The inverse of the map f such that B z = f ( m), is not still not unique, even within the restricted class of problems represented by the dipolar tensor (eqn.D1).However, a choice of appropriate constraints reduces the degrees of freedom, thereby allowing us to estimate m from a DQM scan.In the present case, we set the out-ofplane magnetization component m z = 0.The inversion problem is thus reduced from three to two dimensions, but the reduced dimensionality is still not sufficient to estimate the in-plane magnetization vector.This is clear when considering equation D2, where we show that any divergence-free m xy , including trivial ordering, produces the same field distribution (|B| = 0, section E E6).

G2. Regularization
We adopt several regularisation strategies developed by the Lorentz transmission electron microscopy community to approximate a solution.Closely following the strategy in [53], we replace the ill-posed inversion problem with the following minimization problem: where, F, B z and m are the matrix representation of the map f , the observed B z and the reconstruction of m, respectively.The spatial discretisation, implicit to the minimisation equation and also to a DQM scan, reduces the problem to finite dimensions and limits the maximum resolvable signal frequency to the inverse of the spatial resolution (i.e.< 1  33 GHz in our images).The first term in equation G1 is the squared error needed for the base minimisation.The second term, known as the Tikhonov regularisation, is used to impose a smoothness constraint to the solution m.This is achieved via Γm where Γ = λ∇ is the matrix representation of the gradient ∇ weighted by the constant λ.   to one dimension alone does not guarantee a unique solution as we have a non-trivial null space.Finally, to quicken the minimisation and to avoid local minima, we initialize the problem with a m guess based on a priori information of our system [17].φ m is initialized based on the expected topological texture associated with the characteristic B z distribution shown analytically and via simulations in Fig. 2 in the main text, whereas | m xy | is fixed based on the linear meron ansatz (section E) with an average core size of 150 nm (section F), and the linear ADW ansatz with a similar domain wall width (section E2).In this respect, we only attempt to reconstruct m xy for AFM textures and their immediate neighbourhood with characteristic B z signatures.Additionally, we minimize the problem with B z at

G4. Fitting of (Anti)Merons and Bimerons
Similar to the procedure in section G G3, | m xy | = 1 is assumed everywhere except in the (anti)meron core, which is modelled with a linear ansatz with R M = 150 nm (section F).We show the m xy reconstruction of a simulated a-Bloch AFM meron and an AFM antimeron in Figs.S24 and S25, respectively.In both cases, we start off with ξ a that deviates from the simulated texture, and the corresponding m xy guess is given in panel (c).The resultant minimized solution given in panel (d) closely reproduces the simulated m xy distribution in panel (a).The respective residues are given in panels (e-f).
Finally, we apply the same reconstruction protocol to experimental images, from which we obtain the m xy fits of an anti-clockwise a-Bloch meron (Fig. S26), an antimeron (Figure S27) and a composite structure consisting of one meron and two antimerons (Fig. S28).

G5. Field imaging and Magnetisation Reconstruction of α-Fe2O3
Unique magnetisation reconstruction is a general problem in the stray-field imaging.Nonetheless, the rigorous identification procedure (see E E5), together with the implementation of systematic regularization protocols, has helped constrain the magnetic solution space and provided meaningful canted moment reconstructions.While these reconstructions may not always be unique, they are in practice entirely consistent with AFM textures reported via  X-ray imaging [17].
Although X-ray linearly dichroic imaging is not affected by the same uniqueness problem, it suffers from an intrinsic 180 • ambiguity which would make reconstructions of the canted moment more ambiguous than in the present case.Additionally, the canted magnetisation in α-Fe 2 O 3 is too weak to produce circular dichroic X-ray contrast [17].Therefore, field imaging via DQM provides an invaluable view of the rich magnetic charge phenomenology in α-Fe 2 O 3 .Upon application of strong in-plane magnetic fields, it is possible to drive AFM domain realignment which results in meron-antimeron annihilation [17].To confirm this, two permanent magnet poles were aligned parallel to the sample surface resulting in a homogeneous field across the sample, variable via micrometer stages.Since imaging under high off-NV-axis field is not possible, the external field was applied for a few minutes followed by NV imaging  in zero field, see Fig. S29.Although the scans are approximately in the same area, open-loop scanning requires post processing, wherein the images have been registered using distinct field features in the scans.Several merons and antimerons are initially present.As the magnitude of the external field increases, annihilation processes occur, accompanied by an alteration of the surrounding trivial AFM textures.Due to the limited field magnitude that can be applied in-situ, large-scale annihilation of textures was driven by an ex-situ 400 mT field [17], see Fig. S30.Due to the probabilistic nature of the generation of topological textures via the Kibble-Zurek mechanism, the occurrence of merons with different projected vorticities (red vs blue cores) can vary across the sample [18].It is noted that while Fig. S30 shows predominately anti-clockwise a-Bloch merons (red cores), additional large area scans (see S9) confirm approximately similar distributions of merons with different projected vorticities (red vs blue cores) are obtained due to the absence of chiral interactions in α-Fe 2 O 3 [17].

Figure 1 .
Figure 1.Signatures of emergent magnetic field in hematite α-Fe2O3.(a) Atomic structure of α-Fe2O3 (Fe and O atoms in yellow/green and grey spheres respectively).(b) Discrete representation of the alternating FM sublattice magnetisation M1 (yellow cones) and M2 (green cones) with AFM coupling along the c-axis shown in (a).(c) Illustration of the whirling staggered magnetisation, l (grey cones), forming an anti-clockwise a-Bloch meron, and the resultant canted magnetic moment, m (red cones).The insert shows the relationship between l, m, M1, M2 and the canting angle ∆.(d) A scanning diamond sensor with a single NV centre maps out the B-field generated near the sample surface.The insert shows the energy diagram of the NV ground state (GS) |±1 and |0 sublevels.A microwave field drives the GS spin transition while a 532-nm laser excites the NV to the excited state (ES) (green arrow).The NV then undergoes a radiative decay to GS (red arrow) or a non-radiative and spin selective path via the intersystem crossing (ISC) (blue arrow) enabling ODMR acquisition.(e, f ) ODMR along the fast scan direction, measured on α-Fe2O3 thin film at T = 4 K (e) and across TM at 300 K (f).The fitted f+ (BNV) is plotted as a white line in each panel.(g, h) Bz images retrieved from fitted BNV maps reveal distinct field signatures across TM.Dashed lines in (g) and (h) correspond to fitted BNV traces in (e) and (f), respectively.Scale bars are 1 µm.

Figure 2 .
Figure 2. Classification of topological AFM textures via DQM.(a-c) Topological AFM textures observed below TM.Distinct Bz signature of an ADW simulated (a) and measured (b) above the sample surface.The reconstructed mxy (black arrows) from (b) and its Bz distribution are given in (c).(d-o) Topological AFM textures observed above TM.Simulated (d) and measured (e) Bz field signatures of an anti-clockwise a-Bloch meron and (f) the reconstructed mxy (black arrows) of (e) and its associated Bz.Similarly, the simulated and measured Bz signatures, and the mxy reconstruction of a clockwise a-Bloch meron, an anti-meron and a bimeron is given in panels (g-i), (j-l), (m-o) respectively.Details of mxy reconstruction are given in Ref. (see S7). Scale bars are 200 nm.

Figure 3 .
Figure 3. Emergent magnetic charge distributions.(a-d) Three-dimensional visualisation of B in the volume above an assortment of topological AFM textures.Streamtubes illustrate the magnetic field lines of B above (a) an anti-clockwise a-Bloch meron, (b) a clockwise a-Bloch meron, (c) an ADW and (d) an antimeron.The girth and colour of each streamtube vary with the magnetic field norm |B| and the z-component of the field Bz, respectively.(e-h) Magnetic charge density σm distributions retrieved from downward continuation of (a-d) reveal a magnetic monopolar, anti-monopolar, dipolar and quadrupolar charge character associated to an anti-clockwise a-Bloch meron, a clockwise a-Bloch meron, an ADW and an antimeron respectively.Scale bars are 200 nm.

Figure 4 .
Figure 4. Scaling of two-dimensional integrated magnetic charges.(a-b) Reconstructed magnetic charge distribution, σm, of two experimentally observed antimerons (AM 1 and AM 2) with slightly distorted quadrupolar characters.AM 1 and AM 2 display a respective bias towards positive and negative charges, respectively.Dashed circle in (a) and (b) illustrates the circular integration area S of radius r, centred at the core of the antimeron, to obtain |Qm|.Scale bars are 200 nm.(c) Experimentally retrieved magnitude of the total integrated magnetic charge, |Qm| of multiple merons (|Qm|M, light blue dashed curves) and antimerons (|Qm|AM, light red dashed curves) plotted as a function of integration radius, r.Their average experimental |Qm| profiles, |Qm| avg M and |Qm| avg AM , are represented by dark blue and dark red curves, respectively.Solid black curves plot the theoretically predicted Qm radial dependence based on Eq. 3 for isolated merons and antimerons.

Figure S6 .Figure S7 .
Figure S6.NV orientation characterisation.(a) ODMR with 0.5 mT applied along the NV axis showing the Zeeman-split resonances with additional hyperfine splitting due to the host 15 N. ODMR spectrum as a function of (b) φB with θB fixed at 180 • and (c) θB with φB = 96 • .The maximum splitting in (b) determines φNV = 96 • which is then used in (c) to identify θNV = 120 • .

Figure S8 .
Figure S8.Magnetic field simulation of one-dimensional topological textures.Simulated Bz exhibited by ADWs upon varying ξa.Black arrows illustrate mxy.
Figure S9.Magnetic field simulation of trivial in-plane domain walls.Simulated Bz exhibited by 60 • and 120 • in-plane domain walls with various phase offset.Black arrows illustrate mxy.

Figure S11 .
Figure S11.Linear charge profile merons and IP domain walls.The linear charge profile λM is the integrated charge along the y-axis given as a function of the x-axis of the meron in Fig. S10 and the IP domain walls (IPDW) in Fig. S9.The λM of the meron and IPDWs reveal their respective confined and deconfined nature.

Figure S12 .
Figure S12.Field magnitude of ADWs with different phase ξa.Field magnitude | B| obtained from the magnetic field simulations of ADWs with different ξa.Black arrows illustrate mxy.

Figure S13 .
Figure S13.Magnetic field magnitude of merons with different phase ξa.Field magnitude | B| obtained from the magnetic field simulations of merons with different ξa.Black arrows illustrate mxy.

Figure S14 .
Figure S14.Distortion of a divergence-free texture.Field magnitude | B| exhibited by an a-Néel meron with increasing distortion to the surrounding mxydistribution from (a) to (d).Black arrows illustrate mxy.

Figure S15 .
Figure S15.Magnetic charge density σm of an AFM bimeron with varying phase ξa.Black arrows illustrate mxy.

Figure S16 .
Figure S16.Magnetic field distribution of an AFM bimeron with varying phase ξa.Black arrows illustrate mxy.

Figure S17 .
Figure S17.Magnetic charge density σm distribution of an AFM skyrmion with varying phase ξa.Black arrows illustrate mxy.
) AFM Néel skyrmions (ξ a = 0, π) have vanishing σ m , as they host divergence-free m xy .Other skyrmions with ξ a = 0, π have a non-zero σ m distribution, see Fig. S17.As expected, B z mirrors the σ m trend, Fig. S18.While AFM Néel skyrmions are invisible, DQM can still play a crucial role in imaging AFM skyrmions with mixed or a-Bloch characters, with the added advantage of reading out the projected vorticity via the strength of B z .

Figure S19 .
Figure S19.Charge accumulation of isolated textures.Numerical simulation of charge density σm distribution of (a) an a-Bloch meron, (b) an antimeron and (c) a bimeron.Arrows represent mxy.(d) The total charge |Qm| accumulated in a circular area S of increasing radius r, for each numerical simulation in (a-c), is given by the solid dots.Qm from analytical equations E18 and E19 are represented by solid curves.The analytical Qm for a a-Bloch meron and an antimeron ansatz agree well with the numerically obtained Qm.
We constrain the problem further by fixing the magnitude of the canted magnetization | m xy |.Given that m xy can be written as | m xy |(cos φ m , sin φ m ), one only needs to reconstruct the angle φ m .Note that the reduction of the problem

Figure S20 .
Figure S20.Charge variability of a bimeron.(a-e) Numerical simulations of charge density σm distribution of bimerons with various mxy boundary conditions.Arrows represent mxy.(f) Total charge |Qm| accumulated in a fixed area S of various bimeron simulation with different mxy boundary conditions in (a-e).For convenience, S is defined as the square area covering the entirety of each panel, and is centered at the midpoint between the core of the meron and the antimeron.

Figure S21 .
Figure S21.Average meron core size.The average experimental Bz profile of six merons as a function of radial distance r from the center of the meron core (solid black dots, shaded area is the standard deviation).The experimental profile is compared with the simulated Bz profiles of merons with various core sizes RM (coloured lines).

Figure S22 .
Figure S22.Simulated ADW fit.(a) Simulated mxy of an AFM ADW and the resultant Bz at dNV = 70 nm.(b-c) (b) Fixed | mxy| mask and (c) guess mxy distribution used to initialize the minimisation problem.(d-f ) (d) Minimized mxy solution and respective (e) Bz and (f) mxy residues.Black arrows illustrate mxy.

Figure S24 .
Figure S24.Simulated a-Bloch meron fit.Simulated mxy of an AFM a-Bloch meron and the resultant Bz at dNV = 70 nm.(b-c) (b) Fixed | mxy| mask and (c) guess mxy distribution used to initialize the minimisation problem.(d-f ) (d) Minimized mxy solution and respective (e) Bz and (f) mxy residues.Black arrows illustrate mxy.

Figure S25 .
Figure S25.Simulated antimeron fit.Simulated mxy of an AFM antimeron and the resultant Bz at dNV = 70 nm.(b-c) (b) Fixed | mxy| mask and (c) guess mxy distribution used to initialize the minimisation problem.(d-f ) (d) Minimized mxy solution and respective (e) Bz and (f) mxy residues.Black arrows illustrate mxy.

Figure S26 .
Figure S26.Experimental a-Bloch meron fit.Measured Bz of an AFM right a-Bloch meron at dNV = 70 nm.(b-c) (b) Fixed | mxy| mask and (c) guess mxy distribution used to initialize the minimisation problem.(d-e) (d) Minimized mxy solution and respective (e) Bz residue.Black arrows illustrate mxy.

Figure S28 .
Figure S28.Experimental fit for composite texture.Measured Bz of a composite texture consisting of two antimerons and an a-Bloch meron at dNV = 70 nm.(b-c) (b) Fixed | mxy| mask and (c) guess mxy distribution used to initialize the minimisation problem.(d-e) (d) Minimized mxy solution and respective (e) Bz residue.Black arrows illustrate mxy.

Figure S29 .
Figure S29.Reversal field imaging at remanence.Remanent state of α-Fe2O3 after applying an in-situ in-plane magnetic field with magnitude varying from (a) 150mT, (b) 200mT and (c) 250mT.

Figure S30 .
Figure S30.Annihilation of topological spin entities.Remanent state of α-Fe2O3 after applying an ex-situ in-plane magnetic field of ∼ 400 mT.