Unconventional magnetism mediated by spin-phonon-photon coupling

Magnetic order typically emerges due to the short-range exchange interaction between the constituent electronic spins. Recent discoveries have found a crucial role for spin-phonon coupling in various phenomena from optical ultrafast magnetization switching to dynamical control of the magnetic state. Here, we demonstrate theoretically the emergence of a biquadratic long-range interaction between spins mediated by their coupling to phonons hybridized with vacuum photons into polaritons. The resulting ordered state enabled by the exchange of virtual polaritons between spins is reminiscent of superconductivity mediated by the exchange of virtual phonons. The biquadratic nature of the spin-spin interaction promotes ordering without favoring ferro- or antiferromagnetism. It further makes the phase transition to magnetic order a first-order transition, unlike in conventional magnets. Consequently, a large magnetization develops abruptly on lowering the temperature which could enable magnetic memories admitting ultralow-power thermally-assisted writing while maintaining a high data stability. The role of photons in the phenomenon further enables an in-situ static control over the magnetism. These unique features make our predicted spin-spin interaction and magnetism highly unconventional paving the way for novel scientific and technological opportunities.


Introduction
Magnets host a broad range of intriguing ground states from spin liquids [1] to topological textures [2], such as skyrmions.These play a central role in the various correlated states of electronic matter and subfields of physics, like spintronics [3] and unconventional superconductivity [4].Magnets have also had a tremendous impact on contemporary computing technology via magnetic random access memories [5] and read heads based on the magnetoresistance effects [6].Different kinds of ordering, such as ferromagnetic and antiferromagnetic, emerge primarily due to the short-range exchange interaction between neighboring spins, which get aligned parallel or antiparallel to each other.
While exchange is typically the strongest interaction in magnets, spinlattice or spin-phonon coupling has been found to underlie the transfer of spin angular momentum between the magnetic and lattice degrees of freedom [7].Phenomena such as ultrafast optical switching of magnetic moments [8][9][10][11][12][13] and long-range transport of spin via phonons [14][15][16][17][18][19][20] rely fundamentally on the spin-phonon interaction.Another of its key consequences is mediating a linear [14][15][16][17][18][19][20][21] or nonlinear [21][22][23] coupling between phonons and magnonsthe spin excitations of ordered magnets.In these considerations, a pre-existing exchange interaction underlies the magnetically ordered ground state while the spin-phonon interaction enables mutual coupling and control between the excitations.Spin-phonon coupling has also been exploited to dynamically control the magnetic state via optically driving certain phonon modes [13,[24][25][26], similar to recent schemes employing light as a strong drive to achieve control over magnetic or even nonmagnetic states of matter, such as a superconductor [27].
At the same time, modifying existing ordered states or phase transitions by tuning the equilibrium electromagnetic environment is currently a highly desired and pursued goal [28,29].Along these lines, the linear-in-spin coupling with light has been predicted to mediate quadratic spin-spin interactions mimicking antiferromagnetic exchange and stabilizing spin liquid states [30,31].A paradigmatic work [32] considered a similar linear-in-pseudospin coupling with phonons to demonstrate a quadratic pseudospin-pseudospin interaction examining its effect on the system's dynamical properties.Such a potential linear-in-spin coupling with phonon displacement is forbidden by time-reversal symmetry [14].
Although spin-phonon coupling has recently been established to play an important part in a large number of nonequilibrium spin phenomena, its potential role in determining the fundamental interaction and ground state of a magnet has not been explored.Here, we theoretically demonstrate the emergence of an unconventional long-range, algebraically-decaying interaction between localized spins due to the exchange of virtual phonons coupled to the vacuum photon modes.The basic phenomenon is similar to how electronelectron attraction emerges in a metal from the exchange of virtual phonons, leading to superconductivity.It is also reminiscent of van der Waals interactions emerging from virtual charge density fluctuations [33].On account of the spin-phonon coupling being quadratic-in-spin due to time-reversal invariance, the emergent spin-spin interaction is found to be biquadratic in the spin components, in contrast with the quadratic nature of the conventional exchange interaction.Thus, depending on the sign of the interaction, it enforces order or disorder in the magnet without explicitly favoring parallel or antiparallel configuration.The possibility to promote disorder can help stabilize spin liquid states [1,34] at higher temperatures.Considering the emergence of ferromagnetism i.e., parallel alignment of all spins, the system is found to manifest a first-order phase transition to a large magnetization just below the critical temperature T c .Thus, this unconventional magnet enables a promising possibility for memories which would admit ultralow-power thermally-assisted writing [35,36] of a bit by raising the temperature slightly above T c and cooling in the presence of a weak applied magnetic field thereby obtaining a large magnetization along a desired direction.Such a process becomes ineffective with conventional magnets that manifest a second-order phase transition because cooling slightly below the critical temperature yields a small magnetization.Since the latter determines the energy barrier between the two equal-energy bit states, the small magnetization results in an unstable and unreliable data storage.

Emergent spin-spin interaction
We consider ferromagnetic nanoparticles, each one bearing a large spin with S ≫ 1 due to its magnetically ordered state and an infrared (IR)-active phonon mode with zero wavenumber and THz range frequency confined to the nanoparticle (see Fig. 1).While the spin-phonon interaction couples these two subsystems, the nanoparticles remain independent of each other at this level.Interaction between the different nanoparticles is provided by the electromagnetic photon modes in the system which are delocalized over the entire space and which couple to the IR-active phonon in each of the nanoparticles.The total Hamiltonian for the system thus becomes where H S describes each of the independent ferromagnetic nanoparticle spins and may include contributions from Zeeman coupling or any local magnetic anisotropies.H P captures the phonon along each Cartesian direction on each of the nanoparticles.H EM accounts for the electromagnetic modes of the environment.
The spin-phonon coupling is obtained as where j runs over all the nanoparticles, S j; k is the kth Cartesian component of the nanoparticle j's spin, β jk is the annihilation operator of nanoparticle j's phonon mode polarized along the kth Cartesian component, and b k parametrizes the spin-phonon interaction strength.The form of this coupling has been derived within a simple model in Supplementary Note 1.It is quadratic in spin components due to time-reversal symmetry while the linear coupling to an IR-active phonon necessitates a noncentrosymmetric ferromagnet [37,38], such as Cu 2 OSeO 3 .Equation ( 2) is obtained by quantizing the classical magnetoelastic coupling Hamiltonian [14,16], appropriately generalized to optical phonons (see Supplementary Note 1), in terms of the phonon ladder operators.The effective spin-phonon coupling b k can be expressed in terms of the material parameters (see Supplementary Note 1).Much of the unconventional nature of the ensuing spin-spin interaction and magnetism is a direct consequence of the spin-phonon coupling [see Eq. ( 2)] being invariant under time-reversal, such that it remains the same when S is replaced by −S.As a result, our main results discussed below are independent of the model details.
The nanoparticle phonon modes couple to the electromagnetic photon modes via the dipole-electric field interaction where n runs over all the electromagnetic modes, α n is the annihilation operator for photon mode n and E n (r) is its electric field.The position vector of nanoparticle j is r j , with d j the dipole moment of its phonon modes, and d jk = d j k.We have further employed the rotating wave and longwavelength (see Supplementary Note 7) approximations in considering the coupling between the photons and the zero-wavenumber IR-active phonons.Diagonalization of the phonon plus electromagnetic Hamiltonian yields polaritonic modes (see Supplementary Note 2) that are delocalized over the whole system.Consequently, the nanoparticle spins interact with the common polaritons in the environment.We obtain the following effective spin Hamiltonian that describes the exchange of virtual bosons by integrating out the polariton modes employing the path integral framework and evaluating the canonical partition function (see Methods, Supplementary Note 2, and Supplementary Note 3) where S 2 ≡ S 2 x , S 2 y , S 2 z denotes a vector made by the spin component squares.Λ is a 3 × 3 tensor describing the coupling in units of energy between the different Cartesian components and depends on the composition and coupling of the polaritonic modes to the spins.Treating the full continuum of electromagnetic modes (see Supplementary Note 4), we obtain where Ω is the phonon frequency (assumed to be the same for all nanoparticles), G is the dyadic Green's function of the electromagnetic field, ϵ 0 is the vacuum electric permittivity, c is the vacuum speed of light, and the spinphonon coupling is assumed to be isotropic, i.e., b k = b for all k.This elegant formula is not restricted to a particular geometry, and thus it allows for studying complex structures while enabling the design of optimized systems.The strength of the spin-spin coupling depends on optical, phononic, and magnetoelastic material parameters as well as on the electrostatic (ω = 0) response of the electromagnetic environment.Notably, it already acts in free space, and does not rely on the presence of a cavity of any kind, nor on achieving strong coupling between light and matter resonances [39], nor on other resonant effects.

Spin-Phonon-Photon Coupling
Fig. 1 Emergent spin-spin interaction and ordering mediated by spin-phononphoton coupling.Top panel: Schematic depiction of the system.Localized spins, illustrated as red arrows, are coupled to local phonons, shown as springs.When we consider the phonons to be coupled with global photons (yellow shading) forming polaritonic modes, an exchange of virtual polaritons between the spins causes an effective spin-spin interaction resulting in an in-plane ordering of the spins.Bottom panel: Non-vanishing components of the spin-spin coupling tensor Λj,j ′ between the central nanoparticle located at r j and the one situated at r j ′ in a square array of spheres in the x-y plane.
Equations ( 4) and ( 5) constitute one of our main results and demonstrate an emergent interaction between the spins.For positive (negative) Λ components, energy is minimized by having a large (zero) value of the corresponding spin component squares.Thus, for positive values of the components of Λ, the emergent interaction encourages ordering of the spins without explicitly preferring ferromagnetic or antiferromagnetic configuration.This symmetry and degeneracy between the two kinds of ordering may be lifted by additional interactions not explicitly considered here.Conversely, a negative value of Λ components promotes disorder.This may reinforce effects such as geometrical frustration in spin liquids [1] leading to their higher stability.
Considering a two-dimensional square array of spherical nanoparticles and employing known material parameters (see Supplementary Note 5), we show the various non-vanishing components of Λ as a function of the distance between the spins in Fig. 1.A strongly anisotropic nature of the emergent interaction can be seen.While the in-plane (x-y) components of Λ are positive on average, thereby supporting an ordered state, the out-of-plane (z) component remains negative for all spins.Thus, the emergent spin-spin interaction encourages the spins to remain in the plane giving rise to two-dimensional magnetism.

Mean-field theory of unconventional ferromagnetism
We now examine the emergence of ferromagnetic order due to the spin-spin interaction derived above.To this end, we continue to assume a twodimensional organization of the nanoparticles in either a square or a hexagonal pattern.Further, we consider spherical or cylindrical nanoparticles.In the former, phonons polarized along any spatial direction interact equally with the electric field.In contrast, for cylindrical particles, the phonon mode polarized along the cylinder axis will couple most strongly to the electric field, so that the interaction can be approximated as being due to a single dipole component.Assuming all spins to be in the same state due to translational invariance and that no external magnetic field is applied, the spin-spin interaction [see Eq. ( 4)] results in the following mean-field Hamiltonian (see Supplementary Note 6) with h j = 4(Λ j; x ⟨S x ⟩ 3 x + Λ j; y ⟨S y ⟩ 3 y + Λ j; z ⟨S z ⟩ 3 z)/S 3 the effective magnetic field.Here, ⟨•⟩ denotes the expectation value and Λ j; x ≡ j ′ Λxx j,j ′ and so on.Thus, only the diagonal components of the Λij tensor contribute to the net magnetic order in the whole ensemble.We obtain the self-consistency equation for determining the x component of each spin (see Supplementary Note 6) and so on for the y and z components where B S (x) is the Brillouin function and β = 1/(k B T ) with k B the Boltzmann constant and T the temperature.Figure 2(a) qualitatively shows the graphical solution of the self-consistency Eq. ( 7) assuming non-vanishing coupling only along the x axis for S ≫ 1, i.e., ⟨S x ⟩ = SB s 4βΛ x ⟨S x ⟩ 3 /S 3 .In contrast with the case of conventional ferromagnetism which admits a unique stable solution to the self-consistency equation, here we find two solutions.This is a direct result of the mean-field h components scaling as the third power of the corresponding spin component expectation value, instead of the first power as is the case for conventional magnetism.
The absolute value of the two solutions as a function of temperature is displayed in Fig 2(b).We find that the spin expectation value develops a finite and large value abruptly as the temperature is lowered, which corresponds to a first-order phase transition.In contrast, conventional magnetism corresponds to a second-order phase transition in which the magnetization increases gradually as the temperature is lowered below the Curie temperature.The solution associated with high expectation value is stable and the other is unstable, because they correspond to a minimum and a maximum in the Helmholtz free energy, respectively.When temperature decreases further, the stable solution replaces the trivial solution, ⟨S x ⟩ = 0, as the energetically favorable state since it becomes the global minimum of the free energy.Furthermore, due to the negative value of Λ j; z , the expectation value of the spin z component vanishes and the spins are oriented in the x-y plane.
Figure 2(c) shows the critical temperature of the ferromagnetic state as a function of the lattice constant for the two kinds of nanoparticle arrays considered, in air, thereby providing guidance for achieving a desired critical temperature by choosing the right arrangement of the nanoparticles.The critical temperature is affected by both the density and the configuration of the lattice.As can also be seen in Figs.2(b), the shape of the nanoparticles modifies the expectation value of the spin x components and the critical temperature, and can thus be used as an additional degree of freedom to obtain the desired result.The reason is the geometrical factor appearing in the self-consistency equation and specifically in h.For the case of cylinders, the effective magnetic field h is proportional to Λ x while for the case of spheres it is proportional to √ 2Λ x due to equal contributions from x and y components, yielding a lower critical temperature.

Static control of magnetism
There has been a tremendous interest in and technological need for controlling magnetic ground states via external knobs.In conventional magnets relying on exchange interaction between the electronic spins, this proves to be a daunting task since it requires controlling the internal states of and interactions between the constituent electrons.Nevertheless, transient control over these has been obtained at ultrafast timescales by creating a strong nonequilibrium in the participating electrons [13,[24][25][26][40][41][42][43][44].Since the unconventional magnetism under consideration relies, in part, on the photon modes, it opens the possibility to control the spin interactions and the consequent magnetic order by modifying these modes, which are easily accessible to the outside world.A recent experiment [45] has already observed an ex-situ change in ferromagnetism via certain resonant effects, distinct from our considerations here.We now demonstrate a strong in-situ tunability of the ferromagnet critical temperature via an external knob, such as a gate voltage, that alters the electromagnetic environment of the system.To this end, we consider that the nanoparticles are dispersed in an electrically insulating medium with static permittivity ϵ m (with free space corresponding to ϵ m = 1).Further, the nanoparticles are now deposited on a substrate and placed at distance d from it [see Fig. 3(a)].The photon modes supported by the system are now modified via the electric field boundary condition imposed by the substrate (see Methods).This directly affects the nature of the polaritons mediating the spin-spin interaction and thus the ferromagnetic state.Figure 3(b) shows the corresponding critical temperature as a function of nanoparticle-substrate separation for insulating and perfectly conducting substrates, where a critical temperature decrease is obtained as the nanoparticles approach the substrate.This is understood by considering that the spin-spin coupling depends on the Green's function at zero frequency, which can be decomposed into two contributions: the free space and the scattered.The latter accounts for the effects due to the environment engineering.
In the case we study here (see Methods), the scattered contribution due to the substrate partially cancels the free-space one, thus decreasing the coupling and consequently the strength of the effective mean field, which determines the critical temperature (see Eq. ( 7)).Interestingly, when the nanoparticles approach the conducting substrate, the critical temperature is reduced dramatically and even vanishes in the limit d → 0. Intuitively, the underlying mechanism can be understood as follows.The image dipoles of the particles in the substrate create a field which tends to cancel out the free-space contribution.Since it is possible to tune, for example, a silicon substrate from being effectively an insulator to a conductor using a gate voltage [46], the critical temperature of the ferromagnetic state may be controlled, in-situ and in equilibrium, via the applied voltage [see Fig. 3(a)].Furthermore, choosing a specific operating temperature, the gate voltage enables a complete turning on and off of the magnetic order.Besides its potential application in memories, this functionality could enable computing architectures based on magnets or spin waves.

Discussion
In evaluating the physical observables here, we have employed materials parameters typical of magnetic materials, such as iron garnets and ferrites, as detailed in the Supplementary Note 5. We note that in order to evaluate the coupling strength of the spin-phonon-photon coupling described here, a candidate material would need to be characterized in terms of both its optical phonons and magnetoelastic properties.This set of parameters is not available for most materials, presumably in part because this specific combination had not been previously identified as playing an important role.Therefore, this work implies that an interesting direction for future research would be to look for materials that maximize the effects described here.For instance, recent experiments demonstrate a significantly enhanced magnon-phonon interaction [47] indicating that our choice of parameters is likely conservative, and better materials would be available.
While we have considered configurations with periodically arranged nanoparticles for concreteness, the qualitative results regarding unconventional ferromagnetism remain the same for arbitrary and disordered two-dimensional arrangements.Hence, it is not crucial for experimental realizations to achieve a perfect pattern.At the same time, our proposal is valid in and offers a new dimension to artificial spin liquids [34] by choosing magnets with optimal optical phonon modes and lithographically fabricating the desired configuration.Furthermore, our proposed unconventional magnets based on spin-phonon coupling may enable a more effective and faster means of implementing cooling in nanostructures based on adiabatic demagnetization, especially considering the unconventional first-order nature of the transition.
In conclusion, we theoretically demonstrate an emergent biquadratic spinspin interaction mediated by spin-phonon-photon coupling.The role of photons in this phenomenon lends a strong anisotropy to the spin-spin interaction making the magnetism two-dimensional as well as enabling static in-situ control over the critical temperature of the magnetically ordered state.Since our investigated unconventional spin interaction emerges from an exchange of virtual bosons, much like superconductivity, our work opens new avenues for exploring similarly unconventional spin interactions and magnetism mediated by different bosonic modes available in solid state systems.

Obtaining polariton modes
The H P , H EM , and H EM−P terms constitute a Hamiltonian which can be represented as where β and α are vectors containing the phonon and photon annihilation operators, Ω and ω are diagonal matrices describing the energy of the phononic and photonic modes, and g describes the coupling between the phonons and photons.The eigenvalues of H PP correspond to the energy of the polaritonic modes, while the eigenvectors contain the coefficients relating the polaritonic operators with the phononic and photonic ones (see Supplementary Note 2).

Integrating out the polaritons
Equation ( 2) can be expressed in terms of the polaritons since the phononic operators are related with the polaritonic ones, and thus the total Hamiltonian reduces to H S , the polaritonic modes and their interaction with spins.The partition function of the system is calculated in the basis of the spin and polaritonic states.The effective spin-spin interaction is obtained by calculating the contribution due to the polaritonic modes employing the path integral framework [48,49] (see Supplementary Note 3).
The coherent state of each polariton is considered.Since the polaritons do not interact with each other, their contribution to the partition function is a product of Gaussian integrals (one for each mode), which have analytical solutions.For each polaritonic mode, the jth spin interacts with the j ′ th spin with strength equal to the product of their couplings with the polariton divided by the polariton's energy.By summing over all polaritonic modes, the spinspin coupling can then be related with the inverse of H PP .For a continuum of modes, the latter can be expressed in terms of the dyadic Green's function of the electromagnetic field (see Supplementary Note 4).

Accounting for a substrate
The dyadic Green's function of a multilayered structure can be efficiently calculated [50].For a two-layered structure, where the layers are characterized by static electric permittivities ϵ m and ϵ s and magnetic permeability equal to unity, the coupling reads Λj,j ′ (r j , with r j , r j ′ corresponding to nanoparticles in the ϵ m layer.
A Supplementary Note 1: Spin-phonon coupling In this note, we derive the spin-phonon coupling term employed in the main text within a simplified and general framework.This also offers guidance with respect to the materials relevant for our proposal.We first motivate the spinphonon coupling term relevant to our analysis on symmetry grounds.Then, we briefly review the well-established theory of spin-phonon or magnetoelastic coupling for ferromagnets and acoustic phonons [14,16,51].This analysis is then generalized to the case of lattice with a basis thereby formulating spinphonon coupling for acoustic and optical phonons in a unified framework.This allows us to establish additional symmetry requirements on the materials that may host our considered spin-phonon coupling term.In addition, the developed framework enables an estimation of the spin-optical phonon coupling based on the more widely available measurements of magnetoelastic constants in magnets [51].

A.1 General form on time-reversal symmetry grounds
The term "spin-phonon coupling" has been employed to discuss a broad range of distinct effects in the literature.Hence, we must first clarify the particular term that we are interested in and justify its relevance.To this end, we begin by considering the lowest order terms in the combined potential energy density for the spin and phonon system.This approach forms the foundation of conventional magnetoelasticity theory [14].The lowest order terms are where S k ≡ S k (r) = S k (r)/S represents the direction cosine of the position dependent spin profile with k representing a Cartesian coordinate, and Q is the generalized displacement coordinate representing the phonon mode in question.In the above list, terms 1, 2, 5, and 6 are forbidden by time-reversal symmetry requiring the Hamiltonian to remain the same under the substitution S k → − S k .The terms 5 through 8 can further be considered higher-order since they involve the gradient of the spin direction cosine in a ferromagnetic ground state and are also disregarded, similar to what is done in the conventional magnetoelasticity theory [14].This leaves terms 3 and 4, with the 4th being higher order than the third.Hence, in this work, we focus on the term 3, ∼ S 2 k Q, the effect of which on the ground state of a spin system has not been considered before.On the other hand, and in contrast with term 3, the term 4 gives rise to a shift in the phonon frequency depending on the magnetic state and has been more widely investigated.
Therefore, the time-reversal symmetry alone warrants that in a strong ferromagnet the leading order effect of spin-phonon coupling is captured by a term ∼ S 2 k Q.However, since we focus on a very specific phonon mode in this work -the zero wavenumber infrared (IR) active optical mode -additional symmetry constraints due to the crystal structure may be present.We examine this in more detail below.

A.2 Magnetoelasticity theory and spin-pair model
Besides relying on time-reversal and crystal symmetry arguments, the conventional magnetoelasticity theory can be constructed using Néel's spin-pair model [51,52] depicted in Fig. S1.It considers a pair of aligned spins making an angle θ with the position vector that separates the spins by a distance R. Within a simplified model, it is assumed that the spin-pair interaction energy W depends on R and θ, thereby admitting a Taylor expansion in terms of even spherical harmonics [51]: where g(R) and l(R) parametrize the separation dependence of the various terms, and we only show the first two even spherical harmonics.The odd ones are forbidden by, again, invariance under time-reversal.Here, the term g(R) independent of θ includes contribution from exchange interaction and does not depend on the spin direction.Thus, it does not cause spin-phonon coupling, but instead a renormalization of the elastic forces and equilibrium lattice configuration.We are not interested in this effect here and thus may approximate The physical origin of this spin-lattice coupling term is material-dependent with contributions from spin-orbit interaction causing single-ion anisotropies, magnetic dipolar fields and so on [14,51,53].
Considering how a strain affects the distance R and direction θ associated with the spin-pair, one may evaluate the resulting change in W . Summing over all the spin-pairs in the unit cell, one obtains an expression of the form: for the magnetoelastic energy in a cubic crystal [14,51].Here, u kk ′ ≡ 1/2(∂R k /∂x k ′ + ∂R k ′ /∂x k ) are strain tensor components, b 1,2 parametrize the spin-phonon coupling, and S k ≡ S k (r) are the components of spatially resolved spin.In ferromagnetic nanoparticles, the so-called macrospin approximation is valid due to the internal exchange being strong.Thus, the spin components do not depend on the position within the nanoparticle and the integral in the equation above simply yields a volume factor thereby enabling an adequate description in terms of the total nanoparticle spin or its direction (see Sec. A.6).The parameters b 1,2 can be evaluated in terms of the function l(R) as shown below, but are typically treated as experimentally determined material parameters.
In considering acoustic phonons, one simply replaces the strain components in the equation (S3) above by the relevant phonon displacements [16].While the above has been evaluated for cubic crystals and acoustic phonons, the main form and features of the relations go well-beyond.They also work for optical phonons [54,55] as well as polycrystalline materials [15].This is because the main physics relates to how the spin-pair energies are affected by rotation and extension/compression, which can be treated on a more generic footing [56].Due to similar generality reasons, the same expression Eq. (S3) has also been employed, successfully and extensively, for multisublattice ferrimagnets [17], such as yttrium iron garnet, despite their highly complex unit cell.

A.3 Diatomic chain model
In this section, we develop a description of spin-phonon coupling within a simple diatomic chain model [57] that explicitly accounts for two sublattices and optical phonon modes.A key goal is to use the spin-pair model in addressing the spin-phonon coupling for acoustic and optical phonon branches within a unified framework.Our simplified model is meant to establish and guide general symmetry and phenomenological considerations for multisublattice magnets, which are not easily treated analytically.
The model [57], depicted in Fig. S2, consists of a one-dimensional lattice with a basis comprising two oppositely charged atoms displaced by a distance d.As a result, the effective spring constants between the neighboring atoms are different depending on the bond distance being d or a − d resulting in the Fig. S2 A diatomic chain model for a two-sublattice magnet hosting acoustic and polar optical phonons [57].The equilibrium locations of the two kinds of atoms (circles) are marked by red crosses and squares.Here, a is the lattice constant and d is distance between the two atoms in the basis.The atomic displacements from their respective equilibrium positions are described by un and vn.potential energy: where K and G are the two effective spring constants.Assuming the same mass M for the two atoms for simplicity, we obtain the dynamical equations of motion: Considering periodic boundary conditions, we assume plane wave solutions of the form and similar for v n .Substituting these in the equations of motion (S5) and (S6), we obtain the eigenspectrum and eigenmodes: The lower (upper) sign in the above expressions corresponds to the acoustic (optical) branch that yields in-phase ũ = ṽ (out-of-phase ũ = −ṽ) displacements of the atoms for q = 0 mode.Furthermore, it can be seen that substituting K = G (corresponding to d = a/2) effectively results in a single phonon branch that is continuously connected at the Brillouin zone boundary q = ±π/a.This is expected from a single-sublattice system the eigenmodes of which have been written down using our two-sublattice notation and thus a smaller Brillouin zone.When considering a spatially homogeneous electric field E ẑ, i.e., within the long-wavelength approximation as applicable for nanoparticles smaller than the relevant wavelengths, the equations of motion (S5) and (S6) are modified as follows where Q a is the charge magnitude for both atoms.This results in the electric field coupling to u n − v n and consequently the q = 0 optical mode that corresponds to ũ = −ṽ [Eq.(S9)].Hence, for coupling to IR light within our polar phonon model, we are primarily interested in the q = 0 mode of the optical branch.
Having established the phonon eigenmodes, we turn to evaluating the spinpair energies assuming, as before, that all the atomic spins are aligned.The spin energy [Eq.(S2)] involving the two spin-pairs is where the subscripts K and G label the two bonds within a unit cell.The spin-pairs' distances will be affected by longitudinal phonon modes along the z direction while transverse phonons influence θ K,G .Here, we focus on the longitudinal phonon modes and thus, treat θ K = θ G = θ as static variables resulting in cos θ becoming the third directional cosine of the spin: S z .The spin energy may thus be expressed as: where we have employed S2), Taylor expanded l K,G functions retaining only the first order terms in the phonon displacements, and employed the concise notation Disregarding the constant offset stemming from the equilibrium configuration, we obtain the spin-phonon coupling energy ∆W S as: where we need to express the atomic displacements in terms of the phonons.Let us first consider the acoustic branch and examine phonons with small, but finite, wavenumber q such that qa ≪ 1. Employing the lower sign in

A.4 Consideration of inversion symmetry
In the previous subsection, we have related the strength of spin-optical phonon coupling to its acoustic counterpart based on the assumption |l ′ G | ≪ |l ′ K | and within the diatomic model introduced in Ref. [57] for studying optical phonons.In reality, we need a weaker condition |l ′ G − l ′ K | ∼ l ′ K for this estimate.Such a condition is ensured by the different bond distances d and a − d in our considered model.We notice that this necessarily breaks the inversion symmetry.While in our model the distances d and a − d need to be different for the existence of optical phonon modes, this is not generically true since we could have obtained the second phonon branch by simple assuming different masses for the two atoms.On the other hand, to justify l ′ G ̸ = l ′ K between the two same spins, we necessarily need the two spin bonds to be different thereby breaking inversion symmetry.Thus, we conclude that our considered coupling between the spin and q = 0 IR active phonon is only available in noncentrosymmetric magnets.
Building on these microscopic insights, we now consider a direct argument to arrive at the same conclusion.Assuming spatial inversion symmetry (due to the crystal structure), only a phonon mode that is odd under spatial inversion may couple to the spatially homogeneous electric field associated with IR light.This is because such a field is odd under spatial inversion while we need the light-phonon coupling term to be invariant under spatial inversion to respect the crystal symmetry.Such a phonon mode cannot directly couple to our assumed aligned spin since the latter is invariant under spatial inversion.Any such coupling between the odd-under-inversion phonon and even-under-inversion spin will violate the inversion symmetry of the crystal.

A.5 Materials
Thus, our considered spin-IR active phonon coupling term is available in noncentrosymmetric magnets.We can thus consider materials with a crystal structure that violates inversion symmetry.At the same time, centrosymmetric materials in the bulk can still become noncentrosymmetric under a straingradient [58] that may be deliberate or accidental.For example, the surface of a ferromagnetic nanoparticle is expected to be highly and nonuniformly strained.Along these lines, buckling of thin layers has also been proposed to convert a centrosymmetric material into a noncentrosymmetric one thereby activating additional couplings [21].
Returning to noncentrosymmetric magnets, a large number of candidates exist and some of them have been studied intensely in recent years due to their hosting skyrmions [37,38].For our considered nanoparticles, the magnetic ground state of each will remain uniform and without any skyrmions, in consistence with our assumptions.One prominent material example is Cu 2 OSeO 3 , in which a phonon magnetochiral effect resulting from Dzyaloshinskii-Moriya interaction has already been observed [59].However, the particular spinphonon coupling parameter relevant to our proposal has not been probed, to the best of our knowledge.
Finding all the optical and spin-phonon parameters relevant to our proposal for a single material has proven to be a daunting task since this specific term of spin-phonon coupling has not been studied much.Also, there is little overlap between such magnetoelastic and IR studies.Hence, in our considerations below, we choose typical parameter values for different materials that have been measured in experiments.We hope that our theoretical proposal will motivate first-principles calculations and/or experiments to characterize the relevant noncentrosymmetric magnets with a focus on this proposal.

A.6 Estimating the spin-phonon coupling strength
As discussed and motivated above, Eq. (S3) for acoustic phonons can be employed for their q = 0 optical counterpart since the latter can be visualized as q = π/a acoustic phonons in an extended Brillouin zone scheme, where a is the lattice constant.Below, we consider the optically-active longitudinal phonon mode polarized along the z direction.The differently polarized phonons can be considered in an analogous manner.Further, we disregard any transverse phonon modes here.
As per Ref. [60], the effective longitudinal strain can be related with the canonical position w(r) by u zz (r) = 2w(r)/(aρ), where ρ is the material density.Quantizing the system as detailed in Ref. [61], we directly obtain β(q)e iq•r + β † (q)e −iq•r , (S23) where Ω(q) is the phonon frequency.The operators β(q) and β † (q) pertain to a continuum description.Within a small nanoparticle, only a uniform mode is optically active.We thus express the quantities in terms of discrete q values and keep only the q = 0 mode, obtaining where V is the volume of the nanoparticle.Therefore, the interaction between the spin and the phonon operators with the help of Eqs.(S3) and (S24), reads where b = b 1 2ℏV /(ρΩa 2 ) quantifies the spin-phonon interaction and has units of energy. where ⊗ denotes the dyadic product between two vectors and N cav indicates the total number of electromagnetic modes.We note that Ω has dimensions 3N p × 3N p , g has dimensions 3N p × M , ω has dimensions M ×M , and Λj,j ′ has dimensions 3×3, where M = N cav +3N p .

C Supplementary Note 3: Effective spin-spin Hamiltonian
An effective spin-spin coupling is derived, based on the spin-phonon-photon interaction, by tracing out the polaritonic degrees of freedom.To do so, we employ the theoretical framework detailed in Refs.[48,49,[62][63][64], which is based on the Euclidean path integral formulation.The thermodynamic functions can be calculated from the canonical partition function Z = Tr[exp(−βH)] with β = 1/(k B T ), H the total Hamiltonian of the system and k B the Boltzmann constant.A convenient basis to calculate the trace of the partition function of the bosonic polaritonic modes is formed by the coherent states |α⟩, which are eigenstates of the annihilation operator a, i.e., a|α⟩ = α|α⟩, and form a complete set, and assuming that there are N spins and M polaritonic modes we have where the sum is taken over all spin states.By tracing out the polaritonic degrees of freedom, we wish to obtain an effective Hamiltonian such that (S38) In order to accomplish this tracing out of the bosonic modes, we employ the result [48,49,63,64] and pause to discuss it.Equation (S39) above is tantamount to replacing the operators a i by their expectation values α i in the coherent state.In this sense, it effectively disregards the quantum commutations between operators and is reminiscent of an analogous replacement procedure within the path integral framework for integrating out excitations [49].Strictly speaking, Eq. (S39) yields the exact result for the free energy per spin in the thermodynamic limit, i.e., when there are a large number of spins (N → ∞).It has been justified semi-rigorously for a finite number M of the bosonic modes [48,62].However, since the number of spins is never infinity, one can consider and employ Eq. (S39) as an approximation that introduces a usually small error depending on N .In this spirit, Eq. (S39) has been successfully and widely employed for capturing the essential physics with a good enough accuracy [48,[62][63][64].We also follow this procedure here.Since the modes are not interacting Z α = exp (−βH S ) i Z αi with The integrals have Gaussian form and we get Then, the total partition function Z α takes the form (S42) In the thermodynamic limit [48,63], we obtain the effective Hamiltonian where ω = diag(ℏω 1 , ℏω 2 , . . ., ℏω M ), and γ is a vector of length M with elements γ ji .We note that the generalization for a spin operator with three components is straightforward, with the coupling Λ becoming a 3 × 3 tensor.

D Supplementary Note 4: Polaritonic coupling for continuum of EM modes
In the derivation above, a discrete number of modes has been considered.Here, we generalize the analysis for a continuum of electromagnetic modes.We follow a macroscopic quantum electrodynamic treatment [65,66], in which the quantized electric field is expressed as with f λ (r ′ , ω) the bosonic annihilation operators of the EM modes, λ an index labeling the electric and magnetic contributions and Gλ (r, r ′ , ω) functions related to the (classical) dyadic Green's function of Maxwell's equations.The electromagnetic interaction of the EM fields with the nanoparticle dipoles in Coulomb gauge within the Power-Zienau-Woolley picture [67] and using the long-wavelength approximation is By formally discretizing the continuum of cavity modes and expressing it through a collective index n = {λ, l, r ′ , ω} (where l denotes Cartesian components), Eq. (S34) can be expressed in terms of the Green's functions where we have used the property λ d 3 r ′ Gλ (r j , r ′ , ω)( Gλ (r j ′ , r ′ , ω)) ⋆T = ℏω 2 ℑ G(r j , r j ′ , ω)/(πϵ 0 c 2 ), where ϵ 0 is the vacuum electric permittivity and c the speed of light in vacuum.We note that we assumed that each particle has E Supplementary Note 5: Parameters For our calculations we consider typical values of magnetic materials.Specifically, the magnetoelastic constant of EuIG b 1 = 10 6 J/m 3 [71], density ρ = 5.4 g/cm 3 [72,73] of YIG, typical lattice constant a = 1.1 nm [72,73] of ferrites and radius r = 100 nm yielding b = 0.094 eV.As discussed above in Sec.A, we are not aware of any material for which all the material parameters required to quantify our proposed effect have been measured.Hence, we have taken typical experimental values available in the literature.Also, we assume a dielectric function given by the Drude-Lorentz model with a single oscillator at angular frequency Ω, oscillator strength f p , loss factor f γ , and background permittivity ϵ bg By comparing the polarizability of a spherical particle in the quasistatic approximations with that of a system with background polarizability and a single dipole transition with dipole moment d, we obtain For ϵ bg = 1, f p = 1, f γ = 0.01, and Ω/(2π) = 100 cm −1 , which correspond to typical values of iron garnets and cuprates [74][75][76], the dipole moment of each nanoparticle equals about 1.7 kD.We note that our material parameters are within realistic range and in order to reach convergence, an 1600 × 1600-nanoparticle array is needed.

F Supplementary Note 6: Mean-field framework
According to the mean field approach, each spin fluctuates around its mean value such that we may expand the spin operator as S = ⟨S⟩ + δs, with δs ≪ ⟨S⟩.With this approach, for the k-th component we have (S i; k ) 2 = (⟨S i; k ⟩ + δs i; k ) 2 ≈ ⟨S i; k ⟩ 2 + 2⟨S i; k ⟩δs i; k keeping up to first-order terms in δs i; k .Thus, a single term of Eq. ( 4) in the main text takes the form (S j; k ) 2 Λkk ′ j,j ′ (S j ′ ; k ′ ) 2 ≈ −3⟨S j; k ⟩ 2 Λkk ′ j,j ′ ⟨S j ′ ; k ′ ⟩ 2 + 2⟨S j; k ⟩ 2 Λkk ′ j,j ′ ⟨S j ′ ; k ′ ⟩S j; k ′ + 2S j; k ⟨S j; k ⟩ Λkk ′ j,j ′ ⟨S j ′ ; k ′ ⟩ 2 , (S53) where we have further employed δs i; k = S i; k −⟨S i; k ⟩.Assuming that the mean value of all spins is the same, Eq. ( 4) of the main text can be cast in the form where ⟨S⟩ 2 ≡(⟨S x ⟩ 2 , ⟨S y ⟩ 2 , ⟨S y ⟩ 2 ), Λ = j,j ′ Λj,j ′ , h j = 2σ Λ j + Λ T j ⟨S⟩/S 3 , Λ j = j ′ Λj,j ′ , and σ = diag(⟨S x ⟩ 2 , ⟨S y ⟩ 2 , ⟨S y ⟩ 2 ).We note that H S has been neglected, which can be realized by a Zeeman term for switched off external magnetic field.It can be shown that for the hexagonal/square arrays in the x-y plane, the mean-field coupling tensor is diagonal and then h j = 4 S 3 Λ j; x ⟨S x ⟩ 3 x + Λ j; y ⟨S y ⟩ 3 y + Λ j; z ⟨S z ⟩ 3 z , (S55) Λ k; j = i Λkk j,j ′ for k = x, y, z.Interestingly, the mean-field Hamiltonian, Eq. (S54), has the same form as in the case of conventional ferromagnetism, in which h j can be viewed as an, effective, internal molecular magnetic field.Following the standard procedure of the mean-field theory [72,77], the self-consistent equation for the expectation value of the k-th component of the spin operator for the j-th nanoparticle is where is the Brillouin function and h = 4 k (Λ k ⟨S k ⟩ 3 ) 2 .We note that the j dependence has been dropped because all spins are identical.
In the limit of S → ∞ the Brillouin function reduces to the Langevin function, L(x) = −1/x + coth(x), and the self-consistent equation reads

Fig. 2
Fig. 2 Temperature dependence of the unconventional ferromagnetic state.(a) Graphical solution of the self-consistency equation for the mean value of spin x component, ⟨Sx⟩ /S = SB S 4βΛx⟨Sx⟩ 3 /S 3 , for various coupling strengths and S ≫ 1.A solution is obtained when the solid line crosses the dashed line y = ⟨Sx⟩.k B and T are the Boltzmann constant and temperature, respectively.(b) Temperature dependence of the absolute mean value of the x component of the spin for an array of cylindrical (black line) and spherical (orange line) nanoparticles, with solid (dotted) lines indicating locally stable (unstable) solutions.The thick solid lines indicate that the solution furthermore corresponds to the global minimum of the free energy.(c) Dependence of the critical temperature Tc on the lattice constant a of the square (solid line) and hexagonal (dashed line) array of spherical (orange line) or cylindrical (black line) nanoparticles.

Fig. 3
Fig. 3 Control of magnetism via modulation of the substrate dielectric constant ϵs.(a) Schematic view of a hexagonal array of cylindrical nanoparticles placed at a distance d from the substrate.The semiconducting substrate can be tuned from being an insulator to a conductor via an applied voltage.The corresponding change in its static dielectric constant alters the polariton modes that mediate the spin-spin interaction.(b) Dependence of the critical temperature Tc of the unconventional ferromagnetic state on d when the substrate is an electrical insulator (ϵs = 2.5; orange line, ϵs = 10; blue line) or a conductor (black line).

Fig. S1 A
Fig.S1A pair of aligned spins make an angle θ with the position vector separating them by a distance R.