Emergent tri-criticality in magnetic metamaterials

Metallic discs engineered on the 100 nm scale have an internal magnetic texture which varies from a fully magnetized state to a vortex state with zero moment. The interplay between this internal structure and the inter-disc interactions is studied in magnetic metamaterials made of square arrays of the magnetic discs. The texture is modeled by a mesospin of varying length with O(2) symmetry and the inter-disc interaction by a nearest neighbour coupling between mesospins. The thermodynamic properties of the model are studied numerically and an ordering transition is found which varies from Kosterlitz-Thouless to first order via an apparent tri-critical point. The effective critical exponent characterising the finite size magnetization evolves from the value for the 2D-XY model to less than half this value at the tri-critical point. The consequences for future experiments both in and out of equilibrium are discussed.


I. INTRODUCTION
Universality, phase transitions and emergent magnetic properties are examples of phenomena that have recently been explored in metamaterials [1][2][3][4][5][6][7][8][9][10] .The ability to choose and investigate the effect of a single parameter, such as spin or spatial dimensionality [2][3][4]9 , as well as the possibility to directly observe individual magnetic elements has been a major impetus in this context [11][12][13][14][15][16][17][18] . Tese are indicators that nano-engineered materials can, in analogy with cold atom systems, become simulators of model many body problems, offering clear advantages over traditional condensed matter systems operating on the atomic scale.In this regard, metamaterials made up of magnetostatically interacting mesoscale islands, or mesospins, are highly attractive.The multi-scale nature of the experimental set up allows for the emergence of new degrees of freedom from the internal spin textures, giving rise to rich behavior beyond that of standard magnetic models [19][20][21][22][23][24][25] .Mesoscopic arrays of circular magnetic islands show a vast ensemble of internal magnetic textures, which vary strongly with the local environment.One of the most characteristic textures is a vortex which can progressively unwind from a state with zero magnetic moment into a collinear state with maximal moment 26 .
The change from vortex to collinear states was shown to be driven by a competition between inter and intraisland interactions so that in an emergent description the interactions between mesospins are self-consistently modified by the collective environment 26 .As a consequence, mesospin ordering occurs via an emergent transition that depends on interactions at both the meso and atomic scales.In the initial experiments the transition was shown to be kinetic in nature, although a route towards true thermodynamic phase transitions was also identified.This suggests that the interplay between collective and internal energy scales could indeed open the door to phases and phase transitions that are not at present obtainable in atomic systems.
In this paper we present a simple model which cap-tures the essence of the interplay between the meso and the atomic length scales.The inter-island interactions are allowed to influence the net moment of the elements, which provides the coupling between the length scales involved.This leaves an XY spin model with an internal degree of freedom; the spin length, which can vary with an associated energy scale.We find that, as a function of this internal energy scale, the magnetic phase transition evolves abruptly from Kosterlitz-Thouless (KT) to 1 st order at a point showing a remarkable resemblance to a tri-critical point.That is, despite the absence of true long range magnetic order and the continuously varying spin length, the phase diagram closely resembles that of the S = 1 Blume-Capel model 27,28 .In this case, we observe an effective critical exponent, relevant for finite size systems, that varies continuously from that observed in 2D XY magnets of finite size, β ≈ 0.23, towards a value characteristic of a tri-critical point.

II. MODELING THE MAGNETIC METAMATERIAL
The total magnetic moment of interacting circular islands depends on the intrinsic material properties: internal magnetic texture, geometry, size and separation.For instance, above the inherent ordering temperature of the material there is no magnetic order at any length scale, while below that temperature, both thermal fluctuations and magnetic texture on the meso-scale are essential elements dictating the moment of the islands 26 .As examples a vortex state is a magnetic texture with a zero net in-plane moment, while a collinear inner state of the islands yields the largest net moment.The single vortex state can be characterized by two observables: a continuous change of the moment within the disc and a shift of the position of the vortex core, as illustrated in the top schematic of Fig. 1.In the vortex state with a zero net moment, the vortex core sits at the center of the disc.The moment on the disc increases from zero as the vortex unwinds and approaches the disc edge.It is 1.The black dots represent results from micromagnetic simulations for an isolated island with Ec/Ev ≈ 1.1.The solid (colored) lines represent the energy obtained from S (Eq.2) for an isolated island for two different values of Ec/Ev.The red dashed lines represents the influence of parallel or antiparallel collinear neighbours for Ec/Ev = 1.1.
In this paper we retain the variable moment length as the main manifestation of the evolving magnetic texture, leaving the effects of the evolution in the vorticity for future work.The total moment on disc i then becomes an in-plane vector M i , allowing for the definition of an in-plane, dimensionless mesospin vector of length where M max is the magnitude of the saturated total moment.The islands can thus be viewed as mesospins whose variable length depends on the internal spin texture of the discs.It varies continuously between zero and one, depending on the competition between internal and many body energy scales and a suitable model must include both these features.The energy scale, E, associated with the variation of r of an isolated disc has been studied in detail in previous work 26 .The internal magnetic energy landscape of a single disc, obtained from micromagnetic simulations using MuMax 330 is represented by the black dots in Fig. 1.Here we plot E/E v − 1 vs. r, where E v is the energy of the pure vortex state (r = 0) and E c is the energy of the collinear state (r = 1).The energy landscape is highly asymmetric, with a maximum at approximately r = 0.8, corresponding to a vortex core positioned inside but close to the edge of the island.As the vortex core reaches the edge and moves outside the disc the energy associated with the magnetic texture rapidly decreases so that E c lies well below the maximum energy.The ratio, E c /E v , determines if an isolated disc carries a moment or not in its lowest energy configuration and this can be varied either side of unity by changing the disc radius 26 .
As shown in the figure, the landscape is qualitatively reproduced by the following phenomenological function where with a = 0.035 and r 0 = 0.85.The agreement between the phenomenological function, S/E v − 1 (solid lines) and the micromagnetic simulations (black dots) is found to be good, as seen in the figure.Also included in the figure are the effects of inter-disc interactions (red dashed lines).Here the lines represent interactions of a disc with spin length r with four fully collinear neighbors (r = 1).Reversal of each neighbor from a parallel to an antiparallel state is therefore characterized by hopping from one line to the next in ascending order.The difference in energy is a measure of the many body interactions that one can expect in an array of discs.Of particular interest here is the case of E c /E v = 1.1, which for the isolated disc indicates preference for a vortex state.When including interactions, however, a collinear configuration is instead favoured with the mesospins lying parallel to each other.This precursor illustrates how inter-disc interactions can influence the collective behaviour of an array promising the emergence of rich many body behaviour.
Magnetostatic inter-disc interactions are typically anisotropic in nature.However, here we assume isotropic nearest neighbour interactions.The choice is doubly motivated: simplicity combined with the possibility of designing hybrid metamaterials with the ascribed properties.With this approach, the Hamiltonian describing interactions between discs, placed on a square lattice taking into account both internal texture and many body interactions can be written J m is the interaction between nearest neighboring islands and θ is the in-plane orientation of the mesospin, 0 ≤ θ < 2π.The first term is similar to the 2DXY model, but includes the varying spin length, 0 ≤ r ≤ 1 and the second term is the parameterized r-dependence defined by Eq. ( 2).The proposed model is similar to the vector Blume-Capel model (VBCM) 27,28,31,32 in which vector spins take discrete lengths (r = {0, 1}).The modeling of the emergent mesospin interactions gives us, in addition the continuous variation of r and the phenomenological energy function S(r).
Anticipating the situation where, below the bulk ordering temperature, the array of interacting discs can be thermally equilibrated, or that the non-equilibrium dynamics can be well represented through an emergent effective temperature, we study the thermal properties of the proposed model.The magnetization, M , is divided up into the mesospin density R and and orientation density Θ, defined The susceptibilities are defined by: , and, . Θ and χ Θ correspond to the magnetization and the magnetic susceptibility of the conventional 2D XY model.It is also convenient to define the parameter ∆E = (E c − E v − 2J m )/E v whose sign designates the preference for broken symmetry or zero spin length in the lowest energy configuration.Using this definition, and with J m = 0.2E v , degeneracy of the internal energy of the discs is obtained when ∆E = −0.4 which marks a point of major importance in this work.

III. METHODS
Arrays of L 2 mesospins with L = 32 on a square lattice with periodic boundaries were simulated using the Metropolis algorithm.Each calculation was based upon 40 000 thermalisation full lattice sweeps prior to 400 000 measurement sweeps, to ensure thermal equilibration and statistically robust results 33,34 .A full lattice sweep entails attempting to update both r and θ once for each lattice site.The required thermalisation time scales were established through monitoring the relaxation time scale for M .E c /E v was varied while always keeping J m = 0.2E v .The initial state at each temperature was set to either a random spin configuration with respect to both θ and r in a "hot start", or an ordered, fully magnetized configuration in a "cold start".

IV. RESULTS
The magnetization, M , obtained from simulation is shown in Fig. 2  start.The results reveal transitions from a high temperature disordered phase to a low temperature quasiordered phase.For negative values of ∆E, a ferromagnetic ground state is energetically favourable.For ∆E large and negative the transitions are smooth, with the finite size magnetization resembling that observed through the KT transition of the 2D XY model, or plane rotator model 35,36 .Increasing ∆E causes a sharpening of the transition, up to an apparent tri-critical point with ∆E ≈ −0.3.Increasing beyond this value, the finite size magnetization undergoes a discontinuous jump, as in a first order transition, see for example ∆E = −0.2.For higher values (∆E ≥ −0.1), the transition into an ordered collinear phase does not occur.
The curves in Fig. 2 are fits to the data sets of the form M = M 0 (T − T c ) β , where T c and β are free parameters and M 0 = 1.Given this phenomenology, one should perhaps consider these curves as guides to the eye, although for ∆E large and negative the results are consistent with the zero parameter fitting procedure outlined in ref. [35]  as well as with many experimental observations 37 .However, rather provocatively, the observed effective exponent β does evolve in a way that is perfectly compatible with observations in a finite system as it crosses over from critical to tri-critical behaviour.The best fit exponents for a few values of ∆E are shown in the upper panel of Fig. 3. Starting from the expected value for the finite 2D XY model for ∆E < −0.9, the fitted β decreases continuously to less than half its initial value, with β ≈ 0.1 close to the apparent tri-critical point.This evolution should be compared with mean field theory where the tri-critical exponent, β tri = 1/4, down from β = 1/2 at the regular critical point and with the 2D BCM where β tri = 1/24 38 is only one third of the 2D Ising critical exponent β = 1/8.We note that the effective tri-critical exponent is quite close to the critical exponent for the Ising model, β = 1/8, although it is difficult to incorporate this observation into a tri-critical scenario for quasiordering of the rotors.
In the lower panel of Fig. 3 we show the evolution of the transition temperature as ∆E increases.For ∆E < −0.9, T C slightly overshoots the extrapolated transition temperature value for the plane rotator model T C ≈ 0.898J m 39 .The overshoot is consistent with the expected logarithmic shift in the effective transition temperature with system size 35 .T C then decreases with increasing ∆E reaching approximately half the plane rotator transition temperature at the tri-critical point before dropping discontinuously to zero for ∆E between −0.3 and −0.2.
The change in the nature of the transition from KT to 1 st order is driven by the change in rotor length, as shown in the upper panel of Fig. 4, where we plot R vs T for different ∆E values.As ∆E passes through −2J m = −0.4 the role played in the free energy by the rotor length changes.For ∆E < −0.4,placing a rotor of maximum length leads to a gain in internal energy for both random and correlated spin configurations.As a consequence, the internal energy favours an ordered state with R = 1, while entropic forces drive R below unity, with maximum entropy for R = 0.5.For ∆E con- siderably greater than -0.4,energy costs are such that R → 0 as T → 0 so that entropy drives the growth in R at all finite temperatures.Between these two limits there is a small window of ∆E for which finite R is energetically favorable if symmetry is broken (or almost broken in the case of a KT transition), otherwise it is entropically driven and energetically unfavourable.This phenomenology is well illustrated in Fig. 4 which shows R to be a monotonically increasing function as temperature is reduced for ∆E < −0.6.For greater values of ∆E, R dips below 0.5 and for ∆E = −0.3, the approximate tricritical value, R clearly decreases as T falls to intermediate values before rebounding to large values through the phase transition.It is this "elastic" resistance to large R values at low temperature that drives the transition first order.At ∆E = −0.2 the transition is clearly 1 st order, while for ∆E = −0.1 no symmetry breaking is observed and R decreases monotonically to zero as T → 0.
This phenomenology in which energy gain at large R depends explicitly on symmetry breaking, is generic to tri-critical systems.Similar behaviour is observed for the VBCM in two-dimensions (see Appendix) and in systems with discrete symmetry 40 .However, there are effects specific to the model studied here which allows for continuous variation of the rotor length.Close examination of Fig. 4 shows that R drops below 0.5 at intermediate temperature, even for ∆E = −0.5 and ∆E = −0.4,while in the equivalent figure for the VBCM (see Appendix), R remains greater than 0.5 for all ∆E ≤ −0.4.This difference arises as the system profits from the entropy associated with a continuous spread of rotor lengths, despite the energy gain from placing rotors of fixed length r = 1.This is clearly a non-universal effect depending on the form of S(r) and could be modified in different nano-engineered arrays.
The KT transition is, from a thermodynamic point of view extremely special as there is no true magnetic symmetry breaking and so no order parameter in the thermodynamic sense.At first sight this might suggest that a 1 st order transition signaled by a discontinuous jump in such a parameter should be excluded.However, the day is saved here by the parameter R, which is a scalar measure of the mesospin density and which is a well defined intensive thermodynamic variable at all temperatures.A thermodynamic signal of the first order transition is therefore a discontinuous jump in R.However, the spin density remains coupled to the magnetization and to the development of quasi-long range orientational order for the mesospins through the tri-critical point and into the first order regime.This is illustrated in the lower panel of Fig. 4, where we show the evolution of Θ with temperature for different values of ∆E.The purely orientation order parameter mimics M , with the pseudo-critical range narrowing as the tri-critical point is approached.
The development of tri-critical coupling between spin density and spin rotation degrees of freedom is illustrated by the three susceptibilities shown in Fig. 5.For ∆E < −0.6, which is deep in the KBT region, χ M and χ Θ show a finite size rounded divergence at the same temperature, which can be fitted to the characteristic exponential form for the KT transition (not shown), while χ R shows a rounded maximum at temperatures that are decoupled from the KT transition and two orders of magnitude smaller than the singular susceptibilities.As ∆E increases into the crossover region towards tri-criticality, a sharp peak in χ R emerges.It rapidly locks onto the divergences in χ M and χ Θ which also sharpen so that, on arriving at the apparent tri-critical point the three susceptibilities show the same sharply singular feature.This can be taken as a signature of the coupling of the internal and external degrees of freedom of the mesospins.
Within the first order regime, the three regions; the high temperature entropic regime, the unfavorable intermediate temperature regime and the broken symmetry ordered phase are illustrated in Fig. 6 in the upper panel.The figure shows snapshots for different temperatures, for ∆E = −0.2.As we are in the regime where, in the absence of interactions, the vortex state is favoured, the low energy magnetic state is generated through many body interactions.In region III both orientational disorder and a wide range of rotor lengths can be observed.In region II while the rotors remain disordered their mean length is clearly reduced, reflecting the energy cost of creating full length rotors while remaining in the disordered phase.In the low temperature phase, I, the symmetry breaking allows for an energy gain on generating extended rotors.In the lower panel we show χ R on a linear scale over the same temperature range.Entry into the intermediate range is marked by a broad maximum at around T = 0.5J m signaling a rapid reduction in mean rotor length.Below the peak, χ R decreases until it hits the first order discontinuity at T ≈ 0.25J m below which the rotor length remains more or less fixed near the maximum value.Also shown is the rotor susceptibility for a non-interacting system (J m = 0).The peak at intermediate temperature is lower and broader when interactions are switched off, illustrating that rotor-rotor interactions offset the energy cost of finite spin length, helping to maintain their presence down to lower temperatures.
The data shown in Figs. 2, 4 and 5 is for hot starts for fixed ∆E.Following this protocol, the symmetry breaking transition disappears between ∆E = −0.2 and −0.1 even though the ground state of an ordered configuration with rotors of unit length remains lower than any disordered state for ∆E < 0. The loss of the transition in this range is a non-equilibrium result characteristic of a first order transition and the presence of metastable states.Making runs from cold starts in the ordered configuration exposes hysteresis in M , due to the loss of ergodicity as shown in Fig. 7.For ∆E = −0.1 the ordered state  The hysteresis can also be observed by ramping ∆E in loops from negative values upwards and back again, while holding the temperature fixed 34 .The full phase diagram in the ∆E, T plane from loops ramping ∆E at fixed temperature is shown in Fig. 8.This phase diagram is similar to that observed for the VBCM 34 .The ordered phase resists finite temperature fluctuations up to ∆E ≈ 0.3, considerably above the equilibrium threshold for stability of the ordered phase.The figure therefore shows a finite region of metastability in which the ordered and disordered phases, labeled I and II coexist for simulations of fixed time scale.The zone of metastability closes at the tri-critical point and ramping ∆E gives an alternative measure of its position, which we estimate to be ∆E = −0.26± 0.01, T tri = 0.375 ± 0.025.This is in good qualitative agreement with our estimate from the evolution of the effective critical exponent.The true phase boundary must run close to the line extrapolating between ∆E = 0 at T = 0 and the tri-critical point, as shown in Fig. 8, although we have not attempted to evaluate it in detail.The crossover between regions II and III in which the rotors are confined to short lengths and in which a full spread of lengths appear is also shown.The position of the broad maximum in χ R which characterizes the crossover is marked for ∆E = {−0.2,−0.1, 0.0}.

V. DISCUSSION
We have shown that engineered two-dimensional arrays of magnetic discs on the mesoscale offer interactions that map convincingly onto a model system showing tri-criticality.In this development, we represent inter-disc interactions and internal spin textures by an effective nearest neighbour coupling between magnetic mesospins of varying length.The energy scale fixing the mesospin length is the magnetic vortex core energy, giving a Blume-Capel type model (BCM) 27,28 with both continuous, in-plane rotor orientations 34,41 and continuous rotor length.The tri-criticality observed numerically is rather special in that it marks the evolution from a Berezinski-Kosterlitz-Thouless phase transition to a first order transition 34,41,42 .
Tri-critical systems are the confluence of three phases whose thermodynamics is governed by three independent thermodynamic variables 43 .As a consequence their critical properties are chracterised by three scaling variables 44 and associated critical exponents 45 .This situation is captured most simply by the Blume-Capel model 27,28 , in which temperature and field conjugate to an order parameter with Z 2 symmetry are joined by an energy scale or chemical potential associated with spin creation and annihilation.
Archtypical examples of tri-crtical systems are the merger of the super-fluid transition of 4 He and the critical point of the demixing transition in 4 He-3 He mixtures 46,47 , or the smectic C * -smectic A transition, which evolves through tri-criticality on mixing two species of liquid crystal 48 .Experimental studies of the tri-criticality are complicated by the difficulty in accessing the three intensive thermodynamic variables.For example for 4 He-3 He mixtures, while exquisite temperature control is possible 49 , the field conjugate to the superconducting order parameter is inaccessible.The mixture can be controlled by varying the 3 He mole fraction which serves as a second order parameter, but the true intensive variable, the chemical potential difference, µ = µ3 He − µ4 He is also inaccessible.In principle, µ could be controlled in ultra-thin helium films 50 but the experimental environment is extremely challenging.In other systems, the situation is even more constrained.Arrays of Josephson junctions can be diluted 51 , allowing the approach to tri-criticality, but the procedure introduces quenched site disorder and the complexity associated with it.First order magnetic transitions occur, for example in FeRh films 52 or in spin ice materials 53 .These transitions could be tuned to tri-criticality 54,55 , but this would require control of both the coupling constant and the chemical potential through applied pressure or site dilution.
Our work opens the door to experiments in which all three intensive thermodynamic variables relative to the tri-critical phase diagram can be controlled through the change of disc spacing, radius and thickness 26 and application of an external field.A change of symmetry for rotor orientations from continuous to discrete and a return to the original BCM is also envisageable through changing the disc shape.
BCM and vector-BCM (VBCM) models have been further extended to Blume-Emergy-Griffiths (BEGM) type models 31,32,56,57 to describe the full 4 He-3 He phase diagram.Here, a bi-quadratic interaction between spins is added capturing isotropic interactions.The extra term introduces the possibility of separating the demixing from the ordering of the internal degree of freedom, allowing for liquid-gas like criticality, 4 He-3 He tri-criticality and a triple point between superfluid, normal 4 He rich and 3 He rich phases in both three 31 and in two dimensions 41,42,58 .The extra bi-quadratic term could also be engineered in the nano-arrays by modification of the disc topology, leading to the development of quadrupole interactions.This would allow for the experimental study of models resembling the BEGM and vector-BEGM in twodimensions.
The magnetization of a system of finite size is a prime experimental indicator of the Kosterliz-Thouless phase transition 35 .It changes in a characteristic manner though the transition, with the emergence of an effective magnetic critical exponent β ≈ 0.23 37 .We have shown here that the effective magnetization exponent crosses over towards tri-criticality in an analagous manner to a thermodynamic exponent, with an effective tri-critical value less than half the critical value, as is the case for the two-dimensional BCM 38,59 .Going beyond this phenomenology would require more extensive numerics and a deeper examination of the theory.This paper provides a platform for this in future work, but more importantly for the present, this straightforward approach provides a platform for experiments on arrays of mesoscale discs in which order parameter crossover, effective or otherwise, should be accessible to measurement.
Further, system size will be an independent control parameter of these metamagnetic systems, providing experimental access to finite size scaling.This is a powerful tool for simulation 59 , but is generally outside the realm of experiment in condensed matter systems.Artificial systems such as the mesospin arrays presented here or cold atom platforms 60 could provide future access to this essential phenomenology.

VI. CONCLUSION
The experimental realisation of emergent tri-criticality in magnetic metamaterials with coupled intra and interisland excitations poses serious experimental challenges, notably the creation of an environment showing equilibrium thermodynamics (real or effective) and controlled departures from it.However it offers a test case that prepares the ground for a vast array of possibilities offered by nano-engineered metamagnets, with both fundamental exploration and technological applications in mind.In particular, the identification of mesospins of variable length in controlled out of equilibrium environments invites applications in adaptive matter 61 through the dynamical modification of the many body energy landscape 62 .The self-modification of the energy discussed here is analogous to the interstitial self-trapping of hydrogen in metals, where the hydrogen interstitial and the local strain field form a dynamic quasi particle 63 .Local energy landscape fluctuations within an array of mesospins could also offer local sensing capabilities and long term memory, a further cornerstones of adaptive or intelligent matter 61 .Clever engineering of the metamaterials may therefore offer pathways towards more advanced materials or even analogue logic mesospin components 64 .

3 FIG. 3 .
FIG. 3. Top panel: The effective critical exponent β as a function of ∆E.The upper dashed line shows β = 0.231, the effective exponent of the 2D XY model.The lower dashed line shows β = 1/8, the value for the 2D Ising model.Bottom panel: The finite size ordering temperature, TC.The dashed line shows TKT = 0.898Jm, the extrapolated value for the 2D XY model.The cross marks the estimated tri-critical temperature.

4 FIG. 4 .
FIG. 4. Top panel: The average mesospin length R vs. temperature for different ∆E .Bottom panel: The orientation density Θ vs. temperature for different ∆E.The open symbols indicate values of ∆E in the 1 st order region of the phase diagram.

FIG. 5 .
FIG. 5. Susceptibility vs. temperature for different ∆E.Top panel: The magnetic susceptibility, χM.Middle panel: The mesospin density susceptibility χR.Bottom panel: The orientation density susceptibility χΘ.The vertical dotted lines show the peak positions overlapping in χM and χΘ, but not in χR for the three lowest value of ∆E.

FIG. 6 .
FIG. 6. Upper Panel: snapshots of mesospin configurations at high (III), intermediate (II) and low temperatures (I).In (III) a full range of r values are observed, while in (II) the range is limited to small values.Transformation from (III) to (II) is a crossover.Evolution from (II) to (I) is via a 1 st order transition for Jm = 0. Lower Panel: χR as a function of temperature for ∆E = −0.2, with and without interactions.

FIG. 7 .FIG. 8 .
FIG. 7. M as a function of temperature for different values of ∆E, starting at random (hot start) or collinear (cold start) spin configurations at each temperature.
survives a cold start up to T = 0.15J m and the metastability survives up to around ∆E = −0.3 in the tri-critical region.