Strong-coupling magnetophononics: Self-blocking, phonon-bitriplons, and spin-band engineering

TL;DR

This paper explores how strong spin-phonon coupling in quantum magnets leads to a self-blocking resonance effect, resulting in novel hybrid states and enabling control over spin excitations via ultrafast phonon driving.

Contribution

It introduces the concept of self-blocking in magnetophononics, linking it to collective spin-phonon states and providing a framework for controlling spin spectra in quantum materials.

Findings

Self-blocking limits phonon energy absorption in driven quantum magnets.

Formation of hybrid phonon-triplon states under strong coupling.

Potential for ultrafast control of spin excitations in experiments.

Abstract

Magnetophononics, the modulation of magnetic interactions by driving infrared-active lattice excitations, is emerging as a key mechanism for the ultrafast dynamical control of both semiclassical and quantum spin systems by coherent light. We demonstrate that, in a quantum magnet with strong spin-phonon coupling, resonances between the driven phonon and the spin excitation frequencies exhibit an intrinsic self-blocking effect, whereby only a fraction of the available laser power is absorbed by the phonon. Using the quantum master equations governing the nonequilibrium steady states of the coupled spin-lattice system, we show how self-blocking arises from the self-consistent alteration of the resonance frequencies. We link this to the appearance of mutually repelling collective spin-phonon states, which in the regime of strong hybridization become composites of a phonon and two triplons. We then identify the mechanism and optimal phonon frequencies by which to control a global nonequilibrium renormalization of the lattice-driven spin excitation spectrum and demonstrate that this effect should be observable in ultrafast THz experiments on a number of known quantum magnetic materials.

Figures

Equations

$$H_{\mathrm{s}}=\sum _{i=1}^N\left(J{\vec{S}}_{1,i}\cdot {\vec{S}}_{2,i}+\lambda J{\vec{S}}_{2,i}\cdot {\vec{S}}_{1,i+1}\right),$$

$$H_{\mathrm{p},\mathrm{BZ}}=\sum _{q\in \mathrm{BZ}}{\omega }_{\mathrm{ph}}(q)b_q^{†}{b}_q,$$

$$H_{\mathrm{p}}={\omega }_0{b}_0^{†}{b}_0,$$

$$H_{\mathrm{l}}=\sum _{i=1}^N{E}_0(t)(b_i^{†}+b_i)=N{E}_0(t)\hat{d},$$

$$H_{\mathrm{sp},J}=\sum _{i=1}^{N}g(b_i+b_i^{†})[{\vec{S}}_{1,i}\cdot {\vec{S}}_{2,i}-⟨{\vec{S}}_{1,i}\cdot {\vec{S}}_{2,i}{⟩}_{\mathrm{eq}}],$$

$$H_{\mathrm{sp},J^{\prime }}=\sum _{i=1}^N{g}^{\prime }(b_i+b_i^{†})[{\vec{S}}_{2,i}\cdot {\vec{S}}_{1,i+1}-⟨{\vec{S}}_{2,i}\cdot {\vec{S}}_{1,i+1}{⟩}_{\mathrm{eq}}].$$

$$S_{1(2),i}^{\alpha }=\frac{1}{2}[\pm (s_i^{†}{t}_{\alpha ,i}+t_{\alpha ,i}^{†}{s}_i)-i{ϵ}_{\alpha \beta \gamma }{t}_{\beta ,i}^{†}{t}_{\gamma ,i}],$$

$$S_{1(2),i}^{\alpha }=\pm \frac{1}{2}(t_{\alpha ,i}+t_{\alpha ,i}^{†})+\mathcal{O}(t^{†}t),$$

$$H_{\mathrm{s}}=\sum _{k,\alpha }{\omega }_k{\tilde{t}}_{\alpha ,k}^{†}{\tilde{t}}_{\alpha ,k}$$

$${\omega }_k=J\sqrt{1-\lambda \cos k}.$$

$$y_k$$

$$y_k^{\prime }$$

$$H_{{\rm sp},J}=g\hat{d}\sum_{k,\alpha}\big{[}y_{k}\tilde{t}^{\dagger}_{\alpha,k}\tilde{t}_{\alpha,k}+{\textstyle\frac{1}{2}}y^{\prime}_{k}\big{(}\tilde{t}^{\dagger}_{\alpha,k}\tilde{t}^{\dagger}_{\alpha,-k}+\text{H.c.}\big{)}\big{]}$$

$$H_{{\rm sp},J^{\prime}}=-\frac{g^{\prime}\hat{d}}{2\lambda}\sum_{k,\alpha}y^{\prime}_{k}\big{[}2\tilde{t}^{\dagger}_{\alpha,k}\tilde{t}_{\alpha,k}+\tilde{t}^{\dagger}_{\alpha,k}\tilde{t}^{\dagger}_{\alpha,-k}+\text{H.c.}\big{]}.$$

$$H_{\mathrm{sp},\mathrm{p},J}$$

$$=gN\hat{d}({\mathcal{U}}_J+{\mathcal{V}}_J),$$

$${\mathcal{U}}_J$$

$${\mathcal{V}}_J$$

$$H_{\mathrm{sp},\mathrm{p},J^{\prime }}$$

$$=-(g^{\prime }/\lambda )N\hat{d}({\mathcal{U}}_{J^{\prime }}+{\mathcal{V}}_{J^{\prime }}),$$

$${\mathcal{U}}_{J^{\prime }}=\frac{1}{N}\sum _{k,\alpha }{y}_k^{\prime }⟨{\tilde{t}}_{\alpha ,k}^{†}{\tilde{t}}_{\alpha ,k}⟩;$$

$$\frac{d}{dt}⟨O⟩(t)$$

$$\displaystyle\qquad+\frac{1}{2}\sum_{j}\gamma_{j}\big{\langle}[L_{j}^{\dagger},O(t)]L_{j}+L_{j}^{\dagger}[O(t),L_{j}]\big{\rangle},$$

$${\gamma }_1=\gamma n({\omega }_0),{\gamma }_2=\gamma [1+n({\omega }_0)],$$

$$\frac{d}{dt}q(t)$$

$$\frac{d}{dt}p(t)$$

$$\frac{d}{dt}{n}_{\mathrm{ph}}(t)$$

$$n_{\text{ph}}(t)=\frac{1}{N}⟨b_0^{†}{b}_0⟩$$

$$\tilde{E}(t)=E_0(t)+g[{\mathcal{U}}_J(t)+{\mathcal{V}}_J(t)]$$

$$\tilde{E}(t)=E_0(t)-(g^{\prime }/\lambda )[{\mathcal{U}}_{J^{\prime }}(t)+{\mathcal{V}}_{J^{\prime }}(t)].$$

$${\gamma }_1={\gamma }_{\mathrm{s}}n({\omega }_k),{\gamma }_2={\gamma }_{\mathrm{s}}[1+n({\omega }_k)].$$

$$\frac{d}{dt}{u}_k(t)$$

$$\frac{d}{dt}{v}_k(t)$$

$$\frac{d}{dt}{w}_k(t)$$

$$+2gq(t)y_k^{\prime }[u_k(t)+\frac{3}{2}],$$

Full text

I Introduction

Rapid advances in laser technology Salén et al. (2019) have made it possible not only to probe but also to pump quantum materials in a controlled manner on ultrafast timescales and at all the frequencies relevant to excitations in condensed matter Kampfrath et al. (2013); Zhang et al. (2017); Buzzi et al. (2018). This has led to phenomena ranging from Floquet engineering of electronic band structures Oka and Kitamura (2019) to enhanced superconductivity Cavalleri (2018) and switching of the metal-insulator transition Iwai et al. (2003). A wide range of experimental and theoretical efforts is now under way to extend such ultrafast control to every aspect of strongly correlated materials beyond the charge, including lattice, orbital, spin, nematic, and chiral degrees of freedom de la Torre et al. (2021).

Among these, spin systems offer perhaps the ultimate quantum many-body states due to their intrinsically high entanglement and relatively low energy scales, which lead to rather clean experimental realizations. Ultrafast switching, modulation, transport, and destruction of semiclassical ordered magnetism have been achieved using light of different frequencies Kampfrath et al. (2011); Mikhaylovskiy et al. (2015); Jäckl et al. (2017). However, a direct coupling to a magnetic order parameter is often not appropriate for the dynamical control of quantum magnetic materials, and increasing attention is focused on using the lattice as an intermediary Disa et al. (2020); Afanasiev et al. (2021); Deltenre et al. (2021); Bossini et al. (2021); Disa et al. (2021). While “nonlinear phononics” Först et al. (2011); Subedi et al. (2014); von Hoegen et al. (2018) exploits the anharmonic lattice potential, to date for low-frequency magnetic control Nova et al. (2017), “magnetophononics” Fechner et al. (2018) uses harmonic phonons to effect the highly nonlinear modulation of exchange-type interactions Giorgianni et al. .

The magnetophononic mechanism is ideally suited to the task at hand, namely studying how driving by coherent light can influence the magnetic properties of an insulating low-dimensional quantum spin system. Unless the magnetic interactions are highly anisotropic, the direct coupling of electromagnetic waves to spins is very weak. Using the lattice to mediate this coupling means choosing an infrared (IR) phonon to excite by THz laser radiation so that a coherent lattice oscillation is triggered. Intense irradiation results in a phonon occupation sufficiently high that the (super)exchange couplings between the localized spins undergo a significant alteration Bossini et al. (2021); Giorgianni et al. ; Paris et al. (2021), leading to readily detectable changes in the properties of the magnetic subsystem. While the THz laser can be used to select any IR-active phonon in the spectrum of available lattice excitations, this phonon introduces a frequency that a priori has no direct connection to the intrinsic excitation frequencies of the spin system. Driving a very fast phonon mode ($\omega_{0}$) would put the spin system in the true Floquet regime, where one might seek spin-excitation bands shifted by $\pm n\omega_{0}$ (for small integer $n$). Finding a very slow phonon mode would allow the spin correlations and excitations to be modulated over the course of a single phonon period. Between these two limits, strong excitation of the collective spin modes at their intrinsic frequencies would go beyond these modifications of the existing magnetic states by opening the possibility of creating fundamentally different types of composite collective state, including hybrid spin-spin and spin-phonon composites.

Here we analyze the physics of the magnetophononic mechanism at strong spin-phonon coupling by considering the nonequilibrium steady states (NESS) of a minimal model consisting of an alternating quantum spin chain coupled to a bulk optical phonon mode. When the phonon frequency matches the spectrum of magnetic excitations, we find strong feedback effects between the spin and lattice sectors that produce a number of unconventional phenomena. We demonstrate an intrinsic self-blocking effect, by which a driven phonon in resonance with the peak density of magnetic excitations absorbs little of the driving laser power. We compute the driving-induced mutual renormalization of the lattice and spin excitations, and link the self-blocking to the distinctive phonon-bitriplon hybrid excitations that emerge for phonon frequencies near the spin-band edges. We then demonstrate how all possible phonon frequencies that lie within the spin excitation band can act with varying efficiency to cause a driving-induced global reshaping of the spin spectrum. We discuss the consequences of self-blocking and of this dynamical spectral renormalization for pump-probe experiments on some quantum magnetic materials known to have strong spin-phonon coupling.

The framework for our study is one we have discussed in detail in Ref. Yarmohammadi et al. (2021), where we set out to establish and analyze the equation-of-motion approach to a magnetophononically driven system. In this work we introduced the dimerized chain as a generic model for a gapped quantum spin system, a bulk optical phonon as the most straightforward implementation of the driving mechanism, and the remainder of the phonon spectrum as the dissipative bath. We applied the Lindblad formulation Lindblad (1976) to derive the quantum master equations Breuer and Petruccione (2006) and used these to perform a detailed study of the NESS of the coupled system in the regime where a weak spin-phonon coupling restricted its response largely to linear orders. This analysis revealed the ingredients and parameters of the minimal magnetophononic model, characterized both phonon and spin NESS by their frequency-dependence and wave-vector content, computed the energy flow throughout the driven and dissipative system, related this to the system temperature, and identified the onset of nonlinear feedback effects.

The present study extends the weak-coupling analysis in three key directions. The first is to strong spin-phonon coupling, to identify and investigate the phenomenology of the driven system when the mutual feedback between the spin and lattice sectors becomes strong. Because one fundamental consequence of strong coupling is strong shifts in the characteristic mode energies, the second direction is to perform systematic scans of the driving frequency. Here we comment that the reality of current experiment is somewhat removed from the NESS protocol, using ultrashort and ultra-intense pulses that both contain a broad range of frequencies and produce high-order response processes; as a result, both of these directions constitute essential steps towards an accurate description of experiment. The third direction is that, if the strong drive is used to establish a NESS whose properties differ significantly from those of the equilibrium system, an independent “probe” laser is required to read these properties, preferably by a range of methods sensitive to the magnetic as well as to the lattice sector, and we introduce such a probe.

Our considerations are directly relevant to at least two quantum magnetic materials, CuGeO3 and (VO)2P2O7, which were found over 20 years ago to be quasi-one-dimensional alternating spin chains with extremely strong spin-phonon coupling. In CuGeO3, this coupling is strong enough to drive a spin-Peierls transition into the dimerized phase below $T_{\rm sp}=14$ K Hase et al. (1993), whereas (VO)2P2O7 is intrinsically dimerized and was found by raising the temperature to show strong renormalization of the phonons by the spin sector Grove et al. (2000). While the spectra of spin and lattice excitations were studied in detail in both materials Popović et al. (1995); Regnault et al. (1996); Uhrig (1997); Braden et al. (1998); Werner et al. (1999); Braden et al. (2002); Garrett et al. (1997); Grove et al. (2000); Uhrig and Normand (2001), they have yet to be considered from the standpoint of matching drivable (IR) phonons to particular frequency regimes within their spin spectra.

The structure of this article is as follows. In Sec. II we introduce the two phonon-coupled alternating spin-chain models we study and the equation-of-motion method by which we compute their driven dynamics. In Sec. III we analyze the phenomenon of self-blocking. The properties of the phonon-bitriplon hybrid states are presented in Sec. IV. Section V studies the dynamical modifications of the spin excitation band and thus demonstrates the potential for spin-band engineering in quantum magnetic materials. The relevance of these findings to two very strongly spin-phonon-coupled quantum spin-chain materials, CuGeO3 and (VO)2P2O7, is discussed in Sec. VI and Sec. VII contains a brief summary and conclusion.

II Models and Methods

Following the logic of Ref. Yarmohammadi et al. (2021), we consider a model for a magnetophononically driven quantum spin system that is minimally complex but nevertheless contains all of the components essential for capturing the physics of real materials. We do not focus on long-ranged magnetic order, because this is not generally stable in a truly low-dimensional magnet. Thus we consider an alternating spin chain, i.e. a one-dimensional quantum magnet with gapped spin excitations. Without loss of generality, and in order to make progress with a straightforward perturbative treatment of the spin system, we use a chain with a substantial dimerization, setting the interaction on the weaker bonds ($J^{\prime}$) to half the size of the stronger bonds ($J$). The IR phonon acting as the intermediary between the driving laser and the spin system can couple to the magnetic interactions in a wide variety of different ways. Here we restrict our considerations to the leading (linear) term in a Taylor expansion and analyze the two distinct coupling geometries shown in Fig. 1, where (i) the phonon modulates the strong (intradimer) bond [Fig. 1(a)], to which we refer henceforth as the “$J$-model,” and (ii) the phonon modulates the weak (interdimer) bond [Fig. 1(b)], which we call the “$J^{\prime}$-model” Yarmohammadi (2022). We will show that the two models yield very similar results for certain magnetophononic phenomena but are quite different for other phenomena.

II.1 Hamiltonian dynamics

The Hamiltonian of the spin system takes the form

[TABLE]

where $\lambda=J^{\prime}/J$, $i$ labels the dimer, 1 and 2 denote the two spins within each dimer, and periodic boundary conditions are assumed. The full phonon Hamiltonian is

[TABLE]

where we have omitted an additional quantum number for the different phonon branches. The acoustic phonons make the largest contributions to the damping, both of optical phonons and of spin excitations, and as a result their effects are included in the phenomenological damping coefficients to be introduced below. For the purposes of magnetophononic driving, we use the frequency of the incoming THz laser radiation to select a single, IR-active optical phonon, which without loss of generality can be dispersionless. The only relevant phonon momentum is $q=0$, because of the dispersion relation (extremely high speed) of the incident light, and hence any phonon dispersion plays no role. The only phonon term in the Hamiltonian of the driven system is then

[TABLE]

where we use $b_{0}$ and $\omega_{0}$ as shorthand for $b_{q=0}$ and $\omega_{\text{ph}}(q=0)$. A further Hamiltonian term is required to describe the driving of this phonon by the electric field, $E_{0}(t)$, of the laser,

[TABLE]

where $\hat{d}=(b^{\dagger}_{0}+b_{0})/\sqrt{N}$ is the operator specifying the local atomic displacement. The linear dependence on $N$ indicates clearly that any finite driving induces a finite value of the phonon displacement observable, $q(t)=\langle\hat{d}\rangle$, and hence a macroscopic occupation ($n_{0}\propto N$) of the $q=0$ boson.

To complete the magnetophononic Hamiltonian we specify the two types of spin-phonon coupling shown in Fig. 1. For the $J$-model,

[TABLE]

where the second term denotes the equilibrium value of the spin interaction on the strong bonds of the chain, and its presence ensures that the dimerization, $\lambda$, does not change in absence of the driving term. For the $J^{\prime}$-model,

[TABLE]

The spin-phonon coupling coefficients have units of energy, and for convenience we will normalize $g$ to $J$ and $g^{\prime}$ to $J^{\prime}$. While the two coupling types are dichotomous, in almost any real material one may expect the atomic displacements associated with any phonon mode to include components that alter all of the magnetic interaction terms in the system.

We proceed by diagonalizing the spin system, for which we introduce bond operators expressing the creation and annihilation of singlet and triplet states on each dimer, $i$. In the relevant limit, where small numbers of triplets form the elementary excitations (henceforth “triplons”) above a sea of singlets, the exact identity Sachdev and Bhatt (1990); Normand and Rüegg (2011)

[TABLE]

can be reduced to

[TABLE]

where $\alpha\in\{x,y,z\}$ denotes both the spin component and the triplon flavor. The full expression [Eq. (7)] shows explicily how triplon creation is accompanied by singlet annihilation, and vice versa, ensuring the hard-core property of the triplons. It is also the basis of a systematic perturbative approach that could be used to perform accurate calculations for alternating chains with much weaker dimerization ($\lambda\rightarrow 1$) Knetter and Uhrig (2000); Schmidt and Uhrig (2003). However, for the moderate values of $\lambda$ that we consider ($\lambda\lesssim 1/2$), a reliable description of the elementary magnetic excitations is obtained by using only the first term of Eq. (8) (i.e. restricting the Hamiltonian to bilinear terms) and neglecting the hard-core property of the triplons Yarmohammadi et al. (2021). A Fourier transformation and a Bogoliubov transformation, the latter using the basis of diagonal triplons $\tilde{t}_{\alpha,k}$ Yarmohammadi et al. (2021), brings the spin Hamiltonian to the form

[TABLE]

with dispersion relation

[TABLE]

To apply these transformations to the Hamiltonian describing the spin-phonon coupling, we introduce the wave-vector-dependent coefficients

[TABLE]

With these we express the $J$-model in the form

[TABLE]

and for the $J^{\prime}$-model we obtain

[TABLE]

These two equations allow us to observe the leading differences and similarities of the two models. The most striking difference is that the prefactor $gy_{k}$ of $\tilde{t}^{\dagger}_{\alpha,k}\tilde{t}_{\alpha,k}$ in the $J$-model is changed to $-g^{\prime}y_{k}^{\prime}/\lambda$ in the $J^{\prime}$-model, which amounts to a sizable decrease because $|y^{\prime}_{k}|\ll y_{k}$ for most $k$ values. An intriguing similarity arises in the prefactor of the pair-creation and -annihilation terms, where $gy^{\prime}_{k}$ is changed to $-g^{\prime}y^{\prime}_{k}/\lambda$. The sign of the prefactor matters little, because it can be changed by the unitary transformation $\tilde{t}_{\alpha,k}\to i\tilde{t}_{\alpha,k}$, and thus we anticipate that similar results are to be expected if one compares $J$- and $J^{\prime}$-models with the property $g/J=g^{\prime}/J^{\prime}$ (i.e. $g=g^{\prime}/\lambda$). We will illustrate this situation in Secs. III and IV.

The spin-phonon coupling contains trilinear bosonic terms incorporating two triplon operators and the displacement operator, $\hat{d}$, of the driving IR phonon. We treat these trilinear terms by a dynamical mean-field approach. For the “spin part” of this term, the mean-field procedure consists of replacing $\hat{d}$ by its expectation value, $\langle\hat{d}\rangle=q(t)$, and keeping the action of the triplon operators. For the “phonon part” of this term, we replace the spin part by its expectation value to obtain for the $J$-model

[TABLE]

where we used

[TABLE]

For the $J^{\prime}$-model we obtain the analogous form

[TABLE]

where $\mathcal{V}_{J^{\prime}}=\mathcal{V}_{J}$ and

[TABLE]

we draw attention to the replacement $y_{k}\to y^{\prime}_{k}$ relative to the $J$-model. We stress in addition that the phonon oscillation and the expectation value of the total spin system are both extensive quantities, as a result of which their relative quantum fluctuations tend to zero in the thermodynamic limit ($N\to\infty$), which provides an excellent justification for the mean-field decoupling we employ.

II.2 Quantum master equations

If a real quantum mechanical system is driven continuously, the absorbed energy will cause heating, which will on some timescale push the system beyond its quantum regime, if not also to very high temperatures (with modern laser intensities one may even surpass the melting point). A systematic treatment of the energy flow requires the considerations of an open quantum system, where relaxation and dissipation channels are included in addition to the external drive. For the spin chain pumped by an IR optical phonon (Fig. 1), we showed in Ref. Yarmohammadi et al. (2021) that the dissipation should be included on two levels, specifically the damping of the driven IR phonon and a direct damping of the triplon modes. Both are assumed to have their microscopic origin in the ensemble of phonon modes, particularly the acoustic ones, and both are treated by means of the adjoint Lindblad master equation Breuer and Petruccione (2006)

[TABLE]

where $H$ is the Hamiltonian of the isolated system (the spin sector and the driven phonon), $O(t)$ is an operator in the Heisenberg picture, $\{L_{j}\}$ are Lindblad operators (operators of the isolated system that link it to its environment, the “bath”), and $\{\gamma_{j}\}$ are the corresponding damping coefficients (decay rates). The Lindblad framework requires that the coefficients $\{\gamma_{j}\}$ be relatively weak, but places no constraint on the terms within the isolated system, meaning that it can be applied for all values of $g$ and $g^{\prime}$ [Eqs. (5) and (6)].

To describe the dissipation of the phonon, we follow the conventional choice Breuer and Petruccione (2006) $L_{1}=b_{0}^{\dagger}$, $L_{2}=b_{0}$, and parameterize the decay rates using

[TABLE]

where $n(\varpi)$ is the bosonic occupation number at energy $\hbar\varpi$. The dynamics of the phonon are then that of a driven and damped harmonic oscillator,

[TABLE]

where $p(t)=i\langle b^{\dagger}_{0}-b_{0}\rangle/\sqrt{N}$ is the momentum conjugate to $q(t)$,

[TABLE]

is the number of phonons per dimer, and $\widetilde{E}(t)$, which denotes the effective electric field acting on the phonon in the presence of its coupling to the spin system, is defined for the $J$-model by

[TABLE]

and for the $J^{\prime}$-model by

[TABLE]

To describe the dissipation of the triplons, we proceed in a similar way by choosing the Lindblad operators for each triplon $k$-mode to be $L_{1}=\tilde{t}_{\alpha,k}^{\dagger}$ and $L_{2}=\tilde{t}_{\alpha,k}$, with the corresponding decay rates given by

[TABLE]

Taking $\gamma_{1}$ and $\gamma_{2}$ to be independent of $\alpha$ is a consequence of the isotropy of the spin system, but taking them to be independent of the momentum, $\hbar k$, is a simplifying approximation that we make to avoid overburdening the model with a multitude of parameters. This assumption can be justified by the fact that we consider a one-dimensional spin system whose energy is dissipated into a bath of three-dimensional phonons. In this geometry, for any given wave vector, $k$, in the chain direction there remain two perpendicular directions over which one has to sum in order to capture the full phonon bath, whence one does not expect a strong dependence on $k$ within this continuum of dissipation channels.

Here we remind the reader that the Lindblad operators $\tilde{t}_{\alpha,k}^{\dagger}$ and $\tilde{t}_{\alpha,k}$ correspond to the creation and annihilation of a quasiparticle with spin. As a consequence, these assumed bath processes do not conserve the total spin and thus are not in fact consistent with the form of the spin-phonon coupling assumed in Eqs. (5) and (6), where the isotropy of the spin part means that a phononic damping process cannot change the spin quantum number. As explained in detail in Ref. Yarmohammadi et al. (2021), this inconsistency would be repaired for a spin-isotropic material by assuming a more complex and spin-conserving form for the bath operators (for example $C_{kq}=\tilde{t}_{\alpha,k}^{\dagger}\tilde{t}_{\alpha,q}$) and for a material with anisotropic spin interactions, usually a consequence of finite spin-orbit couplings, by adopting a more complex form for the spin-phonon coupling. However, to make progress in elucidating the phenomenology of strong-coupling magnetophononics in the most transparent way possible, we proceed with the present minimalist formulation of the problem captured by the spin-chain model of Sec. IIA.

The equations of motion for the spin sector of the $J$-model that result from the Lindblad master equation [Eq. (18)] take the form

[TABLE]

where

[TABLE]

Analogously, for the $J^{\prime}$-model we obtain

[TABLE]

The full system of equations of motion then contains three equations for the driving phonon and $3N$ for the triplons (for a system consisting of $N$ dimers). Because every triplon $k$-mode is coupled to the phonon variables, and the latter to sums over all the triplons, the system cannot be split into $N$ separate sets of differential equations. The inversion symmetry of the chain ensures that $y_{k}=y_{-k}$, $y_{k}^{\prime}=y_{-k}^{\prime}$, $u_{k}=u_{-k}$, $z_{k}=z_{-k}$, and $\omega_{k}=\omega_{-k}$, which reduces the $3N$ triplon equations to $3(N+1)/2$ for odd $N$. We solve these coupled differential equations numerically with an adaptive Runge-Kutta solver, which allows long times, $t$, to be accessed reliably. The chain lengths we consider vary between $N=1001$ and $4001$ in order to ensure that finite-size effects are well controlled in all regimes.

In the analyses to follow, we will characterize the system by introducing a number of measures. For the spin system it is convenient to use the number of excited triplons per dimer, which is given by

[TABLE]

If a system is driven at frequency $\omega$, for a time sufficiently long that it has reached a NESS, no other frequency will appear in the expectation values of the observables, regardless of the available fundamental frequencies in the system (notably $\omega_{0}$ when this differs from $\omega$). Only higher harmonics at integer multiples of $\omega$ will appear, and these are expected in any coupled system Yarmohammadi et al. (2021). In order to focus on the important average values of the time-dependent quantities, for any expectation value $X(t)$ we define

[TABLE]

which represents the average of $X$ over one period, $T=2\pi/\omega$. Below we will consider quantities including $n_{\rm ph0}$, $n_{\rm x0}$, and $u_{k,0}$. Only if one considers the transients appearing directly after switching on the drive do the average values, $X_{0}$, acquire a time-dependence, $X_{0}(t)$ Yarmohammadi et al. (2021).

The focus of our calculations is on the NESS of the driven, dissipative system. We consider a representative alternating chain with moderate dimerization, $\lambda=0.5$ in Eq. (10), which places the edges of the two-triplon excitation band at $2\omega_{\rm min}=1.414J$ and $2\omega_{\rm max}=2.449J$ Yarmohammadi et al. (2021). When the spin-phonon coupling is weak, NESS formation requires a timescale of approximately $5/\gamma_{\rm s}$, i.e. five time constants of the spin system Yarmohammadi et al. (2021). However, at large $g$ values one may anticipate strong feedback processes between the phonon and spin sectors, making it necessary to examine the situation in detail, and potentially to wait for significantly longer times to ensure NESS formation. An example of complementary phonon and spin NESS, each characterized by their number operators, is shown in the time domain for a nonresonant phonon frequency in Figs. 2(a) and 2(b) at weak, strong, and very strong values of $g$. Here we have set the driving frequency to the phonon frequency ($\omega=\omega_{0}$), but neither lies within the two-triplon band. It is clear that both time traces contain increasingly complex combinations of harmonics as $g$ is increased, with second-harmonic contributions dominating at the chosen frequency, and also that the amplitude of the oscillatory part of the phonon occupation behaves non-monotonically as a function of $g$, becoming suppressed at very strong $g$. The oscillatory part of the triplon occupation rises very strongly with $g$ on exiting the weak-coupling regime, but is also suppressed at very strong coupling. The corresponding static parts of both occupations are considered in Sec. III below. Concerning the timescale for NESS formation, Fig. 2(c) shows how the initial rise of $n_{\rm ph}$ is truncated by the rise of $n_{\rm x}$ [Fig. 2(d)]. In this nonresonant regime, at $g=0.3J$ there remains a significant time lag between the driving phonon and the following triplon occupations, where the latter limit the former and convergence requires one slow oscillation cycle, whose length is determined by the feedback process. At $g=0.5J$, the lag in response is much shorter and several slow oscillation cycles are required.

We close our technical presentation with a number of comments. The equations of motion are valid at all times from the onset of driving ($t=0$) to infinity and for all applied electric fields, as well as for all phonon occupations up to the Lindemann melting criterion ($n_{\rm ph}\approx 3$). With the present simplified treatment of the spin sector, they are valid up to a triplon occupation of order $n_{\rm x}\approx 0.2$, beyond which a more sophisticated numerical treatment should be used to account for the hard-core occupation constraint. Because the equations of motion are based on a mean-field decoupling of the spin and lattice sectors, the treatment we present becomes more approximate at low phonon frequencies, specifically those below $\omega_{0}=0.2$–$0.3J$ Yarmohammadi et al. (2021). Nevertheless, one may verify by considering the energy flow through the strongly spin-phonon-coupled system that the mean-field approximation remains very accurate at all phonon frequencies close to resonance with the spin system.

Finally, one may question the stability of the alternating chain in the presence of phononic driving, particularly when this is very strong or very slow. In fact the sharp fall in the driven phonon occupation at very small $\omega_{0}$ in Figs. 3(a) and 3(d) below is related to a ground-state instability of the chain, where a stimulated distortion can occur (the average phonon displacement, $q_{0}$, becomes finite) in the presence of sufficiently slow phonons. One may show that the stability criterion takes the form $\omega_{0}^{c}>F(\lambda)g^{2}\lambda^{2}J$, and that for a $J$-model with $\lambda=0.5$ and $g/J=0.5$ this critical value is $\omega_{0}^{c}\simeq 0.07J$, while in a $J^{\prime}$-model with $\lambda=0.5$ and $g^{\prime}/J^{\prime}=0.5$ it is $\omega_{0}^{c}\simeq 0.14J$.

III Self-blocking

III.1 NESS protocol

We consider first the NESS established by steady laser driving at the frequency of the target IR phonon, i.e. $\omega=\omega_{0}$. In Fig. 3(a) we show $n_{\rm ph0}$ at this frequency choice, which we denote ${\overline{n}}_{\rm ph0}$, as $\omega_{0}$ is varied across the full frequency range for the $J$-model. We use a laser electric-field strength ($E_{0}=0.2\gamma$, expressed in energy units with $\hbar=1$) and phonon damping ($\gamma=0.02\omega_{0}$) that we maintain constant for the remainder of the analysis, and refer to these henceforth as standard driving and damping conditions. At small $g$, ${\overline{n}}_{\rm ph0}$ is effectively constant for all $\omega_{0}$, but as $g$ is increased we observe an increasing suppression of ${\overline{n}}_{\rm ph0}$ that sets in precisely where the density of two-triplon excitations is highest, i.e. at $2\omega_{\rm min}=1.414J$ and $2\omega_{\rm max}=2.449J$. This resonant effect becomes gigantic at strong $g$, suppressing the phonon occupation by nearly three orders of magnitude at $2\omega_{\rm min}$.

We have named this effect “self-blocking” because the magnetic system acts to block its own energy uptake by suppressing the driven phonon. This behavior is surprising if one expects stronger energy absorption when more spin excitations coincide with the driving laser frequency. Its explanation lies in the fact Yarmohammadi et al. (2021) that in magnetophononic driving the spin system is not coupled to the light, but only to the driven phonon. In the NESS protocol, the light frequency is fixed but the effective phonon frequency is altered by its hybridization with the spin system, whose dependence on $g$ we analyze in detail below, and thus the laser driving becomes increasingly non-resonant. Analytically, the prefactor of the phonon momentum, $p(t)$, in the master equation for $n_{\rm ph}(t)$ [Eq. (20c)] is not the driving electric field, $E_{0}(t)$, but the quantity ${\widetilde{E}}(t)$ specified in Eqs. (22) and (23). This effective feedback from the spin system is both strongly nonlinear in $g$ and strongly negative, acting to cancel $E_{0}(t)$ almost completely when $\omega_{0}$ is at resonance with the band edges [Fig. 3(a)]. Despite the approximate symmetry of the two-triplon band, self-blocking is weaker by a factor of 10 at $2\omega_{\rm max}$ due to matrix elements within the feedback process.

Turning to the response of the spin system, Fig. 3(b) shows the corresponding average triplon occupancy, ${\overline{n}}_{\rm x0}$. The most striking feature is the strong rounding of the in-band response as $g$ is increased. The band-edge peaks are entirely blunted by the strong suppression of ${\overline{n}}_{\rm ph0}$ [Fig. 3(a)]. We stress that the effective limiting value ${\overline{n}}_{\rm x0}\approx 0.1$ visible in Fig. 3(b) is purely a consequence of the giant self-blocking, and is not connected with the hard-core nature of the triplon excitations, which has not been included in the formalism of Sec. II. This rounding indicates an increasing localization of the spin response, by which the band character of the triplons becomes less relevant under strong driving by the entirely local phonon. Figure 3(b) also displays a somewhat counterintuitive non-monotonic behavior for frequencies close to the band edges, where increasing $g$ leads to a lower number of excited triplons due to the lower number of driving phonons caused by the self-blocking. Normalizing ${\overline{n}}_{\rm x0}$ to ${\overline{n}}_{\rm ph0}$, as shown in Fig. 3(c), reveals a set of near-identical response curves sorted in ascending order of $g$, and hence that larger values of the spin-phonon coupling do indeed lead to larger numbers of excited triplons per excited phonon.

Still one might suspect that self-blocking is a special feature of a strongly dimerized $J$-model, in the sense that a $k=0$ phonon strongly coupled to the intradimer bonds could push the system into a completely local limit of isolated spin dimers. However, to show that self-blocking is a general feature of a magnetophononic model, that occurs also when the geometry does not allow the phonon to cut the spin system into local subsystems, we consider the properties of the $J^{\prime}$-model, illustrated in Fig. 3(d). It is clear on the qualitative level that the self-blocking phenomenon is identical, with strong suppression of the absorbed laser energy, and hence of the driven phonon occupation, setting in at the band edges and rising dramatically with $g^{\prime}$. On the quantitative level, if one compares models with the same values of $g/J$ and $g^{\prime}/J^{\prime}$ the results for ${\overline{n}}_{\rm ph0}$ and ${\overline{n}}_{\rm x0}$ [Fig. 3(e)], and hence for the normalized quantity ${\overline{n}}_{\rm x0}/{\overline{n}}_{\rm ph0}$ [Fig. 3(f)], are similar to the degree that they cannot be distinguished on logarithmic intensity scales. This is a consequence of the close similarities between the spin-phonon coupling terms, and hence between the equations of motion, discussed in Sec. II. We comment also that in both models the dominant self-blocking effects are concentrated around the edges of the two-triplon spectrum of the isolated spin system, indicating that, despite any tendency towards localization favored by the strong phonon coupling, the band character of the spin system is largely preserved.

Away from the two-triplon band, in Figs. 3(a) for the $J$-model and Fig. 3(d) for the $J^{\prime}$-model we also observe a significant suppression of phonon energy entering the system at any frequency $\omega_{0}<2\omega_{\rm min}$. This nonresonant self-blocking is also nonlinear in $g$, exceeding one order of magnitude at $g=0.5J$ and $g^{\prime}=0.5J^{\prime}$. Its appearance only in the low-$\omega_{0}$ regime, but not at $\omega_{0}>2\omega_{\rm max}$, points to an origin in multiple harmonic processes ($2\omega_{\rm min}\leq n\omega_{0}\leq 2\omega_{\rm max)}$) Yarmohammadi et al. (2021). Although only the two-phonon harmonic ($n=2$) at $\omega_{\rm min}$ is visible directly, stronger $g$ distributes the response of the system to a given $n\omega_{0}$ across a broader range of frequencies. By contrast, a driving phonon at the band center ($\omega_{0}=2J$) has vanishing matrix elements ($y^{\prime}_{\pi/2}=0$) with the resonant spin modes, and hence ${\overline{n}}_{\rm ph0}$ recovers almost to its $g=0$ or $g^{\prime}=0$ values for all $g$ or $g^{\prime}$.

To understand the context of these results, we stress again that our results are obtained for an idealized NESS experiment, where, as stated in Sec. II, in the long-time limit there are no frequencies in the system other than $\omega$, which has been selected equal to $\omega_{0}$, and multiples thereof. In this sense the panels of Fig. 3 must be interpreted “vertically,” because there is no possibility of spectral-weight transfer between different frequencies. Although the characteristic resonant frequencies of the spin-phonon-coupled system are shifted as $g$ is increased, and thus the system is simply off-resonance for energy uptake, in the $\omega=\omega_{0}$ NESS protocol there is no other option. Hence the self-blocking observed in Fig. 3 can be sought in experiments constructed to achieve a NESS, and we have found in this limit that it is a truly giant phenomenon.

III.2 Pulsed protocol

For a broader view, however, one does wish to understand the full response of the driven system beyond the NESS protocol. It is well known even in the absence of driving that strong spin-phonon coupling leads to hybridization and anticrossing of the spin and phonon excitations, such that the bare spin and phonon dispersion relations are no longer the characteristic resonant frequencies. Before presenting our results for systems with fixed phonon frequencies, $\omega_{0}$, driven at a range of different frequencies, $\omega$, we note that a conventional ultrafast experiment already introduces a spectrum of driving frequencies within the envelope of each ultrashort pulse, and hence allows automatically for the “horizontal” transfer of spectral weight (i.e. between frequencies). For this reason we use the terminology “pulsed protocol,” although in the remainder of the present work we will compare the NESS obtained in systems driven continuously at variable $\omega$ (with constant laser intensity), leaving the accurate modelling of pulsed driving to a future study.

We focus primarily on the $J$-model and, because self-blocking is strongest at the band edges, in Fig. 4 we consider phonon frequencies, marked by the vertical dashed lines, just below and above each of the band edges. Figure 4(a) makes clear that the phononic response of a system driven at a frequency just below the lower band edge (here $\omega_{0}=1.35J$) has a conventional Lorentzian resonance centered at the bare phonon frequency for small $g$, but is weakened and pushed away from the band edge at stronger $g$. Figure 4(b) shows the analogous result when $\omega_{0}$ lies just inside the spin band, where the phonon peak is damped very strongly with increasing $g$, and also moves away from the band edge by the same level-repulsion effect. Here it is accompanied by the development of a second feature, appearing at $2\omega_{\rm min}$ at $g=0.1J$, which is repelled below the band edge as $g$ increases. The accompanying spin response [Figs. 4(c) and 4(d)], which we analyze in Sec. IV below, shows that the appearance of two mutually repelling peaks is generic, and that the characteristic excitation frequencies are shifted quite significantly away from $\omega_{0}$ and $2\omega_{\rm min}$ at large $g$. The situation for phonon frequencies ($\omega_{0}=2.4J$ and $2.5J$) bracketing the upper band edge is exactly analogous, although with a slightly weaker mutual repulsion [Figs. 4(f) to 4(i)].

In the context of self-blocking, these results show how the phononic response is shifted “horizontally,” losing its overlap with the bare response curve, as $g$ increases. This shift is the physical reason underlying the rapid drop in the quantity $\overline{n}_{\rm ph0}$ (i.e. for driving at $\omega=\omega_{0}$). Nevertheless, we stress again that driving at frequency $\omega$ can only produce a response at frequencies $n\omega$, with no spectral-weight transfer to neighboring frequencies: Fig. 4 was prepared by driving at every individual frequency $\omega$ represented on the $x$-axes and by solving Eqs. (22) and (23) to obtain the NESS at each $\omega$. This does mean that Fig. 4 can be used to pose the question of whether any self-blocking is present if the laser frequency is adjusted to resonance with the peak of the phononic response at each value of $g$. In this situation one observes that phonons lying outside the spin excitation band undergo only a small reduction with increasing $g$ (factors of 2 or less up to $g=0.5J$), whereas phonons lying inside the band are suppressed by 1-2 orders of magnitude due to their hybridization with many two=triplon scattering states.

IV Hybrid phonon-bitriplon states

As the discussion of Sec. IIIB made clear at the qualitative level, the hybridization of phononic and magnetic excitations is fundamental to the properties of any spin-phonon-coupled system, and thus to the phenomenology of magnetophononic driving. To elucidate the nature of the states created by the coupling of the driving phonon to triplon pairs, we proceed to a quantitative analysis of the frequency shifts and consider the extent of hybridization within these composite collective entities.

Returning to the features of Fig. 4, we have already remarked on the appearance of two mutually repelling excitation features when the system is driven at a phonon frequency close to the band edge. Because the second set of excitations is not evident in Figs. 4(a) and 4(g), where the in-band phononic response to an out-of-band $\omega_{0}$ is very weak, in Figs. 5(a) and 5(b) we show the same data on a logarithmic intensity scale, confirming the presence of a second peak that grows with $g$. In Figs. 5(c) and 5(d) we show the analogous data obtained for $J^{\prime}$-models with the same range of $g^{\prime}/J^{\prime}$ values. Clearly, as in Sec. IIIA, the physics of hybrid states in the $J^{\prime}$-model is identical to the $J$-model up to minor quantitative details (that can be traced to the matrix elements of the triplon pair-creation and -annihilation terms in Eqs. (12) and (13)), and thus we do not discuss the $J^{\prime}$-model further in this section.

A key additional observation concerning the second set of excitation features is that they appear for small $g$ at the band-edge frequencies, $2\omega_{\rm min}$ in Figs. 4(a) to 4(d), 5(a), and 5(c) and $2\omega_{\rm max}$ in Figs. 4(f) to 4(i), 5(b), and 5(d), before being repelled further from the bare phonon frequency as $g$ is increased. Thus one obtains a picture of “magnetic” hybrid states being induced within the spin sector by the influence of the driven “phononic” state, before stronger $g$ values cause a strong admixture of lattice and spin character.

One way to confirm the hybrid nature of these states is to begin in the regime where weak hybridization is guaranteed. In Fig. 6(a) we show the phonon and triplon occupations obtained when the driving frequency, $\omega_{0}=0.6J$, lies far below the two-triplon band. While $n_{\rm ph0}(\omega)$ undergoes only minor changes, indicating that this is a rather well localized phononic mode, its hybridization is clearly strong enough to shift its frequency out of the regime covered by Fig. 3(a). $n_{\rm x0}(\omega)$ indicates that a spin response does emerge with $g$ despite the nonresonant self-blocking, and in Fig. 6(b) we show that this magnetic dressing of the phononic mode contains all the $k$-components of $n_{\rm x0}(\omega)$. In Figs. 6(c) and 6(d) we show the analogous results for a driving frequency, $\omega_{0}=3.0J$, lying well above the two-triplon band, where the hybridization effects are qualitatively similar but quantitatively are much weaker. Here we do not attempt to discern the weak effects of these off-resonant driving processes at the band edges.

Returning to the regime of strong hybridization controlled by $g$, driving frequencies near both band edges are shown in Fig. 4. For a complete analysis of the mutual repulsion, we gather the characteristic frequencies of these phonon and spin spectra in Figs. 4(e) and 4(j), which display clearly the $g^{2}$ evolution in frequency shift expected of hybrid excitations. To quantify the admixture of lattice and spin character, at the lower band edge one may define a hybridization parameter $s=g/|\omega_{0}-2\omega_{\rm min}|$, and similarly with $2\omega_{\rm max}$ for the upper band edge. A language of “phononic” and “magnetic” hybrids is useful at $\omega_{0}=0.6J$ and $3.0J$ (Fig. 6), where $s<1$ for all $g$. However, when $s\approx 10$ both hybrids are strongly magnetic and phononic, and indeed the 50:50 weight distribution evident at larger $g$ in Figs. 4(c), 4(d), 4(h), and 4(i) suggests states that are maximally hybridized. For the hybrids repelled outside the band, the coinciding peaks in $n_{\rm ph0}(\omega)$ and $n_{\rm x0}(\omega)$ identify them as a strongly triplon-dressed version of the “phononic” hybrid shown in Fig. 6. The in-band hybrids lie in a continuum of propagating triplon-pair states, and thus manifest themselves as broader peaks, whose maxima lie at slightly different energies in $n_{\rm ph0}(\omega)$ [Figs. 4(b) and 4(f)] and in $n_{\rm x0}(\omega)$ [Figs. 4(d) and 4(h)].

For the driving and damping of our system, all of the strongly hybridized states are to a good approximation “phonon-bitriplons,” in which each phonon hybridizes with one triplon pair (${\tilde{t}}_{k}^{\dagger}{\tilde{t}}_{-k}^{\,{\dagger}}$) of zero net momentum. For specificity we reserve the term “phonon-bitriplon” for the composite collective entity formed at $s\geq 1$. However, it is clear that all of the hybrid states forming in a system with a spin-phonon coupling of the form described by Eqs. 5 and 6 involve the dressing of phonons by bitriplons, and conversely. We comment that the formation of composites from a boson pair and a single boson of a different kind is not common in condensed matter physics. A more common scenario is composite formation involving a boson and a pair of fermions, for example where a photon and an exciton form a polariton. Scattering processes of one boson into two different bosons are somewhat more familiar, with phonon-bimagnon processes being discussed both theoretically Lorenzana and Sawatzky (1995a, b) and experimentally Grüninger et al. (2000); Windt et al. (2001) in the optical spectra of cuprate quantum magnets. Phonon-bimagnon processes at GHz frequencies are now applied increasingly in the field of magnonics Bozhko et al. (2020), while photon-bimagnon processes are engineered in cavity spintronics Harder and Hu (2018). Similar physics could also be realized using ultracold atoms Eckardt (2017), where the optical lattice blurs the distinction between photon and phonon, although we are not aware of an experiment. On a very different energy scale, in particle physics the virtual decay of the Higgs boson into pairs of W or Z bosons Buchmüller and Lüdeling (2006); The CMS Collaboration (2019) is an off-shell process with intermediate $s$, where the level repulsion of Figs. 4(e) and 4(j) is known as a “Higgs mass renormalization.”

V Spin-band engineering

In Sec. IV we have shown how the driven phonon creates an additional hybrid state in the spectrum of the system, but we have not yet shown whether phononic driving can alter the excitation spectrum (i.e. the effective bulk properties) of the spin sytem. At first sight one might think that spin-band modification is not possible while the driven phonon is completely harmonic, and that this requires the anharmonic terms considered in nonlinear phononics. However, we will show that the coupling to the magnetic subsystem introduces an intrinsic non-linearity, and thus that the spin spectrum can indeed be altered to a significant extent by a harmonically driven optical phonon.

V.1 First-order band engineering

The leading order of a Magnus expansion Blanes et al. (2009) for time-dependent Hamilton operators consists of the time-averaged Hamiltonian. If one considers the $J$-model, Eqs. (25b) and (25c) describe oscillations of $v_{k}$ and $w_{k}$ at an average frequency $\omega_{k}+gq_{0}y_{k}$ that differs for each wave vector (neglecting higher-order corrections in $g$), and hence the average value of $v_{k}$ vanishes. Averaging Eq. (20b) implies that

[TABLE]

and hence that

[TABLE]

in the $J$-model, meaning that the driving causes the equilibrium (stationary) position of the phonon to be displaced by a finite amount, $q_{0}$. By similar considerations for the $J^{\prime}$-model we obtain

[TABLE]

In Fig. 7(a) we compute $q_{0}$ in the $J$-model, using a laser frequency resonant with the driving phonon ($\omega=\omega_{0}$, for which we use the notation ${\overline{q}}_{0}$), and in Fig. 7(b) we show our results for ${\overline{q}}^{\prime}_{0}$ in the $J^{\prime}$-model.

In contrast to the results of Secs. III and IV, the stationary displacements differ dramatically between the two models. The values of ${\overline{q}}_{0}$ are always negative and reach 4% of the lattice dimension at their maximum extent, which is one order of magnitude larger than the values of ${\overline{q}}^{\prime}_{0}$. The phonon frequencies most effective in controlling ${\overline{q}}_{0}$ are neither those at the band edges, where giant self-blocking suppresses the displacement almost completely, nor those at the band center, where the driving terms decouple, but the “quarter-band” ones around $k=\pi/4$ and $3\pi/4$. A broader range of mid-band frequencies is effective in modulating ${\overline{q}}^{\prime}_{0}$, where the essential qualitative feature is a change in the sign of the displacement as the driving frequency is moved through the band center. The origin of these differences can be found by a closer inspection of the differential equations. In $\mathcal{U}_{J}$ the matrix element is $y_{k}$, while in $\mathcal{U}_{J^{\prime}}$ it is $y^{\prime}_{k}$, which changes sign at $k=\pi/2$. The evolution of both quantities, shown by the dashed white lines in Fig. 7, indicates that their dominant contributions result from the resonant wave vector, $k_{\text{res}}$, selected by the driving phonon according to $\omega_{0}=2\omega_{k_{\text{res}}}$. The excellent agreement between the two sets of curves in Fig. 7 illustrates clearly the validity of Eqs. (31) and Eqs. (32).

These stationary phonon displacements have a direct effect on the spin bands by modifying the magnetic interactions. Here we consider only the linear-order terms, which in the $J$-model cause the change $J\to{\tilde{J}}=J(1+gq_{0})$, and hence renormalize the triplon dispersion [Eq. (10)] to

[TABLE]

while in the $J^{\prime}$-model $J^{\prime}\to{\tilde{J}}^{\prime}=J^{\prime}(1+g^{\prime}q_{0})$ yields the renormalized dispersion

[TABLE]

At this linear level, the triplons retain a cosinusoidal dispersion and the band-engineering effects of the phonon separate cleanly. A phonon mode coupling strictly to the intradimer bond of an alternating spin chain will renormalize the band center downwards without changing the band width, and a phonon mode coupling strictly to the interdimer bond will renormalize the band width, upwards or downwards, without changing the band center.

With a view to later detection (Sec. V.3), we define the lower and upper edges of the renormalized bands as

[TABLE]

and in Fig. 8 we show the evolution of the two-triplon band extrema with $g$ and $g^{\prime}$ for two phonon frequencies chosen near the lower and upper maxima of $|{\overline{q}}_{0}|$ [Fig. 7(a)]. Qualitatively, in the $J$-model we observe the downwards renormalization of the band center by both phonons contrasting with a weak band-narrowing (for $\omega_{0}=2.2J$) or band-broadening (for $\omega_{0}=1.7J$) in the $J^{\prime}$-model. Quantitatively, the $\omega_{0}=2.2J$ phonon is more effective in the $J$-model and the $\omega_{0}=1.7J$ phonon in the $J^{\prime}$-model, but far the most important observation is that our standard driving leads to percent effects on the band center [$J$-model, Figs. 8(a) and 8(c)] but only hundredths of a percent on the band width [$J^{\prime}$-model, Figs. 8(b) and 8(d)].

In Fig. 9 we investigate the scope for increasing band-engineering effects by increasing the laser electric field (mindful of the increased thermal management this mandates). For the chosen value of $g$ and $g^{\prime}$, we find that the two-triplon band extrema display a quadratic evolution only as far as our standard driving value, $E_{0}=0.2\gamma$, before entering an extended regime of largely linear dependence. Only the band center for driving by a phonon with $\omega_{0}=1.7J$ displays a more complex behavior, due presumably to extreme feedback effects setting in at very strong fields. We remark that such strong fields, used as ultrashort pulses to control the system temperature, can yield 10% effects on the engineered band center and 1% effects on the band width.

We conclude our survey of linear-order spin-band engineering by examining its dependence on the spin damping coefficient, $\gamma_{\rm s}$. To do this, in Fig. 10 we show not the band extrema but the stationary phonon displacements ${\overline{q}}_{0}$ and ${\overline{q}}^{\prime}_{0}$. While the qualitative behavior of the driven stationary displacement is not changed by $\gamma_{\rm s}$, it is clear that the quantitative extent follows an approximate $1/\gamma_{\rm s}$ trend, and hence that the observation of spin-band engineering will require systems in which the intrinsic damping of the spin modes is weak. This result also underlines the fact that the finite driven values of ${\overline{q}}_{0}$ and ${\overline{q}}^{\prime}_{0}$ result from nonlinearities introduced by the coupling to the spin system.

V.2 Second-order band engineering

For completeness we analyze the second-order term in the Magnus expansion. While every order of the expansion can be computed systematically Blanes et al. (2009), the convergence is rapid in the parameter regime of our present considerations and a full understanding of our calculated results can be obtained from the first- and second-order terms only. We consider the minimal Hamiltonian

[TABLE]

where $t_{+}$ creates a triplon of flavor $\alpha$ at $k\neq 0$ and $t_{-}$ a triplon of flavor $\alpha$ at $-k$. The oscillating amplitude of the cross-term is $\mu(t)=gy^{\prime}_{k}q(t)$ in the $J$-model and $\mu(t)=-gy^{\prime}_{k}q(t)/\lambda$ in the $J^{\prime}$-model; we assume that $\mu(t)=\mu_{0}\cos(\omega t)$ for a suitably chosen time offset, which in view of the dominant cosinusoidal behavior of $q(t)$ Yarmohammadi et al. (2021) is well justified. In the interaction picture, the dynamics of the diagonal terms is contained in the operators and the time-dependent Hamiltonian takes the form

[TABLE]

The time-averaged action can be computed by the Magnus expansion as an asymptotic series in $\mu_{0}$, where as in Subsec. V.1 the first order is given by the time average over one period ($T=2\pi/\omega$). However, away from resonance, meaning that the system is not driven near $\omega=2\omega_{k}$, there is no first-order term, and this situation arises if one considers the changes to the band edges ($\omega_{k}=\omega_{\text{min}}$ or $\omega_{k}=\omega_{\text{max}}$) caused by driving the system with a phonon whose frequency lies well inside or outside the band.

The general form of the second-order term is

[TABLE]

Here we insert

[TABLE]

and perform the inner integration over $t_{2}\in[0,t_{1}]$ to obtain

[TABLE]

where $\widehat{B}$ denotes $[t_{+}t_{-},t^{\dagger}_{+}t^{\dagger}_{-}]=t^{\dagger}_{+}t_{+}+t^{\dagger}_{-}t_{-}+1$. The $t_{1}$-integration [Eq. (38)] effects a time average in which only the $\cos^{2}(\omega t_{1})$ term contributes $1/2$, whereas all the other combinations vanish, and thus we derive the second-order correction

[TABLE]

The corresponding energy shift of the two-triplon excitation band is

[TABLE]

which by continuity will extend also to the band edges at $k=0$ and $k=\pi$.

For driving in resonance with an in-band phonon, $2\omega_{\text{min}}<\omega=\omega_{0}<2\omega_{\text{max}}$, we deduce from Eq. (42) that $\delta\omega^{(2)}_{\text{min}}$ is positive whereas $\delta\omega^{(2)}_{\text{max}}$ is negative, and thus that the second-order contribution is a band-narrowing. In Fig. 11 we show this correction as function of the driving laser electric field at a fixed value of $g=0.3=g^{\prime}$. The highly non-linear, and even non-monotonic, behavior is a consequence of the fact that the oscillation amplitude, $q(t)$, is contained within the quantity $\mu_{0}$. Nevertheless, these second-order effects are very small when compared with the first-order ones (Subsec. V.1), and for this reason we do not compute any higher orders. However, we will show next that the second-order corrections are indeed detectable in our calculations.

V.3 Detection of tailored spin bands

We turn now to the question of how to measure driving-induced spin-band renormalization in the magnetophononic protocol. Because the pump electric field, at frequency $\omega$, is required to excite the target phonon(s) creating the desired NESS, a full characterization of the NESS properties will require the introduction of an additional frequency. Thus we introduce a further field component,

[TABLE]

where as before $E_{0}$ is the strong pump drive and now $E_{1}$ is a significantly weaker “probe” drive, represented by the small green waves in Fig. 1. Because we are investigating NESS, both field components are continuous and the time delay used in true pump-probe studies is absent.

In the remainder of this subsection we consider a single driving phonon whose frequency is located within the two-triplon band ($2\omega_{\text{min}}<\omega_{0}<2\omega_{\text{max}}$) and pump it resonantly ($\omega=\omega_{0}$), while scanning the detection frequency, $\Omega$, around the lower and upper band edges where the driving-induced changes are expected to be clearest. As the most sensitive diagnostic of the edges of the modified two-triplon band, we take the liberty of showing not $n_{\rm x0}(\omega_{0},\Omega)$ but the respective components $u_{k=0,0}$ to characterize the lower band edge and $u_{k=\pi,0}$ for the upper. These quantities are readily computed from the equations of motion of Sec. II and we comment below on current experimental developments in the use of coherent light as a direct probe of magnetic phenomena.

Figure 12 shows the value of $u_{k=0,0}$ measured on scanning $\Omega$ through the lower band edge for different values of the spin-phonon coupling in both the $J$- and the $J^{\prime}$-model. In contrast to Figs. 4(c), 4(d), 4(h), and 4(i), which showed two band-edge peaks in $n_{\rm x0}(\omega)$ when changing the driving frequency in the presence of a near-resonant phonon mode, here the system is driven at $\omega=\omega_{0}=1.7J$ and $2.2J$, $s<1$, and we observe a single band-edge peak in the magnetic response. The shift of this peak away from the equilibrium band edge (shown by the red shading) increases strongly with $g$ in the $J$-model [Figs. 12(a) and 12(c)], where in accordance with Fig. 8 it is downward for both phonons (and, given their relative separation from the lower band edge, surprisingly similar in magnitude). By contrast, the shift increases only weakly with $g^{\prime}$ in the $J^{\prime}$-model, to the extent that its expected change of sign is not discernible between Figs. 12(b) and 12(d) because the shift is so weak in the latter case.

For a fully quantitative interpretation of the probe spectrum, we proceed to perform a systematic variation of all relevant parameters. In Fig. 13 we vary the pump and probe amplitudes (maintaining $E_{1}/E_{0}=0.2$) at fixed $g/J=0.5$ and $g^{\prime}/J^{\prime}=0.5$. While the area under the peaks depends quadratically on $E_{0}$ in both $J$- and $J^{\prime}$-models, an approximately quadratic shift of the peak position is clear only in the $J$-model, where this effect is very strong [Figs. 13(a) and 13(c)]. In the $J^{\prime}$-model, it may be possible to discern a very weak quadratic variation of the peak positions with $E_{0}$, but there are clearly other contributions to their shift from the band edge.

In Fig. 14 we vary the spin damping, finding a minor broading of the resonance peaks in line with general expectations. The reduction of the peak shifts with increasing $\gamma_{\rm s}$ are at first sight less intuitive, but results from the fact that the total response of the spin system is reduced by stronger damping, as already observed in the stationary phonon displacements, ${\overline{q}}_{0}$ and ${\overline{q}}^{\prime}_{0}$, shown in Fig. 10.

Finally, in Fig. 15 we vary the probe amplitude. Again we consider both $J$- and $J^{\prime}$-models, but only the $\omega_{0}=2.2J$ phonon, in order to compare the situation with no pump field, $E_{0}=0$, to standard driving. Applied alone [Figs. 15(a) and 15(b)], the weak probe field changes neither the position nor the width of the resonance peak, and only its height depends quadratically on $E_{1}$. However, the peak appears with an additional shift, which we denote simply by $\delta\omega_{\text{s}}$ in each panel of Figs. 15, that results from the weak hybridization of the band-edge states with a phonon whose frequency lies well inside the band. Quite generally, such a shift appears at both ends of the spectrum of the spin-phonon system at equilibrium, as the hybridization assures the addition of a single mode, such that both band edges are shifted outwards. When the phonon frequency is far from the band edges, the band-edge states can be regarded as $s<1$ analogs of the phonon-bitriplon of Sec. IV, but as $\omega_{0}$ is moved systematically towards one or other band edge then a true phonon-bitriplon would form at this edge. The probe field is added to provide the weak driving at band-edge frequencies that is required for a detailed characterization of the equilibrium and nonequilibrium response (meaning without and with $E_{0}$). We comment that the hybridization effect leading to $\delta\omega_{\text{s}}$ is very similar in the $J-$ and $J^{\prime}$-models, as may be expected given the relevant matrix elements, which are $gy^{\prime}_{k}$ and $-g^{\prime}y^{\prime}_{k}/\lambda$. This weak frequency shift is readily accounted for in order to quantify accurately the band-engineering effects produced by applying the pump field [Figs. 15(c) and 15(d)].

Figure 16(a) illustrates, using the example of the upper band edge in a $J$-model, the three effects contributing to the positions of the resonance peaks detected in the pump-probe protocol. It is clear that the first-order shift, which we denote $\delta\omega^{(1)}_{\text{max}}$, is largest in absolute value, and is negative (Subsec. V.1). The next effect by magnitude is the hybridization shift, $\delta\omega_{\text{s}}$, which is positive at the upper band edge. However, analyzing the computed peak positions reveals that these two contributions are not sufficient for an accurate description, whereas the small discrepancy is well accounted for by $\delta\omega^{(2)}_{\text{max}}$.

Figures 16(b) and 16(c) represent the relative signs and sizes of these three contributions at both band edges in the $J$-model, comparing them with an equivalent $J^{\prime}$-model for two different phonon frequencies. Here we introduce additional, self-explanatory notation to distinguish the hybridization shifts of both models at both band edges. This format makes clear that $\delta\omega^{(1)}$ is the largest shift and is always downwards, whereas $\delta\omega^{\prime(1)}$ is smaller and can change sign with the location of the driving phonon in the lower or upper half of the band. By contrast, $\delta\omega_{\text{s}}$ and $\delta\omega^{(2)}$ are always similar to $\delta\omega_{\text{s}}^{\prime}$ and $\delta\omega^{\prime(2)}$, and have the same (opposing) signs in all cases. We remark again that the quantitative results are similar for the upper and lower band edges, which is why we have analyzed $u_{k=\pi,0}$ and $u_{k=0,0}$ interchangeably, instead of doubling the length of our discussion. We remark once more that only $\delta\omega^{(1)}$ and $\delta\omega^{(2)}$ are consequences of the magnetophononic driving, whereas $\delta\omega_{\text{s}}$ is an equilibrium effect of spin-phonon hybridization.

For complementary insight into our pump-probe results, in Fig. 17(a) we show the phonon spectrum, $n_{\rm ph0}(\omega_{0},\Omega)$, matching the spin response of Fig. 15(c), i.e. for probe frequencies around the lower band edge, and in Fig. 17(b) the phononic response for probing around the upper band edge. In the absence of a probe beam, the phonon occupation is essentially flat around the band edges, with no discernible features forming in these regions when the only driving is resonant with the available in-band phonon at $\omega_{0}=2.2J$. However, increasing the probe intensity reveals that the formation of the predominantly magnetic spectral features at the band edge in Figs. 12, 13, 14, and 15 is accompanied by a small dip in $n_{\rm ph0}(\omega_{0},\Omega)$. This weak response indicates the weak phononic character of these hybrid states; the fact that it is negative is a consequence of the removal of phonon energy required in the excitation of the band-edge hybrid.

We conclude our analysis of spin-band engineering by stressing that our focus has been to illustrate all of the qualitative phenomena present when pumping with a light-driven phonon and probing with a separate laser in the same frequency range. For this purpose we have not considered the effects of the phonon coordinate on the magnetic interactions beyond linear order. We have also not dwelled on maximizing the energy shifts we induce, and thus we have shown mostly percent effects. However, we comment that the band shifts exceeding 10% that are found at stronger pump fields (Fig. 13) can for a one-dimensional system mean that well over 50% of the spin spectral weight is shifted completely out of its previous energy range. In the present study we have also used a rather generic spin band, whose gap is approximately equal to its band width, whereas focusing on a particularly narrow-band system would lead to much stronger relative shifts of spectral weight.

VI Strong spin-phonon coupling in quantum magnetic materials

VI.1 CuGeO3

The strong spin-phonon coupling in CuGeO3 was revealed by the fact that it drives a spin-Peierls transition at $T_{\rm sp}=14$ K Hase et al. (1993). While it is clear that the leading physics of this system is a dimerization of the spin chain formed by the strongest (Cu-O-Cu) superexchange bonds, details of the magnon dispersion and of other thermodynamic measurements led to the introduction of both a significant next-neighbor coupling, $J_{2}$, and a non-negligible interchain coupling. A recent ultrafast investigation that used soft x-ray frequencies to probe the low-lying electronic states Paris et al. (2021) also suggested that the observed damping was a consequence of coherently excited phonon modes coupling strongly to short-ranged magnetic correlations.

The crystal structure of CuGeO3 has space group P$mma$, which is reduced to C$mca$ when the spin-Peierls transition enlarges the unit cell (while maintaining the orthorhombic symmetry). IR-active phonons are available over a wide range of energies in the high-temperature structure Popović et al. (1995), and the rather large unit cell of the spin-Peierls phase makes their number significant, although for laser-driving purposes we note that all of the $A_{u}$ modes are silent. Inelastic neutron scattering (INS) has been used to characterize all the phonon modes of the high-temperature phase Braden et al. (1998, 2002), finding the strongest response at frequencies of 3.2 and 6.8 THz, and suggesting that the spin-Peierls transition is of a type occurring without an accompanying soft mode Werner and Gros (1998); Pouget (2001). Here we comment that ultrafast methods appear to offer a qualitatively different approach to the investigation of low-lying phonons around and below the spin-Peierls temperature Spitz et al. . In Fig. 18 we illustrate three phonon modes of the high-temperature structure that have been identified by comparing electronic structure calcuations Spitz et al. with experiment Popović et al. (1995). One of these [Fig. 18(a)] is IR-silent whereas the other two are expected to be promising candidates for coherent laser driving. We note also that most of these phonon modes tend to involve motions of all the atoms in the system, and hence they will have both $J$- and $J^{\prime}$-model character in the language of our two simplifying models (Fig. 1).

The magnetic excitation spectrum, also measured by INS Regnault et al. (1996), shows relatively broad triplon bands, by which is meant that their gap is considerably smaller than their band width. In Figs. 19(a) and 19(b) we show the one-triplon dispersion along and across the chain direction, from which the early two-dimensional fit of Ref. Uhrig (1997) deduced the illustrative superexchange parameters $J=10.7$ meV, which sets the one-triplon band center, $J^{\prime}=8.3$ meV (i.e. $\lambda=0.78$), $J_{2}=0$, and $J_{a}=1.5$ meV (interchain); we remark here that later studies provided a more refined global parameter set Knetter and Uhrig (2001). Figure 19(c) represents the full energy range of the two-triplon excitation spectrum for these parameters and also shows the in-band locations of four IR-symmetric phonon modes of the high-temperature structure, as measured by Ref. Popović et al. (1995). As noted above, the $A_{u}$ mode [Fig. 18(a)] is IR-silent, but the $B_{1u}$, $B_{2u}$, and $B_{3u}$ modes all have similar oscillator strengths Popović et al. (1995). We therefore propose the $B_{1u}$ mode shown in Fig. 18(b) as a good candidate for driving in the lower half of the two-triplon band and the $B_{3u}$ mode of Fig. 18(c) as a good candidate for driving in the upper half of the band, while we do not show the atomic motions in the $B_{2u}$ mode because it is rather close to the two-triplon band center. We stress that, because ultrafast driving experiments are also performed with ultra-intense electric fields, the choice of drivable phonons is by no means restricted to the modes depicted in Fig. 19, and our qualitative message is rather that CuGeO3 remains an excellent candidate material for observing the phenomena we analyze (Secs. III, IV, and V). However, we note also that the relevant experiments do need to be at performed at low temperatures, $T\ll T_{\rm sp}=14$ K, in order to preserve the dimerized state of the system.

VI.2 (VO)2P2O7

As a low-dimensional $S=1/2$ quantum spin system with an excitation gap, (VO)2P2O7 first attracted attention as a candidate to realize the two-leg ladder geometry. However, INS measurements of the triplon dispersion Garrett et al. (1997) soon revealed its nature as a quasi-one-dimensional alternating spin chain, and nuclear magnetic resonance (NMR) revealed a large and complex structural unit cell with alternation in all three lattice directions Kikuchi et al. (1999). Further theoretical analysis then deduced the presence of frustrated interchain coupling, leading to a magnetic model containing two species of dimerized spin chain, lying on alternating planes and with an effective coupling that is weak as a consequence of interchain frustration Uhrig and Normand (2001). For experimental purposes, the intrinsic dimerization of the dominant spin chains avoids the need for temperatures as low as those required in CuGeO3.

The crystal structure of (VO)2P2O7 also has orthorhombic symmetry, with space group P$ca$21, and a very large unit cell (104 atoms). This structure lacks inversion symmetry, and hence all of the phonons may have either IR or Raman character, depending on the polarization of the light. Thus an extremely large number of IR-active phonons is available over the full range of energies. Initial studies of the phonon spectrum by Raman scattering Grove et al. (2000) found by raising the temperature that a very strong renormalization of the phonons by the spin sector takes place. Theoretical fits Uhrig and Normand (2001) suggested that the spin-phonon coupling should be very strong, $g\simeq 0.5J$, which was the basis for our extending the analyses of Secs. III, IV, and V to this value of $g$ and $g^{\prime}$.

Figure 20 shows the atomic displacements in three phonon excitations with relatively large oscillator strengths, which were found in exploratory electronic structure calculations of (VO)2P2O7 Colonna . These calculations used Quantum Espresso Giannozzi et al. (2017) to obtain a stable and insulating structural solution by assuming antiferromagnetic order, but were neither optimized for correlation effects on the V ions nor compared to any experiments. Nonetheless, they do illustrate the key features of (i) a very large number of phonon modes and (ii) a general increase in the average oscillator strength of these phonons with frequency, arising because modes with larger dipole moments are stiffer. Despite the significant quantitative mismatch in frequencies, it is not unrealistic to suggest that the 121 cm*-1* (3.63 THz) $B_{1}$ mode shown in Fig. 20(a) and the 167 cm*-1* (5.01 THz) $B_{1}$ mode shown in Fig. 20(b), which has a weaker $A_{1}$ partner at 162 cm*-1* (4.85 THz) in the calculation, are candidates for the 70 and 123 cm*-1* modes [2.10 and 3.69 THz, the latter with a partner at 118 cm*-1* (3.54 THz)] that appear most strongly in experiment Grove et al. (2000). One may anticipate that these two modes are the most suitable low-frequency candidates for coherent laser driving, but we stress again that many suitable modes are present at higher frequencies, including the 212 cm*-1* (6.35 THz) $A_{1}$ mode shown in Fig. 20(c).

Following the fit of Ref. Uhrig and Normand (2001), Figs. 21(a) and 21(b) show the one-triplon dispersions along and across the direction of the alternating chains in the two inequivalent planes (which we denote $A$ and $B$). Again both triplon bands are broad, with gaps rather smaller than their band widths Garrett et al. (1997), and with fitting parameters $J_{A}=12.3$ meV, $J^{\prime}_{A}=8.1$ meV (i.e. $\lambda_{A}=0.66$), $J_{aA}=1.0$ meV, $J_{bA}=1.4$ meV (the latter pair mutually frustrating interchain interactions) and $J_{B}=10.4$ meV, $J^{\prime}_{B}=8.0$ meV (i.e. $\lambda_{B}=0.77$), $J_{aB}=1.1$ meV, $J_{bB}=1.6$ meV. Figure 21(c) represents the full energy range of the two two-triplon excitation spectra, assuming only processes involving pairs of A and pairs of B triplons, and again shows the locations of the three phonons whose normal modes are displayed in Fig. 20. We conclude that (VO)2P2O7 should offer an excellent materials platform for realizing all three of the magnetophononic phenomena revealed in our study.

VI.3 Ultrafast band-engineering experiments

Experiments designed to follow the theoretical protocal of creating true spin NESS, at frequencies resonant with the spin spectrum and in bulk-driven quantum magnets, require a thin-film geometry and very efficient thermal transfer in order to maintain a low sample temperature Yarmohammadi et al. (2021). In principle, self-blocking allows some relaxation of the constraints on pump intensity, driving time, and sample thickness, although in practice strong electromagnetic driving can induce heating by a variety of channels. An ultrafast pulsed protocol avoids extreme heating problems through the very short driving time, but usually involves strong electric fields and samples of $\mu$m up to mm thicknesses, and thus the pulse repetition time should remain long. In the context of modifying the spin excitation spectrum, we comment that ${\tilde{J}}(q_{0})$ and ${\tilde{J}}^{\prime}(q_{0})$ are in general highly nonlinear functions of $q_{0}$, and where conventional experimental probes usually require only a low-order expansion, coherent laser driving can produce very large $q_{0}$ values Giorgianni et al. .

A final issue concerns the optimal type of experiment to perform. Coherent light remains a rather insensitive direct probe of magnetism, and experiments performed to date, such as absorption, reflection, and polarization rotation (birefringence), probe only some effects of the lattice that reflect the spin-phonon coupling. To identify other probes of novel magnetic states applicable to a complete analysis of a driven CuGeO3 or (VO)2P2O7 system, we will continue to review the rapidly evolving technological developments in the measurement of quantities such as magnetic circular dichroism, the magneto-optic Kerr effect, or second-harmonic polarimetry as time-resolved variants become available across an increasing spectrum of probing frequencies (up to and including x-rays). The formalism of Sec. II remains fully applicable to the expectation values measured by these more direct probes of magnetic order, correlations, and excitations.

VII Discussion and Conclusion

Ultrafast laser technology has enabled qualitative advances in the study of nonequilibrium phenomena in correlated condensed matter. Extending the reach of ultrafast driving methods to the rich variety of complex many-body states available in quantum magnetic materials requires overcoming the fact that the direct coupling of light to spin is generally rather weak, and thus inefficient. For this purpose we investigate the magnetophononic channel, in which the driving laser couples to an infrared-active optical phonon and the associated lattice displacements modulate the magnetic interactions. This approach offers highly frequency-specific driving possibilities by exploiting the resonances both between laser and phonon and between the driven phonon and the excitations of the quantum spin system. Intense driving electric fields and strong spin-phonon coupling then allow one to probe the properties of a correlated quantum magnet driven far from its equilibrium state.

The characteristic energy scales in quantum magnetic materials are typically rather low, making their quantum many-body states very sensitive to heating, as a result of which no serious analysis can avoid taking the energy flow into account. To model the problem of a quantum magnet with both driving and dissipation through the lattice phonons, we adopt the Lindblad treatment of open quantum systems and analyze the non-equilibrium stationary state (NESS) of system subjected to continuous laser driving. Within this framework we consider a very straightforward example of a gapped spin system, the dimerized spin chain, whose elementary excitations are triplons that can be treated as conventional (rather than hard-core) bosons for sufficiently strong dimerization. We adopt a similarly minimal model for the driven phonon, namely an infrared-active optical mode coupled to only to the strong bonds ($J$-model) or only to the weak bonds ($J^{\prime}$-model) of the dimerized chain.

Having previously used this minimal magnetophononic model to establish the framework for the weak-coupling, or linear-response, regime of the $J$-model Yarmohammadi et al. (2021), the primary focus of our present study is the regime of strong feedback, or back-action, of the driven spin system on the driving phonon. Particularly when the phonon frequency is chosen close to the band edges of the two-spin excitation spectrum, a strong spin-phonon coupling causes strong hybridization into composite collective states whose characteristic frequencies differ significantly from those of their constituents. In the model we study, these collective states are phonon-bitriplons, a somewhat rare example of a composite formed from three bosons of two different types.

In a NESS experiment, the shift of characteristic frequencies causes a dramatic “self-blocking” effect, by which the spin system acts as a strong negative feedback on the driven phonon, pushing its resonance off the driving frequency and thus drastically suppressing the effective driving field. Only in an experiment with a range of driving frequencies, such as those present in an ultrashort pulse, would one observe the shifts of spectral weight associated with the level-repulsion caused by the spin-phonon hybridization. Even in this situation, however, the self-blocking caused by strong mixing with off-resonant spin levels remains significant. We comment here that our present study retained the NESS protocol for all possible driving frequencies, and did not include the intense and instantaneous nature of an ultrashort pulse.

While driving phonons resonant with the two-spin band edges is an excellent way to create composite collective hybrid states, an important consequence of self-blocking is that it is not a very efficient way to engineer the bulk properties of a quantum magnet. Optical control is the only technique available to switch these properties on the characteristic timescales of the spin system, and our analysis reveals the important insight that the frequencies most efficient for this purpose lie in specific regions within the two-spin excitation spectrum. The dominant band-engineering effects arise at linear order as a consequence of a stationary displacement of the driven phonon, which results from the steady population of excited triplons created by its action on the spin system. While the second-order contribution is weak, it is also detectable for the typical parameters of a phonon-driven quantum magnet, and thus is required for a quantitative description of experiment. To detect spin-band engineering, we introduce a weak “probe” beam, which in the NESS protocol is actually a further continuous driving term, applied in addition to the pump but at a completely independent frequency. (We remind the reader that the focus of our present study is not on conventional pump-probe physics, which use a time delay between pump and probe pulses to investigate transient phenomena at switch-on.) This technique yields clear additional signals in the phonon and spin response at the band-edge frequencies of the renormalized, or optically engineered, band. Applying a band-edge probe electric field to a system with a driven mid-band phonon reveals an extra intrinsic frequency shift in the detected signal, caused by off-resonant spin-phonon hybridization, that should also be included in any quantitative analysis.

Our minimal quantum magnetic model contains only two types of bond, within ($J$) and between ($J^{\prime}$) the spin dimers, and coupling the driven phonon to each bond type separately yields some valuable insight. Certain aspects of the response, which with strong spin-phonon coupling is dominated by the feedback between the two sectors, are very similar, in particular the self-blocking and the formation of composite collective states at the band edges. This can be traced to the fact that the matrix elements for the driven phonon to create spin excitation are the same, up to a sign, when the relative couplings ($g/J$ and $g^{\prime}/J^{\prime}$) are the same. By contrast, the band-engineering effect of the driven phonon is completely different between the two situations, which is a consequence of how the back-action from the spin system modifies the phonon. Only in the $J$-model does the driving produce a rather large stationary shift of the equilibrium atomic displacement, whereas in the $J^{\prime}$-model this effect is at least an order of magnitude weaker; once again this behavior can be traced to the matrix elements in the relevant equations of motion. Although it is easy to conclude that phonons coupling to the strong bonds are better suited for spin-band engineering, in fact the difference is also qualitative, in that modulating the intradimer bond changes the band center whereas modulating the bonds between dimers can be used to alter the band width. Although the effect on the band center is always a reduction, the band width can be renormalized both upwards or downwards by choosing the frequency of the driving phonon.

In the present work we have focused on a spin chain as a representative quantum magnet. However, the generic features of the phenomenology we unveil are not restricted to spin chains, and could also be found in spin ladders and in valence-bond states in two and three dimensions. One important ingredient of our model is the sharp peaks in the density of states at both band edges, which concentrates the strongest response of the spin system to two narrow ranges of frequency, and this property can be found in many quantum magnets with a spin gap or with narrow bands, such as those induced by frustrated couplings. Even for magnetically ordered states in three dimensions, there is a large jump in the density of states at the upper band edge Deltenre et al. (2021). Band engineering, which relies on mid-band rather than band-edge phonons, is less dependent on the structure of the density of states, and is a more direct consequence of strong spin-phonon coupling at the strongest bonds of the system. The concept of phononic driving is of course applicable throughout condensed matter, and future work in quantum magnetism can be expected to investigate its manifestations in itinerant as well as in localized systems, in ordered as well as in non-ordered magnets, and in systems with the wealth of complex forms of order found only in magnetic materials, including $3Q$ textures, quadrupolar order, chiral order, nematic order, and still others.

Here we have restricted our considerations to linear magnetophononics, in that we consider only a single driven phonon to explore the leading nonequilibrium phenomena. At this level, the dominant effects are produced by the $q_{0}$-linear correction to the magnetic interactions, and we do not consider higher-order terms in $J(q_{0})$ (the second-order correction we discuss appears in $q(t)$, generated by the equations of motion). Nonlinear magnetophononics Giorgianni et al. considers the simultaneous effects of two or more driving phonons, and in this situation the quadratic correction to $J(q_{a},q_{b})$ contains terms modulating the magnetic interactions at frequencies of $2\omega_{a}$, $2\omega_{b}$, $\omega_{a}+\omega_{b}$, and $\omega_{a}-\omega_{b}$. The sum and difference frequencies enlarge very considerably the range of possibilities available for matching the driving frequency to the characteristic energies of the spin system. In particular, phonon difference frequencies are the key to magnetophononic driving in systems where the spin energy scale is very small, which is the situation relevant to a high percentage of quantum magnetic materials. Observing the signals of such nonlinear driving require one or both of strong spin-phonon coupling and ultra-intense electric fields.

Finally, we have discussed two materials matching our models, to show that the phenomena we identify should be detectable in CuGeO3 and especially in (VO)2P2O7. However, a primary task of theoretical physics is to build models matching the materials on which experiments are performed. For this task, the minimal magnetophononic model we have constructed is extremely versatile, in that more complex spin systems (for example higher-dimensional, ordered, or with anisotropic interactions), more complex phonons (for example dispersive, multiple, or coupling to both $J$ and $J^{\prime}$ bonds), and also more complex dissipative processes (in particular spin-conserving ones) are readily accommodated within the Lindblad formulation to generate qualitatively similar equations of motion. As noted above, although we have considered only the expectation values of the most fundamental phonon and triplon operators, all more complex observables are formed from these and hence our framework is easily adapted to compute the quantities probed by experiment, and in particular those obtained from direct optical probes of magnetic correlations.

Acknowledgements.

We are indebted to N. Colonna, S. Das, and L. Spitz for performing DFT and phonon spectrum calculations, and for their assistance in producing Figs. 18 and 20. We thank S. Behrensmeier for assistance with Figs. 19 and 21, and F. B. Anders, D. Bossini, K. Deltenre, F. Giorgianni, Ch. Rüegg, L. Spitz, and R. Valentí for helpful discussions. We are grateful to the German Research Foundation (DFG) for financial support through projects UH 90-13/1 and B8 of ICRC 160, as well as to the Mercator Research Center Ruhr for support through the Mercur Kooperation Ko-2021-0027. We acknowledge the National Science Foundation for financial support through award numbers DMR-1945529, PHY-1607611, and PHY-1748958 and the Welch Foundation for support through award number AT-2036-20200401. This project was partially funded by The University of Texas at Dallas Office of Research and Innovation through the SPIRe program.