Two-Dimensional Ferronematics, Canonical Harmonic Maps and Minimal Connections

We study a variational model for ferronematics in two-dimensional domains, in the “super-dilute” regime. The free energy functional consists of a reduced Landau-de Gennes energy for the nematic order parameter, a Ginzburg–Landau type energy for the spontaneous magnetisation, and a coupling term that favours the co-alignment of the nematic director and the magnetisation. In a suitable asymptotic regime, we prove that the nematic order parameter converges to a canonical harmonic map with non-orientable point defects, while the magnetisation converges to a singular vector field, with line defects that connect the non-orientable point defects in pairs, along a minimal connection.


Introduction
Nematic liquid crystals (NLCs) are classical examples of mesophases or liquid crystalline phases that combine fluidity with the directionality of solids [24]. The nematic molecules are typically asymmetric in shape e.g. rod-shaped, wedge-shaped etc., and these molecules tend to align along certain locally preferred directions in space, referred to as directors. Consequently, NLCs have a direction-dependent response to external stimuli such as electric fields, magnetic fields, temperature and incident light. Notably, the directionality or anisotropy of NLC physical and mechanical responses make them the working material of choice for a range of electro-optic applications [36]. However, the magnetic susceptibility of NLCs is much weaker than their dielectric anisotropy, typically by several orders of magnitude [19]. Hence, NLCs exhibit a much stronger response to applied electric fields than their magnetic counterparts and as a result, NLC devices are mainly driven by electric fields. This naturally raises a question as to whether we can enhance the magneto-nematic coupling and induce a spontaneous magnetisation by the introduction of magnetic nanoparticles (nanoparticles with magnetic moments) in nematic media, even without external magnetic fields. If implemented successfully, these magneto-nematic systems would have a much stronger response to applied magnetic fields, compared to conventional nematic systems, rendering the possibility of magnetic-field driven NLC systems in the physical sciences and engineering.
This idea was first introduced in 1970 by Brochard and de Gennes in their pioneering work on ferronematics [19] and these composite systems of magnetic nanoparticle (MNP)-dispersed nematic media are referred to as ferronematics in the literature [20,21,19]. The system has two order parameters -the Landau-de Gennes (LdG) Q-tensor order parameter to describe the nematic orientational anisotropy and the spontaneous magnetisation, M, induced by the suspended MNPs. Brochard and de Gennes suggested that the nematic directors, denoted by n, can be controlled by the surface-induced mechanical coupling between NLCs and MNPs. Equally, the spontaneous magnetisation, M profiles can be tailored by the nematic anisotropy through the MNP-NLC interactions, and this two-way coupling can stabilise exotic morphologies and defect patterns.
We work with dilute ferronematic suspensions relevant for a uniform suspension of MNPs in a nematic medium, such that the distance between pairs of MNPs is much larger than the individual MNP sizes and the volume fraction of the MNPs is small, building on the models introduced in [20,21] and then in [15,14]. In these dilute systems, the MNP-MNP interactions and the MNP-NLC interactions are absorbed by an empirical magneto-nematic coupling energy. These coupling energies can also be rigorously derived from homogenisation principles, as elucidated in the recent work [22]. We work with two-dimensional, simply-connected and smooth domains Ω, in a reduced LdG framework for which the Q-tensor order parameter is a symmetric, traceless 2 × 2 matrix and M is a two-dimensional vector field. This reduced approach can be rigorously justified using Γ-Convergence techniques (see [31] since in three dimensions, the LdG Q-tensor order parameter is a symmetric, traceless 3 × 3 matrix with five degrees of freedom). We use the effective re-scaled free energy for ferronematics, inspired by the experiments and results in [41] and proposed in [15,14]. This energy has three components -a reduced LdG free energy for NLCs, a Ginzburg-Landau free energy for the magnetization and a homogenised magnetonematic coupling term: In two dimensions, we have f (Q, M) := 1 4 (|Q| 2 − 1) 2 + ξ 4 (|M| 2 − 1) 2 − c 0 QM · M.
We work with a dimensionless model where ε 2 is interpreted as a material-dependent, geometrydependent and temperature-dependent positive elastic constant, ξ is a ratio of the relative strength of the magnetic and NLC energies and c 0 is a coupling parameter. ξ is necessarily positive, positive c 0 coerces co-alignment of n and M whereas c 0 < 0 coerces n to be perpendicular to M [14]. We only consider positive c 0 in this paper.
For dilute suspensions, ε and ξ are necessarily small. In [14], the authors study stable critical points of this effective ferronematic free energy on square domains, with Dirichlet boundary conditions for both Q and M. Their work is entirely numerical but does exhibit a plethora of exotic morphologies for different choices of ε, ξ and c 0 . They demonstrate stable nematic point defects accompanied by both line defects and point defects in M, and there is considerable freedom to manipulate the locations, multiplicity and dimensionality of defect profiles by simply tuning the values of ξ and c 0 . In particular, the numerical results clearly show that line defects or jump sets are observed in stable M-profiles for small ξ and c 0 , whereas orientable point defects are stabilised in M for relatively large ξ and c 0 . Motivated by these numerical results, we study a special limit of the effective free energy in (1.1), for which both ξ and c 0 are proportional to ε and we study the profile of the corresponding energy minimizers in the ε → 0 limit, subject to Dirichlet boundary conditions for Q and M. This can be interpreted as a "super-dilute" limit of the ferronematic free energy for which the magnetic energy is substantially weaker than the NLC energy, and the magneto-nematic coupling is weak. In the "super-dilute" limit, "ε" is the only model parameter and ξ, c 0 are defined by the constants of proportionality which are fixed, and hence ε → 0 is the relevant asymptotic limit. Our main result shows that in this distinguished limit, the minimizing Q-profiles are essentially canonical harmonic maps with a set of non-orientable nematic point defects, dictated by the topological degree of the Dirichlet boundary datum. This is consistent with previous powerful work in [9] in the context of the LdG theory and unsurprising since the LdG energy is the dominant energy. The minimizing M-profiles are governed by a Modica-Mortola type of problem, quite specific to this super-dilute limit [29]. They exhibit short line defects connecting pairs of the non-orientable nematic defects, consistent with the numerical results in [14]. These line defects or jump sets in M are minimal connections between the nematic defects, and the location of the defects is determined by a modified renormalisation energy, which is the sum of a Ginzburg-Landau type renormalisation energy and a minimal connection energy. The modified renormalisation energy delicately captures the coupled nature of our problem, which makes it distinct and technically more complex than the usual LdG counterpart.
We complement our theoretical results with some numerical results for stable critical points of the ferronematic free energy, on square domains with topologically non-trivial Dirichlet boundary conditions for Q and M. The converged numerical solutions are locally stable, and we expect multiple stable critical points for given choices of ε, ξ and c 0 .The numerical results are sensitive to the choices of ε and c 0 , but there is evidence that the numerically computed stable solutions do indeed converge to a canonical harmonic Q-map and a M-profile closely tailored by the corresponding Q-profile. The Q-profile has a discrete set of non-orientable nematic defects and the M-profile exhibits line defects connecting these nematic defects, in the ε → 0 limit. Whilst the practical relevance of such studies remains uncertain, it is clear that strong theoretical underpinnings are much needed for systematic scientific progress in this field, and our work is a first powerful step in an exhaustive study of ferronematic solution landscapes [47] (also see recent work in [23], [40]).
The paper is organised as follows. In Section 2, we set up our problem and state our main result, recalling the key notions of a canonical harmonic map and a minimal connection. In Section 3, we state and prove some key technical preliminary results. In Section 4, we prove the six parts of our main theorem, including convergence results for the energy-minimizing Q and M-profiles in different function spaces, and the convergence of the jump set of the energyminimizing M to a minimal connection between pairs of non-orientable nematic defects, in the ε → 0 limit. The defect locations are captured in terms of minimizers of a modified renormalized energy, which is the sum of the Ginzburg-Landau renormalized energy and a minimal connection energy. The modified renormalized energy is derived from sharp lower and upper bounds for the energy minimizers in the ε → 0 limit, in Sections 4.4.1 and 4.4.2. In Section 5, we present some numerical results and conclude with some perspectives in Section 6.

Statement of the main result
Let S 2×2 0 be the set of 2 × 2, real, symmetric, trace-free matrices, equipped with the scalar product Q · P := tr(QP) = Q ij P ij and the induced norm |Q| 2 := tr(Q 2 ) = Q ij Q ij . Let Ω ⊆ R 2 be a bounded, Lipschitz, simply connected domain. The "super-dilute" limit of the ferronematic free energy is defined by ξ = ε; c 0 = βε where β, ε are positive parameters. For Q : Ω → S 2×2 0 and M : Ω → R 2 , we define the functional where the potential f ε is given by and κ ε ∈ R is a constant, uniquely determined by imposing that inf f ε = 0. We consider minimisers of (2.1) subject to the Dirichlet boundary condition at any point of ∂Ω. Here I is the 2 × 2 identity matrix. The assumption (2.4) implies that the potential f ε , evaluated on the boundary datum (Q bd , M bd ), takes nonzero but small valuesthat is, we have f ε (Q bd , M bd ) > 0 for ε > 0 but f ε (Q bd , M bd ) → 0 as ε → 0. (For details of this computation, see Lemma B.3 in Appendix B.) Throughout this paper, we will denote by (Q * ε , M * ε ) a minimiser of the functional (2.1) subject to the boundary conditions (2.3). By routine arguments, minimisers exist and they satisfy the Euler-Lagrange system of equations We denote as N the unit circle in the space of Q-tensors, that is, Equivalently, N may be described as is a real vector space of dimension 2, the set N is a smooth manifold, diffeomorphic to the unit circle S 1 ⊆ C. A diffeomorphism is given explicitely by By assumption, the domain Ω ⊆ R 2 is bounded and convex, so its boundary ∂Ω is parametrised by a simple, closed, Lipschitz curve -in particular, ∂Ω is homeomorphic to the circle S 1 . Therefore, the boundary data (Q bd , M bd ) carries a well-defined topological degree In principle, for a continuous map Q : ∂Ω → N , the degree may be a half-integer, that is deg(Q, ∂Ω) ∈ 1 2 Z. However, the boundary datum Q bd is orientable, by assumption (2.4) -in fact, it is oriented by M bd . This explains why d, in our case, is an integer.
Remark 2.1. The results in this paper -in particular, our main result, Theorem 2.1 belowremain true for slightly different choices of the boundary conditions. For instance, we could consider minimisers of the functional (2.1) in the class of maps Q ∈ W 1,2 (Ω, S 2×2 0 ) that satisfy Q = Q bd on ∂Ω, where the boundary datum Q bd takes the form and deg(n bd , ∂Ω) = d, but we do not impose any relation between n bd and the value of M at the boundary. In this case, minimisers of the functional will satisfy the natural (Neumann) boundary condition ∂ ν M ε = 0 on ∂Ω for the M-component, where ∂ ν is the outer normal derivative. The arguments carry over to this case, with no essential change (see also Remark 4.2).
The canonical harmonic map and the renormalised energy. In order to state our main result, we recall some terminology introduced by Bethuel, Brezis and Hélein [11]. Although the results in [11] are stated in terms of complex-valued maps, as opposed to Q-tensors, they do extend to our setting, due to the change of variable (2.9). Let a 1 , . . . , a 2|d| be distinct points in Ω (with d given by (2.10)). We say that a map Q * : Ω → N is a canonical harmonic map with singularities at (a 1 , . . . , a 2|d| ) and boundary datum Q bd if the following conditions hold: . . , a 2|d| } and Q * = Q bd on ∂Ω; (ii) for any σ > 0 small enough and any j ∈ {1, . . . , 2 |d|}, we have in the sense of distributions in the whole of Ω. (Here and in what follows, we adopt Einstein's notation for the sum).
If B ⊆ Ω \ {a 1 , . . . , a 2|d| } is a ball that does not contain any singular point of Q * , then Q * can written in the form where θ * : B → R is a smooth function. (Equation (2.12) follows from (2.7), by classical lifting results in topology.) Then, the equation (iii) above can be written in the form In other words, a canonical harmonic map can be written locally, away from its singularities, in terms of a harmonic function. The canonical harmonic map with singularities at (a 1 , . . . , a 2|d| ) and boundary datum Q bd exists and is unique, see [11, Theorem I.5, Remark I.1]. The canonical harmonic map satisfies Q * ∈ W 1,p (Ω, N ) for any p ∈ [1, 2), but Q * / ∈ W 1,2 (Ω, N ). Nevertheless, the limit exists and is finite (see [11,Theorem I.8]). Following the terminology in [11], the function W is called the renormalised energy.
Minimal connections between singular points. Given distinct points a 1 , a 2 , . . . , a 2|d| in R 2 , we define a connection for {a 1 , . . . , a 2|d| } as a finite collection of straight line segments {L 1 , . . . , L |d| } such that each a j is an endpoint of exactly one of the segments L k . In other words, the line segments L j connects the points a i in pairs. We define (2.15) Here and throughout the paper, H 1 denotes the 1-dimensional Hausdorff measure (i.e., length). We say that a connection {L 1 , . . . , L d } is minimal if it is a minimiser for the right-hand side of (2.15). A notion of minimal connection, similar to (2.15), was already introduced in [17,2]. However, the minimal connection was defined in [17] by taking the orientation into account -that is, half of the points a 1 , . . . , a 2|d| were given positive multiplicity 1, the other half were given negative multiplicity −1, and the segments L j were required to match points with opposite multiplicity. By constrast, here we do not distinguish between positive and negative multiplicity for the points a i and any segment of endpoints a i , a k is allowed. (In the language of Geometric Measure Theory, the minimal connection was defined in [17] as the solution of a 1-dimensional Plateau problem with integer multiplicity, while (2.15) is a 1-dimensional Plateau problem modulo 2.) The main result. We prove a convergence result for minimisers (Q * ε , M * ε ) of (2.1), subject to the boundary conditions (2.3)-(2.4), in the limit as ε → 0. We denote by SBV(Ω, R 2 ) the space of maps M = (M 1 , M 2 ) : Ω → R 2 whose components M 1 , M 2 are special functions of bounded variation, as defined by De Giorgi and Ambrosio [25]. The distributional derivative DM of a map M ∈ SBV(Ω, R 2 ) can be decomposed as where ∇M : Ω → R 2×2 is the absolutely continuous component of DM, L 2 (dx) is the Lebesgue measure on R 2 , S M is the jump set of M, M + , M − are the traces of M on either side of the jump and ν M is the unit normal to the jump set. (See e.g. [3] for more details). (i) Q * ε → Q * strongly in W 1,p (Ω) for any p < 2; (ii) Q * is the canonical harmonic map with singularities (a * 1 , . . . , a * 2|d| ) and boundary datum Q bd ; (iii) M * ε → M * strongly in L p (Ω) for any p < +∞; (iv) M * ∈ SBV(Ω, R 2 ) and it satisfies among all the (2 |d|)-uples (a 1 , . . . , a 2|d| ) of distinct points in Ω.
Remark 2.2. Theorem 2.1 implies that M * is a locally harmonic map, away from the closure of its jump set, into the circle of radius ( √ 2β + 1) 1/2 . In other words, if B is a ball that does not intersect the closure of the jump set of M * , then M * can locally be written in the form M * = ( √ 2β + 1) 1/2 (cos ϕ * , sin ϕ * ) for some scalar function ϕ * : B → R that satisfies −∆ϕ * = 0 in B. See Proposition 4.12 for the details. Remark 2.3. Let us discuss the extremal case of the renormalized energy W β (a 1 , . . . , a 2|d| ). When β → +∞, the function W β would be minimized by choosing (L * 1 , . . . , L * |d| ) to be zero, meaning that the singular points will move toward each other. In the case where β = 0 instead, the coupling term in the potential wouldn't be present. Therefore, we would have two decoupled Ginzburg-Landau problems.
Remark 2.4. Point defects and line defects connecting point defects do appear for energy minimizers in other variational models e.g. continuum models for a complex-valued map in [30] or for discrete models in [6,5]. However, the mathematics is substantially different to our model problem for which we have two order parameters Q and M, and a non-trivial coupling energy, which introduces substantive technical challenges.

Preliminaries
First, we state a few properties of the potential f ε , defined in (2.2). We define Lemma 3.1. The potential f ε satisfies the following properties.
The proof of Lemma 3.1 is contained in Appendix B.
Remark 3.1. The vector fields n ε , m ε are determined by Q ε only up to their sign -Equation (3.5) still holds if we replace n ε by −n ε or m ε by −m ε . Therefore, the unit vector u ε is uniquely determined by Q ε , M ε only up to the sign of its components (u ε ) 1 , (u ε ) 2 . However, the quantity h(u ε ) is is well-defined, irrespective of the choice of the orientations for n ε , m ε , because (3.4). Let u ε be defined as in (3.6). Then, we have In other words, the change of variables (3.6) transforms the functional into a sum of two decoupled terms, which can be studied independently, and a remainder term, which is small compared to the other ones. Before we proceed with the proof of Proposition 3.2, we state some properties of the functions g ε , h defined in (3.7), (3.8) respectively. These properties are elementary, but will be useful later on.
Proof of Proposition 3.2. For simplicity of notation, we omit the subscript ε from all the variables.
Step 1. Let k ∈ {1, 2}. We have M = u 1 n + u 2 m and hence, We raise to the square both sides of (3.10). We apply the identities (3.11) n · ∂ k n = m · ∂ k m = 0, n · ∂ k m + m · ∂ k n = 0, ∂ k n · ∂ k m = 0 which follow by differentiating the orthonormality conditions |n| 2 = |m| 2 = 1, n · m = 0. (In particular, the first identity in (3.10) implies that ∂ k n is parallel to m and ∂ k m is parallel to n, so ∂ k n · ∂ k m = 0.) We obtain We consider the potential term f ε (Q, M). Since (n, m) is an orthonormal basis of R 2 , we have By substituting (3.13) into the definition (2. 2) of f ε , and recalling (3.7), (3.8), we obtain Combining (3.12) with (3.14), we obtain where |G| denotes the area of G. We estimate separately the terms in the right-hand side of (3.15).
Step 2. In view of the identity n ⊗ n + m ⊗ m = I, Equation (3.5) can be written as We differentiate both sides of (3.16) and compute the squared norm of the derivative. Recalling the assumption (3.4), after routine computations we obtain Thanks to (3.17), we can estimate By our assumptions (3.3), (3.4), the L ∞ -norm of u is bounded and the L 2 -norm of ∇Q is of order |log ε| 1/2 at most. Therefore, we obtain Step 3. By Lemma 3.4, the function h has two strict, non-degenerate minima at the points u ± := (±( √ 2β + 1) 1/2 , 0). As a consequence, for any u ∈ R 2 such that |u| ≤ A (where A > 0 is the constant from (3.4)), we must have for some constant C A that depends only on A and β. Then, for any u ∈ R 2 with |u| ≤ A we have for some (possibly) different constant C ′ A , still depending on A and β only. The assumption (3.4) and the property (3.13) guarantee that u satisfies |u| ≤ A almost everywhere in G. Therefore, we can apply (3.21) to estimatê The elementary inequality (x − 1) 2 ≤ (x 2 − 1) 2 , which applies to any x ≥ 0, implieŝ The proposition follows by ( In this section, we prove that the Q * ε -component of the minimisers converges to a limit, up to extraction of subsequences. The results in this section are largely based on the analysis in [11]. Throughout the paper, we denote by (Q * ε , M * ε ) a minimiser of the functional (2.1), subject to the boundary condition (2.3). We recall that the boundary data are of class C 1 and satisfy the assumption (2.4). Routine arguments show that minimisers exist and that they satisfy the Euler-Lagrange equations (2.5)-(2.6).
Proof. Elliptic regularity theory implies, via a bootstrap argument, that (Q * ε , M * ε ) is smooth in the interior of Ω and continuous up to the boundary. Now, we prove (4.1). We take the scalar product of both sides of (2.5) with Q * ε : In a similar way, by taking the scalar product of (2.6) with M * ε , we obtain By adding (4.3) and (4.4), and rearranging terms, we deduce The right-hand side of (4.5) is strictly positive if |Q * ε | 2 + |M * ε | 2 ≥ C, for some (sufficiently large) constant C that depends on β but not on ε. Therefore, (4.1) follows from the maximum principle. The inequality (4.2) follows by [10, Lemma A.1 and Lemma A.2].

Proposition 4.2. Minimisers
where d ∈ Z is the degree of M * ε and C is a constant that depends only on Ω, Q bd , M bd (not on ε).
Proof. We first consider the case d = 1. Consider balls B 1 := B(a 1 , R), B 2 := B(a 2 , R), of centres a 1 , a 2 and radius R > 0, that are mutually disjoint. Since we have assumed that the degree of the boundary datum Q bd is d = 1, there exists a mapQ : Ω \ (B 1 ∪ B 2 ) → N that is smooth (up to the boundary of Ω \ (B 1 ∪ B 2 )), satisfiesQ = Q bd on ∂Ω and has degree 1/2 on ∂B 1 and ∂B 2 . We define a comparison map Q ε as follows: where Q 1 ε , Q 2 ε are given as for some constant C that does not depend on ε. Indeed, sinceQ is regular on Ω \ (B 1 ∪ B 2 ) and takes values in the manifold N , the energy of Q ε on Ω \ (B 1 ∪ B 2 ) is an ε-independent constant, whereas the contribution of Q 1 ε , Q 2 ε is reminiscent of the Ginzburg-Landau functional and can be computed explicitely.
Next, we construct the component M ε . Let Λ be the straight line segment of endpoints a 1 , a 2 . Thanks to Lemma A.3 in Appendix A, there exists a vector fieldM ε ∈ SBV(Ω, R 2 ) such that a.e. in Ω and, moreover, satisfies SM ε = Λ, up to negligible sets. In particular,M ε is smooth in a neighbourhood of ∂Ω. By comparing (2.4) with (4.7), it follows that eitherM ε = M bd on ∂Ω orM ε = −M bd on ∂Ω. Up to a change of sign, we will assume without loss of generality thatM ε = M bd on ∂Ω. In order to define our competitor M ε , we need to regulariseM ε near its jump set. We define For ε small enough, we have M ε =M ε = M bd on ∂Ω. The absolutely continuous part of gradient ∇M ε can be estimated by differentiating both sides of (4.7), by the BV-chain rule; it turns out that |∇M ε | = c |∇Q ε |, up to an (explicit) constant factor c that does not depend on ε.
By explicit computation, we have Finally, we compute the potential. We need to consider three different contributions. At a point due to (4.7) and Lemma 3.
is bounded by a constant that does not depend on ε. Therefore, we have Together, (4.6), (4.8) and (4.9) imply for some constant C that does not depend on ε. The proof in case d ̸ = 1 is similar, except that in the definition ofQ, we need to consider 2 |d| pairwise disjoint balls B 1 , B 2 , . . . B 2|d| , each of them carrying a topological degree of sign(d)/2. The set Λ is defined as a union of segments that connects the centres of the balls B 1 , B 2 , . . . B 2|d| (for instance, a minimal connection -see Appendix A).
The following estimate is well-known estimate in the Ginzburg-Landau literature [26].

Lemma 4.3.
There exists an ε-independent constant C such that Lemma 4.3 is a direct consequence of Theorem 1.1 in [26]. A compactness result for the Q * εcomponent of minimisers can also be obtained by appealing to results in the Ginzburg-Landau theory. Given a (closed) ballB ρ (a) ⊆ Ω such that |Q * ε | ≥ 1/2 on ∂B ρ (a), the map is well-defined and continuous and hence, its topological degree is well-defined as an element of 1 2 Z. We denote the topological degree of ). We recall that d is the degree of the boundary datum, as given in (2.10).
in Ω and a (non-relabelled) subsequence such that the following statement holds. For any σ > 0 sufficiently small there exists ε 0 (σ) > 0 such that, if 0 < ε ≤ ε 0 (σ), then for any j ∈ {1, . . . , 2 |d|}, any k ∈ {1, . . . , K}. Moreover, for any σ sufficiently small and any 0 < ε ≤ ε 0 (σ), there holds where C is a positive constant C that does not depend on ε, σ. Finally, there exists a limit map Q * : Ω → N such that Proof. The analysis of the Q * ε -component can be recast in the classical Ginzburg-Landau setting, by means of a change of variables. We define q * ε : Ω → R 2 as The terms at the right-hand side can be bounded by Proposition 4.2 and Lemma 4.1, respectively. We obtain where C is an ε-independent constant. Moreover, due to the boundary condition (2.3) and (2.4), q * ε restricted to the boundary ∂Ω coincides with an ε-independent map of class C 1 . More precisely, if we identify vectors in R 2 with complex numbers so that M bd is identified with a complex number, M bd = M bd1 + iM bd2 , then a routine computation shows (the square is taken in the sense of complex numbers). In particular, |q * ε | = 1 on ∂Ω and for any ε sufficiently small. Then, (4.13) follows from (4.12) and (4.17), by means of a compactness argument.
In order to complete the proof of Statement (i) in Theorem 2.1, it only remains to show that the convergence Q * ε → Q * is not only weak, but also strong in W 1,p (Ω). The proof of this fact relies on an auxiliary lemma. We consider the function g ε : for some constant C that may depend on the radius r, but not on ε.
The proof of Lemma 4.5 is given in Appendix C.
be an open ball. We have |Q * ε | ≥ 1/2 in B, so we can apply the change of variables described in Section 3. We consider the vector field u * ε : B → R 2 defined as in (3.6) -that is, we write where (n * ε , m * ε ) is an orthonormal basis of eigenvectors for Q ε , and we define ( where the functions g ε and h are defined in (3.7) and (3.8), respectively. (The remainder term R ε , given by Proposition 3.2, tends to zero as ε → 0, due to (3.9) and the energy bound (4.12)). By Lemma 4.4, we know that F ε (Q * ε , M * ε ; B) ≤ C for some constant C that depends on the ball B, but not on ε. By Fubini theorem, and possibly up to extraction of a subsequence, we find a radius r ∈ (R/2, R) such that with C that does not depend on ε. Moreover, without loss of generality we can assume that Thanks to (4.25), we can write where (n ε , m ε ) is an orthonormal basis of eigenvectors for Q ε . We define As we know already that Q * ε ⇀ Q * weakly in W 1,2 (B ′ ) (by Lemma 4.4), from (4.27) we deduce that Q * ε → Q * strongly in W 1,2 (B ′ ) and (4.21) follows. We turn to the proof of (4.22). Let p, q be such that 1 ≤ p < q < 2. Let σ > 0 be a small parameter, and let By the Hölder inequality, we obtain Thanks to Lemma 4.4 and (4.21), we deduce and, as σ may be taken arbitrarily small, (4.22) follows.
Remark 4.1. As a byproduct of the estimate (4.27), we deduce that g ε (Q * ε ) → 0 strongly in We state an additional convergence property for Q * ε , which will be useful later on. We recall that the vector product of two vectors u ∈ R 2 , v ∈ R 2 can be identified with a scalar, In a similar way, we define the vector product of two matrices Q ∈ S 2×2 0 , If q 1 , q 2 (respectively, p 1 , p 2 ) are the columns of Q (respectively, P), then Alternatively, the vector product Q × P can be expressed in terms of the commutator [Q, P] := QP − PQ, as if we take the derivatives in the sense of distributions. If Q is smooth, then J(Q) is the Jacobian determinant of q := ( √ 2Q 11 , √ 2Q 12 ): More generally, for any Q ∈ (L ∞ ∩W 1,1 )(Ω, R 2 ), J(Q) coincides with the distributional Jacobian of q (see e.g. [35] and the references therein).

Proof of Statement (ii): Q * is a canonical harmonic map
Next, we show that Q * is the canonical harmonic map with singularities at (a * 1 , . . . , a * 2|d| ) and boundary datum Q bd , as defined in Section 2. The proof relies on an auxiliary lemma.
Proof. For ease of notation, we drop the subscript ε and the superscript * from all the variables. We consider the Euler-Lagrange equation for Q, Equation (2.5), and take the vector product with Q: We have Now, we consider the Euler-Lagrange equation for M, Equation (2.5), and take the vector product with M: , so (4.37) can be written as Proof. First, we show that Q * satisfies in the sense of distributions in Ω. To this end, we pass to the limit in both sides of Lemma 4.8.
Let p ∈ (1, 2). By Lemma 4.4, we have Q * ε ⇀ Q * weakly in W 1,p (Ω) and, up to extraction of subsequences, pointwise a.e. As Q * ε is bounded in L ∞ (Ω) by Lemma 4.1, Lebesgue's dominated convergence theorem implies that Q * ε → Q * strongly in L q (Ω) for any q < +∞. As a consequence, we have On the other hand, Proposition 4.2 implies Combining ( Therefore, θ * is smooth in G and so is Q * . In case G touches the boundary of Ω, θ * is continuous up to ∂Ω and hence Q * is.

Proof of Statements (iii) and (iv): compactness for M * ε
In this section, we prove a compactness result for the component M * ε of a sequence of minimisers. The proof relies on the change of variables we introduced in Section 3.
We recall that in Lemma 4.4, we found a finite number of points a * 1 , . . . , a * 2|d| , b * 1 , . . . , b * K such that |Q * ε | is uniformly bounded away from zero, except for some small balls of radius σ around these points. Let . Therefore, we are in position to apply the results from Section 3. We define the vector field u * ε : G → R 2 as in (3.6) -that is, we write where (n * ε , m * ε ) is an orthonormal set of eigenvectors for Q * with n * ε ∈ W 1,2 (G, S 1 ), m * ε ∈ W 1,2 (G, S 1 ), and we define The next lemma is key to prove compactness of the sequence u * ε and, hence, of M * ε .

Lemma 4.10.
Let h be the function defined by (3.8). For any simply connected domain G ⊂⊂ where C is a positive constant that depends only on Ω, β and the boundary datum (in particular, it is independent of ε, G).
for some constant C that depends only on Ω and the boundary datum Q bd . The results in [34,43] extend to our setting due to change of variables Q * ε → q * ε , given by (4.14). The coefficient 2π |d| in the right-hand side of (4.46) depends on this change of variables, which transforms the boundary condition of degree d for Q * ε into a boundary condition of degree 2d for q * ε -see (4.16 for some constant C that depends only on the domain and the boundary data. Now, we apply Proposition 3.2: We have used (3.9) and the elementary inequality ab ≤ a 2 /2 + b 2 /2 to estimate the remainder term R ε . From (4.49), we obtain a.e. on Ω.
The vector field M * is well-defined and does not depend on the choice of the orientation for n * ε , m * ε (so long as the orientation is chosen consistently as ε → 0, in such a way that (4.53) is satisfied). Indeed, if we replace n * ε by −n * ε , then also (u * ε ) 1 will change its sign and the product at the right-hand side of (4.57) will remain unaffected. Therefore, by letting G vary in Ω \ {a * 1 , . . . , a * 2|d| , b * 1 , . . . , b * K }, we can define M * almost everywhere in Ω. An explicit computation, based on (4.54) and (4.57), shows that M * satisfies (4.51) and (4.52). Moreover, due to (4.53) and (4.55), we have M * ε → M * a.e. in G. As the sequence M ε is uniformly bounded in L ∞ (Ω) (by Lemma 4.1), Lebesgue's dominated convergence theorem implies that M * ε → M * in L p (Ω) for any p < +∞.
As we have seen, u * ∈ SBV(G, R 2 ) for any G ⊂⊂ Ω \ {a * 1 , . . . , a * 2|d| , b * 1 , . . . , b * K }. Therefore, by applying the BV-chain rule (see e.g. [3,Theorem 3.96]) to (4.59), and letting G vary, we obtain Moreover, we claim that and hence, due to Lemma 4.4. The total variation of the jump part of DM * is uniformly bounded, too, because of (4.56) (the constant at the right-hand side of (4.56) does not depend on G, so we may take the limit as G ↘ Ω). Then, (4.59) follows.
In order to complete the proof, it only remains to show that M ∈ SBV(Ω, R 2 ). Let φ ∈ C ∞ c (Ω) be a test function, and let σ > 0 be fixed. We define We choose a smooth cut-off function ψ σ such that 0 ≤ ψ σ ≤ 1 in Ω, ψ σ = 0 in Ω \ U σ , ψ σ = 1 in a neighbourhood of each point a 1 , . . . , a 2|d| , b 1 , . . . , b K , and ∥∇ψ σ ∥ L ∞ (Ω) ≤ Cσ for some constant C that does not depend on σ. Then, for j ∈ {1, 2}, we havê We bound the first term in the right-hand side by applying (4.59). To estimate the second term, we observe that the integrand is bounded and supported in U σ . Therefore, we obtain By taking the limit as σ → 0, we deduce that M * ∈ BV(Ω, R 2 ). In fact, we must have M * ∈ SBV(Ω, R 2 ), because the Cantor part of DM * cannot be supported on a finite number of points, a 1 , . . . , a 2|d| , b 1 , . . . , b K . This completes the proof.
We conclude this section by stating a regularity property of M * . We recall that a harmonic map M on a domain U ⊆ R 2 with values in a circle of radius R > 0 is a map that can be written in the form M = (R cos ϕ, R sin ϕ) for some harmonic function ϕ : U → R. Let S M * be the closure of the jump set of M * . Proof. Let B ⊆ Ω be an open ball that does not intersect S M * nor {a * 1 , . . . , a * 2|d| , b * 1 , . . . , b * K }. Then, we have M * ∈ W 1,2 (B, R 2 ), by construction (see, in particular, (4.53) and (4.57)). By lifting results (see e.g. [13,12,8]), M * can be written in the form M * = ( √ 2β +1) 1/2 (cos ϕ * , sin ϕ * ), for some scalar function ϕ * ∈ W 1,2 (B, R). On the other hand, the condition (4.52) shows that ϕ * is uniquely determined by Q * , up to constant multiples of π. In particular, we must have ϕ * = θ * /2 + kπ, where θ * is the function given by (4.42) and k is a constant. Then, Equation (4.43) implies that −∆ϕ * = 0 in B and hence, M * is a harmonic map on B with values in the circle of radius (2 √ β + 1) 1/2 . Now, let B be an open ball that does not intersect S M * nor {a * 1 , . . . , a * 2|d| }, although it may contain one of the points b k . Say, for simplicity, that B contains exactly one of the points b k . We claim that M * is harmonic in B, too. Indeed, since b k is a singularity of degree zero (see (4.11)), we can repeat the arguments above and write M * = ( √ 2β + 1) 1/2 (cos ϕ * , sin ϕ * ) in B \ {b k }, for some harmonic function ϕ * : B \ {b k } → R. By the chain rule, |∇ϕ * | coincides with |∇Q * | up to a constant factor (see (4.60)). The map Q * is smooth in a neighbourhood of b k , because it is canonical harmonic with singularities at {a 1 , . . . , a 2|d| }. Therefore, ∇ϕ * is bounded in B \ {b k }. As a consequence, b k is a removable singularity for ϕ * and, by possibly modifying the value of M * at b k , M * is harmonic in B.
To conclude the proof, it only remains to show that the points {a * 1 , . . . , a * 2|d| } are contained in S M * . If any of the points a j did not belong to S M * , then M * would be locally harmonic (and hence, smooth) in a sufficiently small neighbourhood of a j , except at the point a j . This is impossible, because a j is a non-orientable singularity of Q * (see (4.11)) and there cannot be a map M * that satisfies (4.51), (4.52) and is continuous in a punctured neighbourhood of a j . Therefore, a j ∈ S M * .

Proof of Statements (v) and (vi): sharp energy estimates
In this section, we complete the proof of Theorem 2.1, by describing the structure of the jump set of M * and characterising the optimal position of the defects of Q * (in case the domain Ω is convex). As a byproduct of our arguments, we will also show a refined energy estimate for the minimisers (Q * ε , M * ε ), i.e. Proposition 4.13 below. First, we set some notations. We let  (a 1 , . . . , a 2|d| ) where W, L are, respectively, the Ginzburg-Landau renormalised energy (defined in (2.14)) and the length of a minimal connection (defined in (2.15)). We also recall the definition of the Ginzburg-Landau core energy, which was introduced in [11]. Let B 1 ⊆ R 2 be the unit disk. For any ε > 0, let It can be proved (see [11,Lemma III.3]) that the function ε → γ(ε) − π |log ε| is finite in (0, 1) and non-decreasing. Therefore, the limit (4.64) γ * := lim ε→0 (γ(ε) − π |log ε|) > 0 exists and is finite. The number γ * is the so-called core energy. In this section, we will prove the following result: We will prove the lower and upper inequality in (4.65) separately. From now on, we alwasy assume that the domain Ω is convex.

Sharp lower bounds for the energy of minimisers
The aim of this section is to prove a sharp lower bound for F ε (Q * ε , M * ε ). We know from previous results (Lemma 4.4, Proposition 4.11), that, up to extraction of a subsequence, we have Q Due to Lemma 4.3, we may further assume that (4.66) |Q * ε | − 1 ε ⇀ ξ * weakly in L 2 (Ω).
The length of the jump set S M * can be further bounded from below, in terms of the singular points a * 1 , . . . , a * 2|d| . We recall from Section 2 that a connection for a * 1 , . . . , a * 2|d| is a finite collection of straight line segments L 1 , . . . , L |d| that connects the points a * j in pairs, and that L(a * 1 , . . . , a * 2|d| ) is the minimal length of a connection for the points a * j (see (2.15)). Given two sets A, B, we denote by A∆B their symmetric difference, i.e. A∆B :  We will give the proof of Proposition 4.15 in Appendix A. Here, instead, we focus on the proof of Proposition 4.14.
Proof. We make a change of variable, as introduced in Section 3. This is possible, because we have assumed that Ω is simply connected. Let u * ε be the vector field defined in (3.6). By Proposition 3.2, we have and the remainder term R ε satisfies Lemma 4.10 implies that R ε → 0 as ε → 0. We estimate separately the other terms in the right-hand side of (4.71). The weak convergence Q * ε ⇀ Q * in W 1,2 (G) implies We claim that Let δ > 0 be a small parameter. By Lemma 4.4, we have |Q * ε | → |Q * | = 1 a.e. in Ω and hence, ζ ε → 0 a.e. in G. Therefore, by the Severini-Egoroff theorem, there exists a Borel setG ⊆ G such that |G \G| ≤ δ and ζ ε → 0 uniformly inG as ε → 0. Now, we havê The integral of ξ 2 ε ζ * onG tends to zero, because ξ ε is bounded in L 2 (G) (by Lemma 4.3) and ζ ε → 0 uniformly inG. As ξ ε ⇀ ξ * weakly in L 2 (G) (see (4.66)), we deduce The area of G \G may be taken arbitrarily small, so (4.74) follows.
Finally, for the term in u * ε , we apply classical Γ-convergence results for the vectorial Modica-Mortola functional (see e.g. [7,28]), as well as Lemma 4.16: where u * is the L 1 (G)-limit of the sequence u * ε , as in (4.55). By (4.57), we have S u * = S M * and hence, Combining where γ * is the constant given by (4.64) and C is a constant that does not depend on ε, σ.
Proof. Take σ is so small that the ball B σ (a * j ) does not contain any other singular point a * k , with k ̸ = j. We consider the function J(Q * ε ) defined in (4.30). By Lemma 4.7, we have ). Then, we can apply pre-existing Γ-convergence results for the Ginzburg-Landau functional -for instance, [1, Theorem 5.3]. We obtain a (sharp) lower bound for the Ginzburg-Landau energy of Q * ε : On the other hand, Lemma 3.1 gives As M * ε is uniformly bounded in L ∞ (Ω) (by Lemma 4.1), we obtain via the Hölder inequality The constant C here depends only on β. Finally, Lemma 4.3 implies Combining (4.78) with (4.79), the lemma follows.
We can now complete the proof of Proposition 4.14.
Proof of Proposition 4.14. Let σ > 0 be small enough that the balls B σ (a * j ), B σ (b * k ) are pairwise disjoint. We define We construct open sets G 1 , . . . , G N with the following properties: (i) the sets G i are pairwise disjoint; (ii) their closures, G i , cover all of Ω σ ; (iii) each G j is simply connected; (iv) H 1 (S M * ∩ ∂G i ∩ Ω σ ) = 0 for any j.
For instance, we can partition Ω σ by considering a grid, consisting of finitely many vertical and horizontal lines. Since H 1 (S M * ) < +∞ by Proposition 4.11, we have for all but countably many values of c ∈ R, d ∈ R. We choose numbers that satisfy (4.80), in such a way that Ω ⊆ (c 0 , c N 1 ) × (d 0 , d N 2 ). For a suitable choice of c h , d ℓ , we can make sure that no ball B σ (a * j ) or B σ (b * k ) is entirely contained in a rectangle of the form (c h , c h+1 ) × (d ℓ , d ℓ+1 ), and that any rectangle (c h , c h+1 ) × (d ℓ , d ℓ+1 ) intersects at most one of the balls. Then, the sets G h,ℓ := ((c h , c h+1 ) × (d ℓ , d ℓ+1 )) ∩ Ω σ are all simply connected and satisfy the properties (i)-(iv) above. We relabel the G h,ℓ 's as G i . We apply Lemma 4.17 on each G i , then sum over all the indices i. We obtain On the other hand, Lemma 4.18 implies for any j ∈ {1, . . . , 2 |d|}. Combining (4.81) with (4.82), we obtain By Proposition 4.9, Q * is the canonical harmonic map with singularities at (a * 1 , . . . , a * 2|d| ) and boundary datum Q bd . Then, we can write the right-hand side of (4.83) in terms of the renormalised energy, W, defined in (2.14). First, we observe that ). Then, from (4.83), (4.84) and (2.14) we deduce Now we pass to the limit in both sides of (4.85), first as ε → 0, then as σ → 0. The proposition follows.
The proof of Proposition 4.19 is based on a rather explicit construction. For the component Q ε , we follow classical arguments from the Ginzburg-Landau literature (see e.g. [11,1]), with minor modifications. For the component M ε , we first construct a vector fieldM ε : Ω → R 2 of constant norm, such thatM ε (x) is an eigenvector of Q ε (x) at each point x ∈ Ω. As Q ε has non-orientable singularities at the points a j , there is no smooth vector fieldM ε with this property. However, we can construct a BV-vector fieldM ε , which jumps along finitely many line segments that join the points a j along a minimal connection (see Appendix A). Then, we define M ε by regularisingM ε in a small neighbourhood of the jump set. The regularisation procedure is reminiscent of the optimal profile problem for the Modica-Mortola functional [42].

Numerics
In this section, we numerically compute some stable critical points of the ferronematic free energy, on square domains with topologically non-trivial Dirichlet boundary conditions for Q and M. These numerical results do not directly support our main results on global energy minimizers of (2.1) in the ϵ → 0 limit, since the numerically computed critical points need not be global energy minimizers, and we expect multiple local and global energy minimizers of (2.1) for ϵ > 0. Instead of solving the Euler-Lagrange equations directly, we solve a L 2 -gradient flow associated with the effective re-scaled free energy for ferronematics (2.1), given by Here η 1 > 0 and η 2 > 0 are arbitrary friction coefficients. Due to limited physical data, we do not comment on physically relevant values of ε, β and the friction coefficients. The system of L 2 -gradient flow equations for Q 11 , Q 12 and the components, M 1 , M 2 of the magnetisation vector, can be written as The stationary time-independent or equilibrium solutions of the L 2 -gradient flow satisfy the original Euler-Lagrange equations of (2.1). For non-convex free energies as in (2.1), there are multiple critical points, with many of them being unstable saddle points [47]. One can efficiently compute stable critical points of such free energies by considering the L 2 -gradient flow associated with the non-convex free energies and these gradient flows converge to a stable critical point, for a given initial condition, thus avoiding the unstable saddle points. From a numerical standpoint, the L 2 -gradient flow can be more straightforward to solve than the nonlinear coupled Euler-Lagrange equations, primarily due to the inclusion of time relaxation in the L 2 -gradient flow.
In the following simulations, we take η 1 = 1 and η 2 = ε and do not offer rigorous justifications for these choices, except as numerical experiments to qualitatively support out theoretical results. We impose the continuous degree +k boundary condition and atan2(y, x) is the 2-argument arctangent that computes the principal value of the argument function applied to the complex number x + iy. So −π ≤ atan2(y, x) ≤ π. For example, if x > 0, then atan2(y, x) = arctan y x . The initial condition is prescribed to be We solve the L 2 -gradient flow equation using standard central finite difference methods [33]. For the temporal discretization, we employ a second-order Crank-Nicolson method [33]. The grid size and temporal step size are denoted by h and τ , respectively. In all our computations, we set h = 1/50 and τ = 1/1000.
In Figure 1, we plot the dynamical evolution of the solutions of the gradient flow equations, for k = 1 boundary conditions, with the initial condition (5.5). The time-dependent solutions converge for t ≥ 1, and we treat the numerical solution at t = 1 to the converged equilibrium state. We cannot conclusively argue that the converged solution is an energy minimizer but it is locally stable, the converged Q-profile has two non-orientable defects and the corresponding Mprofile has a jump set composed of a straight line connecting the nematic defect pair, consistent with our theoretical results on global energy minimizers. We consider two different values of ε and it is clear that the Q-defects and the jump set in M become more localised as ε becomes smaller, as expected from the theoretical-results. We have also investigated the effects of β on the converged solutions -the defects become closer as β increases. This is expected since the  cost of the minimal connection between the nematic defects increases as β increases, and hence the shorter connections require the defects to be closer to each other (at least in a pairwise sense).
In Figure 2, we plot the dynamical evolution of the solutions of the gradient flow equations, for k = 2 boundary conditions, with the initial condition (5.5), and we treat the numerical solution at t = 1 to be the converged equilibrium state. Again, the converged solution is locally stable, the Q-profile has four non-orientable defects,the M-profile has two distinct jump sets connecting two pairs of non-orientable nematic defects, and the jump sets are indeed approximately straight lines. Smaller values of ϵ correspond to the sharp interface limit which induces more localised defects for Q, straighter line defects for M and larger values of β push the defects closer together, all in qualitative agreement with our theoretical results. where the white bars represent nematic field n (the eigenvector of Q associated with the largest eigenvalue) and the color represents tr Q 2 = 2(Q 2 11 + Q 2 12 ); the magnetic configuration is shown in the right panel, where the white bars represent magnetic field M and the color represents Theorem 2.1 is restricted to global minimizers of (2.1) in the ϵ → 0 limit, but the numerical illustrations in Figures 1 and 2 suggest that Theorem 2.1 may also partially apply to local energy minimizers of (2.1). In other words, locally energy minimizing pairs, (Q ε , M ε ), may also converge to a pair (Q * , M * ), for which Q * is a canonical harmonic map with non-orientable point defects and M * has a jump set connecting the non-orientable point defects of Q * , with the location of the defects being prescribed by the critical point(s) of the normalization energy in Theorem 2.1. The numerical illustrations in Figures 1 and 2 cannot be directly related to Theorem 2.1, since we have only considered two small and non-zero values of ϵ and for a fixed β > 0, there maybe multiple local and global energy minimizers with different jump sets in M i.e. different choices of the minimal connection of equal length, or different connections of different lengths between the nematic defect pairs. For example, it is conceivable that a locally stable M-profile also connects the nematic defects by means of straight lines, but this connection is not minimal. There may also be non energy-minimising critical points with orientable point defects in M tailored by the nonorientable nematic defects. Similarly, there may be non energy-minimising critical points with non-orientable and orientable nematic defects, whose locations are not minimisers but critical points of the modified renormalised energy in Theorem 2.1. We defer these interesting questions to future work.

Conclusions
We study a simplified model for ferronematics in two-dimensional domains, with Dirichlet boundary conditions, building on previous work in [14]. The model is only valid for dilute ferronematic suspensions and we do not expect quantitative agreement with experiments. Further, the experimentally relevant choices for the boundary conditions for M are not well established and our methods can be adapted to other choices of boundary conditions e.g. Neumann conditions for the magnetisation vector. Similarly, it is not clear if topologically non-trivial Dirichlet conditions can be imposed on the nematic directors, for physically relevant experimental scenarios. Having said that, our model problem is a fascinating mathematical problem because of the tremendous complexity of ferronematic solution landscapes, the multiplicity of the energy minimizers and non energy-minimizing critical points, and the multitude of admissible coupled defect profiles for the nematic and magnetic profiles. There are several forward research directions, some of which could facilitate experimental observations of the theoretically predicted morphologies in this manuscript. For example, one could study the experimentally relevant generalisation of our model problem with Dirichlet conditions for Q and Neumann conditions for M, or study different asymptotic limits of the ferronematic free energy in (1.1), a prime candidate being the ε → 0 limit for fixed ξ and c 0 (independent of ε). This limit, although relevant for dilute suspensions, would significantly change the vacuum manifold N in the ε → 0 limit. In fact, we expect to observe stable point defects in the energy-minimizing M-profiles for this limit, where ξ and c 0 are independent of ε, as ε → 0. Further, there is the interesting question of how this ferronematic model can be generalised to non-dilute suspensions or to propose a catalogue of magneto-nematic coupling energies for different kinds of MNP-MNP interactions and MNP-NLC interactions. The physics of ferronematics is complex, and it is challenging to translate the physics to tractable mathematical problems with multiple order parameters, and we hope that our work is solid progress in this direction with bright interdisciplinary prospects.
Taxonomy: GC, BS and AM conceived the project based on a model developed by AM and her ex-collaborators. GC and BS led the analysis, followed by AM. YW performed the numerical simulations, as advised by AM and GC. All authors contributed to the scientific writing.
GC gratefully acknowledges the hospitality provided by the University Federico II (Naples) under the PRIN project 2017TEXA3H, and AM, BS and GC gratefully acknowledge support from the Erwin Schrodinger Institute in Vienna in December 2019, all of which facilitated this collaboration. AM acknowledges support from the Leverhulme Trust and the University of Strathclyde New Professor's Fund. We thank the referee for their careful reading of the manuscript and comments.
Data Availability Statement: Data sharing not applicable to this article as no datasets were generated or analysed during the current study.
Conflict of interest statement. The authors have no competing interests to declare that are relevant to the content of this article.

A Lifting of a map with non-orientable singularities
The aim of this section is to prove Proposition 4.15. We reformulate the problem in a slightly more general setting.
Let a ∈ R 2 , and let Q ∈ W 1,2 loc (R 2 \ {a}, N ). By Fubini theorem and Sobolev embedding, the restriction of Q on the circle ∂B ρ (a) is well-defined and continuous for a.e. ρ > 0. Therefore, it makes sense to define the topological degree of Q on ∂B ρ (a) as an half-integer, deg(Q, a) ∈ 1 2 Z. As the notation suggests, the degree is independent of the choice of ρ: for a.e. 0 < ρ 1 < ρ 2 , the degrees of Q on ∂B ρ 1 (a) and ∂B ρ 2 (a) are the same. If Q is smooth, this is a consequence of the homotopy lifting property; for more general Q ∈ W 1,2 loc (R 2 \ {a}, N ), this follows from an approximation argument (based on [44,Proposition p. 267]). We will say that a is a nonorientable singularity of Q if deg(Q, a) ∈ 1 2 Z \ Z. Given an open set Ω ⊆ R 2 , a map Q : Ω → N and a unit vector field M : Ω → S 1 , we say that M is a lifting for Q if Any map Q ∈ BV(Ω, N ) admits a lifting M ∈ BV(Ω, S 1 ) (see e.g. [32]). The vector field M * given by Theorem 2.1 is not a lifting of Q * , according to the definition above, because |M * | ̸ = 1. However, |M * | is still a positive constant (see Proposition 4.11), so we can construct a lifting of unit-norm simply by rescaling.
We focus on properties of the lifting for Q-tensors of a particular form, namely, we assume that Q has an even number of non orientable singularities at distinct points a 1 , . . . , a 2d . We recall that a connection for {a 1 , . . . , a 2d } as a finite collection of straight line segments {L 1 , . . . , L d }, with endpoints in {a 1 , . . . , a 2d }, such that each a i is an endpoint of one of the segments L j . We recall that  A minimal connection for {a 1 , . . . , a 2d } is a connection that attains the minimum in the righthand side of (A.2). Given two sets A, B, we denote their symmetric difference as A∆B := (A \ B) ∪ (B \ A). Proof. Suppose, towards a contradiction, that {L 1 , . . . , L d } is a minimal connection with L 1 ∩ L 2 ̸ = ∅. The intersection L 1 ∩L 2 must be either a non-degenerate sub-segment of both L 1 and L 2 or a point. If L 1 ∩ L 2 is non-degenerate, then (L 1 ∪ L 2 ) \ (L 1 ∩ L 2 ) can be written as the disjoint union of two straight line segments, K 1 and K 2 , and This contradicts the minimality of {L 1 , . . . , L d }. Now, suppose that L 1 ∩ L 2 is a point. By the pigeon-hole principle, L 1 ∩L 2 cannot be an endpoint for either L 1 or L 2 . Say, for instance, that L 1 is the segment of endpoints a 1 , a 2 , while L 2 is the segment of endpoints a 3 , a 4 . Let H 1 , H 2 be the segments of endpoints (a 1 , a 3 ), (a 2 , a 4 ) respectively. Then, by the triangular inequality, which contradicts again the minimality of {L 1 , . . . , L d }.

Lemma A.3.
Let Ω ⊆ R 2 be a bounded, convex domain and let a 1 , . . . , a 2d be distinct points in Ω. Let Q ∈ W 1,1 (Ω, N ) ∩ W 1,2 loc (Ω \ {a 1 , . . . , a 2d }, N ) be a map with a non-orientable singularity at each a j . If {L 1 , . . . , L d } is a minimal connection for {a 1 , . . . , a 2d }, then there exists a lifting M * ∈ SBV(Ω, Proof. For any ρ > 0 and j ∈ {1, . . . , d}, we define and Since Ω is convex, L j ⊆ Ω for any j and hence, U j,ρ ⊆ Ω for any j and ρ small enough. Each U j ρ is a simply connected domain with piecewise smooth boundary. Moreover, for ρ fixed and small, the sets U j,ρ are pairwise disjoint, because the L j 's are pairwise disjoint (Lemma A.2). The trace of Q on ∂U j,ρ is orientable, because ∂U j,ρ contains exactly two non-orientable singularities of Q.
By construction, we have S M * ⊆ ∪ j L j . Therefore, it only remains to prove that S M * contains H 1 -almost all of ∪ j L j . Consider, for instance, the segment L 1 ; up to a rotation and traslation, we can assume that L 1 = [0, b]×{0} for some b > 0. Given a small parameter ρ > 0 and t ∈ (0, b), we define K ρ,t := (−ρ, t) × (−ρ, ρ). Fubini theorem implies that, for a.e. ρ and t, Q restricted to ∂K ρ,t belongs to W 1,2 (∂K ρ,t , N ) and hence, by Sobolev embedding, is continuous. Since the segments L j are pairwise disjoint by Lemma A.2, for ρ small enough there is exactly one nonorientable singularity of Q inside K ρ,t . Therefore, Q is non-orientable on ∂K ρ,t for a.e. t ∈ (0, b) and a.e. ρ > 0 small enough; in particular, there is no continuous lifting of Q on ∂K ρ,t . Since M * is continuous on ∂K ρ,t \ L 1 for a.e. ρ and t, we conclude that S M * contains H 1 -almost all of L 1 .
Given a countably 1-rectifiable set Σ ⊆ R 2 and a H 1 -measurable unit vector field τ : Σ → S 1 , we say that τ is an orientation for Σ if τ (x) spans the (approximate) tangent line of Σ at x, for H 1 -a.e. x ∈ Σ. In case Σ is the jump set of an SBV-map M, τ : S M → S 1 is an orientation for S M if and only if τ (x) · ν M (x) = 0 for H 1 -a.e. x ∈ S M . Lemma A.4. Let Ω ⊆ R 2 be a bounded, convex domain and let a 1 , . . . , a 2d be distinct points in Ω. Let Q ∈ W 1,1 (Ω, N ) ∩ W 1,2 loc (Ω \ {a 1 , . . . , a 2d }, N ) be a map with a non-orientable singularity at each a j . Let {L 1 , . . . , L d } be a minimal connection for {a 1 , . . . , a 2d }. Up to relabelling, we assume that L j is the segment of endpoints a 2j−1 , a 2j , for any j ∈ {1, . . . , d}. Let M ∈ SBV(Ω, S 1 ) be a lifting for Q such that S M ⊂⊂ Ω. Then, there exist H 1 -measurable sets T j ⊆ L j and an orientation τ M for S M such that, for any φ ∈ C ∞ c (R 2 ), there holdŝ Proof. Let M * ∈ SBV(Ω, S 1 ) be the lifting of Q given by Lemma A.3. By construction, S M * coincides with ∪ j L j ⊂⊂ Ω up to H 1 -negligible sets. Since we have assumed that S M ⊂⊂ Ω, there exists a neighbourhood U ⊆ Ω of ∂Ω in Ω such that M ∈ W 1,1 (U, S 1 ), M * ∈ W 1,1 (U, S 1 ). A map that belongs to W 1,1 (U, N ) has at most two different liftings in W 1,1 (U, S 1 ), which differ only for the sign [ where ∂ * A is the reduced boundary of A and τ A is an orientation for ∂ * A. Up to H 1 -negligible sets, ∂ * A coincides with S M·M * (see e.g. (1 + τ A · τ j )∇φ · τ j dH 1 = 0, (A.5) On H 1 -almost all of L j \ S M , both τ j and τ A are tangent to L j . Therefore, for H 1 -a.e. x ∈ L j \ S M we have τ A (x) · τ j (x) ∈ {−1, 1}. If we define T j := {x ∈ L j \ S M : τ A (x) · τ j (x) = 1}, then the lemma follows from (A.5).
Lemma A.4 can be reformulated in terms of currents. We recall a few basic definitions in the theory of currents, because they will be useful to complete the proof of Proposition A.1. Actually, we will only work with currents of dimension 0 or 1. We refer to, e.g., [27,45] for more details.
A 0-dimensional current, or 0-current, in R 2 is just a distribution on R 2 , i.e. an element of the topological dual of C ∞ c (R 2 ) (where C ∞ c (R 2 ) is given a suitable topology). A 1-dimensional current, or 1-current, in R 2 is an element of the topological dual of C ∞ c (R 2 ; (R 2 ) ′ ), where (R 2 ) ′ denotes the dual of R 2 and C ∞ c (R 2 ; (R 2 ) ′ ) is given a suitable topology, in much the same way as C ∞ c (R 2 ). In other words, a 1-dimensional current is an R 2 -valued distribution. The boundary of a 1-current T is the 0-current ∂T defined by ⟨∂T, φ⟩ := ⟨T, dφ⟩ for any φ ∈ C ∞ c (R 2 ).
The set of rectifiable 0-currents, respectively rectifiable 1-currents, is denoted by R 0 (R 2 ), respectively R 1 (R 2 ). Given a Lipschitz, injective map f : [0, 1] → R 2 , we denote by f # I the rectifiable 1-current carried by f ([0, 1]), with unit multiplicity and orientation given by f ′ . The mass of f # I is the length of the curve parametrised by f and ∂(f # I) = δ f (1) − δ f (0) ; in particular, ∂(f # I) = 0 if f (1) = f (0). The assumption that f is injective can be relaxed; for instance, if the curve parametrised by f has only a finite number of self-intersections, then f # I is still well-defined and the properties above remain valid.

Lemma A.5.
Let Ω, Q be as above. Let M ∈ SBV(Ω, S 1 ) be a lifting of Q with S M ⊂⊂ Ω. Then, there exist countably may Lipschitz functions f j : [0, 1] → R 2 , with finitely many selfintersections, a rectifiable 1-current R ∈ R 1 (R 2 ) and a permutation σ of the indices {1, . . . , 2d} such that the following properties hold: with P := ∂(R + Q) (and Q as in (A.9)). The current 2P = ∂T − 2d i=1 δ a i is rectifiable, so M(P ) < +∞. Moreover, P isthe boundary of a rectifiable 1-current. Then, Federer's closure theorem [27, 4.2.16] implies that P itself is rectifiable. As a consequence, we can re-write (A.14) as for some integers n k and some distinct points b k ∈ R 2 . By applying [27, 4.2.25], we find countably many Lipschitz, injective maps g j : [0, 1] → R such that For any j, we have either ∂(g j,# I) = 0 (if g j,# (1) = g j,# (0)) or M(∂(g j,# I)) = 2 (otherwise). Therefore, by (A. 16), there are only finitely many indices j such that g j,# (1) ̸ = g j,# (0). Up to a relabelling of the g j 's, we assume that there is an integer q such that g j,# (1) ̸ = g j,# (0) if and only if j ≤ q. Now, the problem reduces to a combinatorial, or graph-theoretical, one. We consider the finite (multi-)graph G whose edges are the curves parametrised by g 1 , . . . , g q , and whose vertices are the endpoints of such curves. There can be two or more edges that join the same pair of vertices. However, we can disregard the orientation of the edges: changing the orientation of the curve parametrised by g j corresponds to passing from the current g j,# I to the current −g j,# I; the difference g j,# I − (−g j,# I) = 2g j,# I can be absorbed into the term 2R that appears in (A.10).
We would like to partition the set of edges of G into d disjoint subsets E 1 , . . . E d , where each E j is a trail (i.e., a sequence of distinct edges such that each edge is adjacent to the next one) and, for a suitable permutation σ of {1, . . . , 2d}, the trail E j connects a σ(2j−1) with a σ(2j) . If we do so, then we can define f j : [0, 1] → R 2 for j ∈ {1, . . . , d} as a Lipschitz map that parameterises the trail E j , with suitable orientations of each edge; for j ≥ d + 1, we define f j := g q+j−d . With this choice of f j , the lemma follows. It is possible to find E 1 , . . . E d as required because the graph G has the following property: any a i is an endpoint of an odd number of edges of G ; conversely, any vertex of G other than the a i 's is an endpoint of an even number of edges of G . This property follows from (A.15). Then, we can construct E 1 , . . . E d by reasoning along the lines of, e.g., [16,Theorem 12].

Proof of Proposition
The equality can only be attained if there are exactly d maps f j and each of them parametrises a straight line segment.

B Properties of f ε
The aim of this section is to prove Lemma 3.1. We first of all, we characterise the zero-set of the potential f ε , in terms of the (unique) solution to an algebraic system depending on ε and β.