1 INTRODUCTION

The efficient generation of coherent spin waves in antiferromagnets is a relevant aim for modern spintronics and magnonics. Short laser pulses with durations about 100 fs and shorter provide one of the most efficient methods for the generation of coherent magnons in a wide frequency range from gigahertzs to terahertzs. The most universal mechanism of the laser-pulse excitation of coherent magnons is based on inverse magneto-optical (optomagnetic) Faraday and Cotton–Mouton effects, where the polarized light pulse generates an effective magnetic field pulse in magnetic materials and thus induces homogeneous spin oscillations, which correspond to coherent magnons in the center of the Brillouin zone [1]. The maximum frequency of the magnon that can be excited due to these mechanisms is equal to the maximally possible difference between the Fourier frequencies of the effective magnetic field pulse.

The parametric excitation of magnons is a fundamentally different mechanism, where a pump pulse with the frequency ωp excites modes with the half frequency ωp/2. Such an excitation is a quite efficient mechanism for the generation of coherent magnons [2, 3], phonons [46], and silent phonons, which do not interact with light in another way [7]. The theory of parametric effects in ferromagnets is quite developed [812]. The parametric excitation of spin waves in antiferromagnets was also studied in [2, 1318]. It is of interest that the parametric pumping condition allows the excitation of pairs of magnons with opposite nonzero wave vectors. As a result, unlike the magneto-optical Faraday and Cotton–Mouton effects, the parametric pumping allows the excitation of nonhomogeneous spin oscillations, which correspond to coherent magnons deep in the Brillouin zone.

The possibility of exciting coherent magnons with a nonzero wave vectors can lead to significant differences in the description of the spin dynamics. An antiferromagnet is usually describe as two antiferromagnetically coupled ferromagnets with the magnetizations M1 and M2. Then, the magnetic order can be described using the antiferromagnetic vector L = \({{{\mathbf{M}}}_{1}} - {{{\mathbf{M}}}_{2}}\) and the magnetization M = M1 + M2 (M = 0). The spin dynamics that corresponds to magnons in the center of the Brillouin zone and has a small amplitude is successfully described in terms of \({\mathbf{L}}\) in the approximation |L| = const. The frequency of extremely possible spin oscillations, which are determined by the exchange interaction, increases with the wave vector. The modulation of the exchange interaction by the electric field of the laser pulse ensures parametric instability conditions, whereas the inverse Cotton–Mouton effect produces an external torque that acts on spins and results in the canting of magnetic sublattices M1 and M2 from the equilibrium position. This mean that |L| can no longer be considered as constant. In this case (i.e., at |L| ≠ const), the dimension of the system of nonlinear spin dynamics equations for the antiferromagnet is doubled from 3 to 6 or from 2 to 4 for the spherical coordinates. In this work, we theoretically analyze the spin dynamics with |L| ≠ const, which is induced by a short laser pulse due to the parametric excitation of magnons with a nonzero wave vector.

2 TWO-SUBLATTICE HEISENBERG FERROMAGNET

We consider a two-sublattice antiferromagnetic crystal with a cubic symmetry (point group \(m\overline 3 m\)), which is described by the Hamiltonian

$${{H}_{0}} = J\sum\limits_{\langle i,j\rangle } {{{\mathbf{S}}}_{i}}{{{\mathbf{S}}}_{j}} + K\sum\limits_i (S_{y}^{2}S_{z}^{2} + S_{x}^{2}S_{z}^{2} + S_{x}^{2}S_{y}^{2}),$$
(1)

where J > 0 is the exchange interaction energy; K < 0 is the magnetic anisotropy energy, which is negligibly low in the case of RbMnF3 [19]; and S is the spin (for RbMnF3, S = 5/2).

In the spherical coordinate system, the magnetizations of the sublattices of the antiferromagnet are defined as \({{{\mathbf{M}}}_{i}} = M(\sin {{\vartheta }_{i}}\cos {{\varphi }_{i}},\sin {{\vartheta }_{i}}\sin {{\varphi }_{i}},\cos {{\vartheta }_{i}})\), where \(i = 1,2\) denote oppositely directed magnetic sublattices, \(\vartheta \in [0,\pi ]\) is the polar angle measured from the z axis, and \(\varphi \in [0,2\pi ]\) is the azimuth angle in the xy plane measured from the x axis, as shown in Fig. 1a. The minimization of the Hamiltonian (1) determines the ground state for which the spherical angles are \(\vartheta _{1}^{{(0)}} = \arccos (1{\text{/}}\sqrt 3 )\) and \(\varphi _{1}^{{(0)}} = \pi {\text{/}}4\) and \(\vartheta _{2}^{{(0)}}\) = π – \(\arccos (1{\text{/}}\sqrt 3 )\) and \(\varphi _{2}^{{(0)}} = 3\pi {\text{/}}4\).

Fig. 1.
figure 1

(Color online) (а) Antiferromagnetic vector L directed along the \([111]\) axis in the spherical coordinate system and its projection on the xy plane. (b) Canting of the oppositely directed magnetizations of the sublattices M1 and M2 at the angle ρ to the z axis. (c) Orthogonality coefficient \(\eta = {\text{|}}\langle {{{\mathbf{n}}}_{j}}{{{\mathbf{n}}}_{{j + \delta }}}\rangle {\text{|}} = \left| {{{{\left( {\frac{{ - {{\rho }^{2}} - {{\beta }^{2}}{{{\sin }}^{2}}\vartheta }}{2}} \right)}}^{S}}} \right| \ll 1\) versus the canting angle of the sublattices ρ at \(\beta = 0\). (d) Color (hφ, hρ) map of the parametric amplification of the inverse Cotton–Mouton effect, where hφ and hρ are the parametric coefficients in Eqs. (12). The boundary of the black region corresponds to the threshold of the parametric instability.

3 LIGHT-INDUCED SPIN DYNAMICS IN TERMS OF SPIN CORRELATION FUNCTIONS

The collinear magnetizations of the sublattices in a nonequilibrium state are canted due to the dynamic modulation of the exchange interaction by the electric field \({\mathbf{E}} = {{E}_{0}}(\cos \alpha ,\sin \alpha ,0)\) of short laser pulses propagating in the antiferromagnetic crystal along the z axis, where α is the angle between the electric field E in the xy plane and the x axis. Below, we consider the angles of the sublattices

$$\begin{array}{*{20}{c}} {{{\vartheta }_{1}} = \vartheta - \rho ,} \\ {{{\vartheta }_{2}} = \pi - \vartheta - \rho ,} \\ {{{\varphi }_{1}} = \varphi + \beta ,} \\ {{{\varphi }_{2}} = \varphi + \pi - \beta ,} \end{array}$$
(2)

which are canted by small angles \(\rho \ll 1\) and \(\beta \ll 1\). It will be shown below that the ellipticity induced by the pump pulse is sensitive only to the angle \(\rho \) shown in Fig. 1b.

To describe the spin dynamics, it is necessary to know the energy of the interaction of the magnetic subsystem of the antiferromagnet with the electric field of the laser pulse (it is a Raman process). This energy appears in the Lagrangian of the system (see below) and can be represented in the form

$${{U}_{{{\text{sp}}}}}(\vartheta ,\varphi ,\rho ,\beta ) = \mathop {\sum \,'}\limits_i \,\sum\limits_{\delta = 1}^6 {{\chi }_{{nmkl}}}E_{n}^{{(1)}}E_{m}^{{(2)}}\langle S_{k}^{i}S_{l}^{{i + \delta }}\rangle ,$$
(3)

where summation is performed over all magnetic sites (i) and nearest neighbors (δ); \({{\chi }_{{nmkl}}}\) is the fourth-rank tensor for the cubic symmetry (\(m\overline 3 m\)) specified by only three nonzero parameters \({{a}_{1}} = {{\chi }_{{11}}} = {{\chi }_{{22}}} = {{\chi }_{{33}}}\), a2 = χ12 = \({{\chi }_{{21}}} = {{\chi }_{{13}}} = {{\chi }_{{31}}} = {{\chi }_{{23}}} = {{\chi }_{{32}}}\), and a3 = χ44 = \({{\chi }_{{55}}} = {{\chi }_{{66}}}\), where a short notation is used for indices of the tensors; and \({{{\mathbf{E}}}^{{(1)}}}(t)\) and \({{{\mathbf{E}}}^{{(2)}}}(t)\) are the electric fields of ultrafast laser pulses, which lie in the xy plane and have frequencies ω1 and ω2, respectively, whose difference ωp = ω1 – ω2 is in parametric resonance with spin dynamics frequencies.

To calculate spin correlation functions in Eq. (3), we use the method of spin coherent states [20]:

$$\langle S_{k}^{i}S_{l}^{{i + \delta }}\rangle = {{S}^{2}}\langle n_{k}^{i}n_{l}^{{i + \delta }}\rangle ,$$
(4)

where \({\text{|}}{\mathbf{n}}\rangle = {{e}^{{i\vartheta ({\mathbf{n}} \times z){\mathbf{S}}}}}{\text{|}}z\rangle = (1 + \,{\text{|}}\mu {{{\text{|}}}^{2}}){{e}^{{{{\mu }_{j}}{{S}_{ - }}}}}{\text{|}}z\rangle \) and μ = \({{e}^{{i\varphi }}}{\text{tan}}\frac{\vartheta }{2}\). It is known that spin coherent states \({\text{|}}{{{\mathbf{n}}}_{1}}\rangle \) and \({\text{|}}{{{\mathbf{n}}}_{2}}\rangle \) are nonorthogonal and \(\langle {{{\mathbf{n}}}_{1}}{{{\mathbf{n}}}_{2}}\rangle = {{\left( {\frac{{1 + {{{\mathbf{n}}}_{1}}{{{\mathbf{n}}}_{2}}}}{2}} \right)}^{S}}\); this property complicates the calculations of matrix elements. However, in our case, \(\eta = {\text{|}}\langle {{{\mathbf{n}}}_{j}}{{{\mathbf{n}}}_{{j + \delta }}}\rangle {\text{|}}\) = \(\left| {{{{\left( {\frac{{ - {{\rho }^{2}} - {{\beta }^{2}}{\text{si}}{{{\text{n}}}^{2}}\vartheta }}{2}} \right)}}^{S}}} \right| \ll 1\) at small angles \(\rho ,\beta \ll 1\), as shown in Fig. 1с for the case \(\beta = 0\). It is important that only the matrix elements for the nearest neighbors i and \(i + \delta \) are needed to calculate the exchange contribution to \({{U}_{{{\text{sp}}}}}\) (see Eq. (3)). To illustrate, we present two matrix elements for relevant correlation functions \(\langle S_{i}^{x}S_{{i + \delta }}^{x}\rangle ,\langle S_{i}^{x}S_{{i + \delta }}^{y}\rangle \):

$$\begin{array}{*{20}{c}} {\langle {{{\mathbf{n}}}_{i}},{{{\mathbf{n}}}_{{i + \delta }}}{\text{|}}n_{i}^{y}n_{{i + \delta }}^{y}{\text{|}}{{{\mathbf{n}}}_{i}},{{{\mathbf{n}}}_{{i + \delta }}}\rangle } \\ { = - {{S}^{2}}\sin ({{\vartheta }_{i}} - {{\rho }_{i}})\sin ({{\vartheta }_{{i + \delta }}} + {{\rho }_{{i + \delta }}})} \\ { \times \;\sin ({{\varphi }_{i}} + {{\beta }_{i}})\sin ({{\varphi }_{{i + \delta }}} - {{\beta }_{{i + \delta }}}),} \\ {\langle {{{\mathbf{n}}}_{i}},{{{\mathbf{n}}}_{{i + \delta }}}{\text{|}}n_{i}^{x}n_{{i + \delta }}^{y} + n_{i}^{y}n_{{i + \delta }}^{x}{\text{|}}{{{\mathbf{n}}}_{i}},{{{\mathbf{n}}}_{{i + \delta }}}\rangle } \\ { = - {{S}^{2}}\sin ({{\vartheta }_{i}} - {{\rho }_{i}})\sin ({{\vartheta }_{{i + \delta }}} + {{\rho }_{{i + \delta }}})\sin 2{{\varphi }_{i}}.} \end{array}$$
(5)

4 LAGRANGIAN

The spin dynamics of the considered antiferromagnet is described by the Lagrangian (for clarity, it is given in the site representation)

$$\mathcal{L}({{\varphi }_{1}},{{\vartheta }_{1}},\rho ,\beta ) = {{T}_{{{\text{WZ}}}}} - U,$$
(6)

where

$$\begin{gathered} {{T}_{{{\text{WZ}}}}} = - S\hbar \sum\limits_{i = 1}^N \sin {{\vartheta }_{0}}({{\rho }_{i}}{{{\dot {\varphi }}}_{{1i}}} + {{\beta }_{i}}{{{\dot {\vartheta }}}_{{1i}}}) = \frac{{ - S\hbar \sin {{\vartheta }_{0}}}}{{2N}} \\ \times \,\sum\limits_{\mathbf{k}} ({{\rho }_{{\mathbf{k}}}}{{{\dot {\varphi }}}_{{1, - {\mathbf{k}}}}} + {{\beta }_{{\mathbf{k}}}}{{{\dot {\vartheta }}}_{{1, - {\mathbf{k}}}}} + {{\rho }_{{ - {\mathbf{k}}}}}{{{\dot {\varphi }}}_{{1,{\mathbf{k}}}}} + {{\beta }_{{ - {\mathbf{k}}}}}{{{\dot {\vartheta }}}_{{1,{\mathbf{k}}}}}), \\ \end{gathered} $$
(7)
$$U = {{U}_{{{\text{ex}}}}} + {{U}_{{\text{A}}}} + \Delta U_{{{\text{ex}}}}^{{{\text{sp}}}}.$$
(8)

The kinetic energy of the spin system (Wess–Zumino factor) is defined in quantum mechanics as a Berry phase: \({{\gamma }_{{\text{B}}}} = \sum\nolimits_{i = 1}^2 {{{{\mathbf{A}}}_{i}}{{{{\mathbf{\dot {n}}}}}_{i}}} \), where Ai is the Berry connection; although the Berry connection is not gauge invariant, the Berry phase itself is gauge invariant [21].Footnote 1 We use the well-known gauge in which \({{{\mathbf{A}}}_{i}}{{{\mathbf{\dot {n}}}}_{i}}\, = \,(1 - {\text{cos}}{{\vartheta }_{i}}){{\dot {\varphi }}_{i}}\). The substitution of \({{\vartheta }_{i}}\) and \({{\varphi }_{i}}\) from Eqs. (2) into this formula gives Eq. (7) for the kinetic energy of the antiferromagnet in the first approximation in \(\rho \) and \(\beta \).

If \(\Delta U_{{{\text{ex}}}}^{{{\text{sp}}}}\) is omitted in the Lagrangian (6), optimization in \(\rho \) and \(\beta \) specifies the known formula T = \(\frac{{{{\chi }_{ \bot }}}}{{2{{\gamma }^{2}}}}{{\left( {\frac{{d{\kern 1pt} {\mathbf{l}}}}{{dt}}} \right)}^{2}}\) for the kinetic energy of the antiferromagnet, where l = L/|L| is the unit antiferromagnetic vector and \({{\chi }_{ \bot }}\) is the perpendicular susceptibility of the antiferromagnet. This expression coincides with the known Andreev–Marchenko formula for the kinetic energy of the antiferromagnet [23].Footnote 2

The equations of the spin dynamics with such a term describing the kinetic energy and \(\Delta U_{{{\text{ex}}}}^{{{\text{sp}}}} = 0\) are very close to the equation of the mathematical pendulum. The inclusion of the term \(\Delta U_{{{\text{ex}}}}^{{{\text{sp}}}}\) introduces intriguing features of the famous problem of the mathematical pendulum with a vibrating suspension point, which has very diverse bifurcations, instabilities, and phase transitions.

Expressions for \({{U}_{{{\text{ex}}}}}\), \({{U}_{{\text{A}}}}\), and \(\Delta U_{{{\text{ex}}}}^{{{\text{sp}}}}\) are obtained from Eqs. (1) and (3) by averaging over the ground state of the Hamiltonian taking into account that the coherent spin states are nonorthogonal, as mentioned above. The second of two pairs of conjugate variables (\(\rho \), \(\varphi \)) and (\(\vartheta \), \(\beta \)) is omitted below because it does not affect the Cotton–Mouton effect and does not contribute to birefringence at \({\mathbf{E}} = ({{E}_{x}},{{E}_{y}},0)\), as will be shown below.

Thus, substituting ni = (sinϑicosφi, \({\text{sin}}{{\vartheta }_{i}}{\text{sin}}{{\varphi }_{i}}\), cosϑi) with \(\vartheta = {{\vartheta }_{0}} + {{\vartheta }_{1}}\) and \(\varphi = {{\varphi }_{0}} + {{\varphi }_{1}}\), where \({{\vartheta }_{1}},{{\varphi }_{1}} \ll 1\), into Eqs. (1) and (3), retaining only terms up to the second order in \({{\vartheta }_{1}}\), \({{\varphi }_{1}}\), \(\rho \), and \(\beta \), omitting harmonic terms as usual in the Raman process with the frequencies \({{\omega }_{1}} + {{\omega }_{2}}\), and taking the standard Fourier transform from the site representation to the quasimomentum one, we obtain

$$\begin{gathered} {{U}_{{{\text{ex}}}}} = J{{S}^{2}}\frac{Z}{{2N}} \\ \times \,\sum\limits_{\mathbf{k}} \text{[}(1 + {{\gamma }_{{\mathbf{k}}}}){{\rho }_{{\mathbf{k}}}}{{\rho }_{{ - {\mathbf{k}}}}} + {{\sin }^{2}}{{\vartheta }_{0}}(1 - {{\gamma }_{{\mathbf{k}}}}){{\varphi }_{{1,{\mathbf{k}}}}}{{\varphi }_{{1, - {\mathbf{k}}}}}], \\ \end{gathered} $$
(9)
$${{U}_{{\text{A}}}} = \frac{1}{N}\sum\limits_{\mathbf{k}} {\text{|}}K{\text{|}}( - {\kern 1pt} \cos 2{{\vartheta }_{0}}{{\rho }_{{\mathbf{k}}}}{{\rho }_{{ - {\mathbf{k}}}}} + {\text{si}}{{{\text{n}}}^{4}}{{\vartheta }_{0}}{{\varphi }_{{1,{\mathbf{k}}}}}{{\varphi }_{{1, - {\mathbf{k}}}}}),$$
(10)
$$\Delta U_{{{\text{ex}}}}^{{{\text{sp}}}} = {{S}^{2}}E_{0}^{2}\frac{Z}{{2N}}\sum\limits_{\mathbf{k}} \left\{ {\left[ {{{w}_{{\mathbf{k}}}}{{\rho }_{{\mathbf{k}}}}{{\rho }_{{ - {\mathbf{k}}}}}_{{_{{_{{_{{_{{_{{_{{}}}}}}}}}}}}}}} \right.} \right.$$
(11)
$$ + \,{\text{si}}{{{\text{n}}}^{2}}{{\vartheta }_{0}}\frac{{{{a}_{1}} + {{a}_{2}}}}{2}(1 - {{\gamma }_{{\mathbf{k}}}}){{\varphi }_{{1,{\mathbf{k}}}}}{{\varphi }_{{1, - {\mathbf{k}}}}}$$
$$\left. { \pm \frac{{{{a}_{1}} - {{a}_{2}}}}{2}{\text{si}}{{{\text{n}}}^{2}}{{\vartheta }_{0}}(1 - {{\gamma }_{{\mathbf{k}}}}){\text{cos}}2\alpha \frac{{{{\varphi }_{{1,{\mathbf{k}}}}} + {{\varphi }_{{1, - {\mathbf{k}}}}}}}{2}} \right]$$
$$\left. { \times \;f(t)\cos ({{\omega }_{{\text{p}}}}t + \zeta ) \mp \frac{{{{a}_{3}}}}{2}{{\rho }_{{\mathbf{k}}}}{{\rho }_{{ - {\mathbf{k}}}}}\sin 2\alpha f(t)\sin ({{\omega }_{{\text{p}}}}t + \zeta )} \right\},$$

where Z  is the number of nearest neighbors, N is the number of magnetic sites in the sample, wk = \(\frac{{{{a}_{1}}\, + \,{{a}_{2}}}}{2}{\text{si}}{{{\text{n}}}^{2}}{{\vartheta }_{0}}\) + a2cos2ϑ0 + \({{\gamma }_{{\mathbf{k}}}}\left( {\frac{{{{a}_{1}}\, + \,{{a}_{2}}}}{2}{\text{co}}{{{\text{s}}}^{2}}{{\vartheta }_{0}}\, + \,{{a}_{2}}{\text{si}}{{{\text{n}}}^{2}}{{\vartheta }_{0}}} \right)\), \(f(t)\) is the function describing the duration and shape of the laser pulse, \(\zeta \) is the phase depending on the time onset, and \({{\gamma }_{{\mathbf{k}}}} = \frac{1}{Z}\sum\limits_{{\mathbf{\delta }} = 1}^6 {{e}^{{i{\mathbf{k\delta }}}}}\).

5 EULER–LAGRANGE EQUATIONS

To describe the spin dynamics of the two-sublattice antiferromagnet, we use the Euler–Lagrange equation for the Lagrangian (6) with the terms given by Eqs. (7), (9), and (10), which have the form

$$\begin{array}{*{20}{c}} {{{{\dot {\rho }}}_{{\mathbf{k}}}} - {{\omega }_{\varphi }}[1 + {{h}_{\varphi }}\psi (t)]{{\varphi }_{{1,{\mathbf{k}}}}} = \pm T_{{\mathbf{k}}}^{\varphi },} \\ {{{{\dot {\varphi }}}_{{1,{\mathbf{k}}}}} + {{\omega }_{\rho }}[1 + {{h}_{\rho }}\psi (t)]{{\rho }_{{\mathbf{k}}}} = 0,} \end{array}$$
(12)

where

$${{\omega }_{\varphi }} = \frac{{JSZ(1 - {{\gamma }_{{\mathbf{k}}}})}}{{2\hbar \sin {{\vartheta }_{0}}}} + \frac{{4{\text{|}}K{\text{|}}}}{{3S\hbar \sin {{\vartheta }_{0}}}},$$
$${{\omega }_{\rho }} = \frac{{JSZ(1 + {{\gamma }_{{\mathbf{k}}}})}}{{2\hbar \sin {{\vartheta }_{0}}}} - \,{\text{|}}K{\text{|}}\frac{{\cos 2{{\vartheta }_{0}}}}{{S\hbar \sin {{\vartheta }_{0}}}},$$
$$T_{{\mathbf{k}}}^{\varphi } = \frac{{E_{0}^{2}S\sin {{\vartheta }_{0}}}}{\hbar }\frac{{{{a}_{1}} - {{a}_{2}}}}{2}\cos (2\alpha )(1 - {{\gamma }_{{\mathbf{k}}}})\psi (t),$$
$${{h}_{\varphi }} = \frac{{E_{0}^{2}{\text{si}}{{{\text{n}}}^{2}}{{\vartheta }_{0}}}}{J}\frac{{{{a}_{1}} + {{a}_{2}}}}{2},$$
$${{h}_{\rho }} = \frac{{E_{0}^{2}{{w}_{{\mathbf{k}}}}}}{{J(1 + {{\gamma }_{{\mathbf{k}}}})}},$$

and \(\psi (t) = f(t)\cos ({{\omega }_{{\text{p}}}}t + \zeta )\).

Equations (12) describe the dynamics of a magnon with the frequency \({{\omega }_{{{\mathbf{k}},{\text{M}}}}} = \sqrt {{{\omega }_{\varphi }}{{\omega }_{\rho }}} \), which is excited simultaneously by the external torque (right-hand side of Eq. (12)), which is due to the inverse Cotton–Mouton effect, and by the internal parametric torque. The latter is manifested mainly in the instability region, i.e., in the parametric resonance region. The presence of the parametric instability significantly enhances the response of the system to the inverse Cotton–Mouton effect (Fig. 2a). According to the Floquet theory, the variables \({{\rho }_{{\mathbf{k}}}}\) and \({{\varphi }_{{1,{\mathbf{k}}}}}\) are represented in the form

$$\left( {\begin{array}{*{20}{c}} {{{\rho }_{{\mathbf{k}}}}} \\ {{{\varphi }_{{1,{\mathbf{k}}}}}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{{c}_{{11}}}} \\ {{{c}_{{21}}}} \end{array}} \right){{e}^{{i{{\omega }_{{\text{p}}}}t/2}}} + \left( {\begin{array}{*{20}{c}} {{{c}_{{12}}}} \\ {{{c}_{{22}}}} \end{array}} \right){{e}^{{ - i{{\omega }_{{\text{p}}}}t/2}}},$$
(13)

where \({{c}_{{i1}}}(t)\) and \({{c}_{{i2}}}(t)\) are the real periodic functions of the time with the pump period Tp = 2π/ωp; it is obvious that \(c_{{ij,{\mathbf{k}}}}^{*} = {{c}_{{ij, - {\mathbf{k}}}}}\) and the frequency ω is an eigenvalue of the uniform system of Eqs. (12). The functions \({{c}_{{i1}}}(t)\) and \({{c}_{{i2}}}(t)\) are similar to Bloch functions in the theory of crystals and the studied system of Eqs. (12) and similar parametric systems are recently called “time crystals” [25]. In this case, since we are interested only in the threshold of the appearance of the instability near the main parametric resonance, we suggest \({\text{|}}\omega - {{\omega }_{{\text{p}}}}{\text{/}}2{\text{|}} \ll 1\).

Fig. 2.
figure 2

(Color online) (a, b) Amplitude of the induced ellipticity obtained by solving Eq. (17) (a) versus the pump frequency ωp at pump intensities h = (from bottom to top) 0.2, 0.4, 0.6, and 0.8 and (b) versus the pump intensity h at the pump frequency ωp = (blue line) ωM and (red line) 2ωM, where ωM is the magnetic dynamics frequency. The vertical dashed line marks the threshold pump intensity for the parametric excitation. (c) Color (ωp, ωM) map of the amplitude of the induced ellipticity in a logarithmic scale at the pump intensity above the threshold of the parametric instability.

Substituting Eq. (13) into Eqs. (12), omitting nonresonance terms of the form exp(ipt), as usual in the Floquet theory, assuming \(f(t) \equiv 1\) for fairly long pulses, and performing simple algebra, we reduce the eigenvalue problem to the solution of four linear differential equations with constant coefficients. Then, the eigenvalue problem for the resulting differential equations is reduced to the algebraic equation expressing the zero determinant of the 4 × 4 matrix of coefficients of Eqs. (12). The solution of this equation is lengthy and rather unclear. It is significantly simplified in the case \({{h}_{\varphi }} = 0\), holding the main features of the general solution and can be represented in the known form [11]

$${{\rho }_{{\mathbf{k}}}} = - {{\nu }_{{\mathbf{k}}}} + \sqrt { - {{{(\omega - {{\omega }_{{\mathbf{k}}}})}}^{2}} + \frac{{h_{\rho }^{2}\omega _{{\text{p}}}^{2}}}{4}} ,$$
(14)

where ρk is the imaginary part of the eigenfrequency ω near the threshold of the parametric instability and determines the growth rate of the amplitude of generated magnons and νk is the damping frequency of magnons. This formula immediately yields the known expression for the threshold of the parametric instability at the pump frequency ωp ≈ 2ωk [8]; in this case, the intensity of the pump field should be \({\text{|}}{{h}_{\rho }}{{\omega }_{{\mathbf{k}}}}{\text{|}} = {{\nu }_{{\mathbf{k}}}}\). This threshold of the parametric instability is apparently the lowest. Other unstable modes with different thresholds and frequencies can also be excited in the system under consideration (see, e.g., Fig. 3c).

Fig. 3.
figure 3

(Color online) (a) Waveforms of pump pulses \(\psi (t)\) with various indicated durations Δτ calculated by Eq. (18). (b, c) Time dependences of the ellipticity induced by laser pump pulses with the polarization \({\mathbf{E}}\parallel x\) (shown in panel (a) in the corresponding colors) obtained by solving Eq. (17) at ωp = 2ωM and (b) zero and (c) nonzero h values with indicated scaling factors. (d) Time dependences of the ellipticity induced by laser pump pulses with the duration Δτ = 6.8 at h = 0 and h above the threshold of the parametric instability with indicated scaling factors. (e) Amplitude of the induced ellipticity versus the duration of laser pump pulses Δτ at h = 0 and h values below and above the threshold of the parametric instability.

The presented formula for the frequency threshold of the parametric instability specifies the region (surface) in the k space that contains coherent pairs of k and –k magnons interacting with each other through the pump field. Such pairs of magnons are similar to Cooper pairs in superconductors [11], but unlike the latter, they are formed only under pumping although they can exist during a certain time after the end of pumping. Some analogy can also be indicated with the production of pairs of k and –k electrons in nonlinear quantum optics. Such a coherent paired state of magnons is quantum-mechanically entangled state and opens interesting experimental prospects.

6 COTTON–MOUTON EFFECT

We recall that the Cotton–Mouton effect (or the Voigt effect or the magnetic linear birefringence) is a magneto-optical effect according to which linearly polarized light with the plane of polarization tilted to the magnetic order parameter (magnetization or antiferromagnetic vector), passing through a medium becomes elliptically polarized; this effect is quadratic in the magnetic order parameter [26, 27]. The inverse Cotton–Mouton effect is a change in the magnetic order parameter under the action of linearly polarized light [1].

Pump-induced generation of spin waves, which is described by the functions ρk and \(\frac{{{{\varphi }_{{1,{\mathbf{k}}}}} + {{\varphi }_{{1, - {\mathbf{k}}}}}}}{2}\), can be detected by probe pulses through the magnetic linear birefringence (i.e., the direct Cotton–Mouton effect). Differentiating \(\Delta U_{{{\text{ex}}}}^{{{\text{sp}}}}\) (11) twice with respect to Ex and twice with respect to \({{E}_{y}}\), we obtain

$$\begin{gathered} {{\varepsilon }_{{xx}}} - {{\varepsilon }_{{yy}}} = \mp J{{S}^{2}}Z{{\sin }^{2}}{{\vartheta }_{0}}\frac{{{{a}_{1}} - {{a}_{2}}}}{2}\cos (2\alpha ) \\ \times \;\sum\limits_{\mathbf{k}} (1 - {{\gamma }_{{\mathbf{k}}}})\frac{{{{\varphi }_{{1,{\mathbf{k}}}}} + {{\varphi }_{{1, - {\mathbf{k}}}}}}}{2}. \\ \end{gathered} $$
(15)

The substitution of Eq. (15) into the Fresnel equation at k = (0, 0, 1) yields the formula for the birefringence effect [26]

$$\begin{gathered} \frac{{\Delta n}}{{{{n}_{0}}}} = {{n}_{1}} - {{n}_{2}} \approx \pm \frac{2}{3}\frac{\varkappa }{{n_{0}^{2}}}\cos (2\alpha ) \\ \times \;\sum\limits_{\mathbf{k}} (1 - {{\gamma }_{{\mathbf{k}}}})\frac{{{{\varphi }_{{1,{\mathbf{k}}}}} + {{\varphi }_{{1, - {\mathbf{k}}}}}}}{2}, \\ \end{gathered} $$
(16)

where ϰ is the birefringence index, which is related to the permittivity as \({{\varepsilon }_{{ii}}} = {{\varepsilon }_{0}} + \varkappa l_{i}^{2}\) due to symmetry and is measured by standard methods; in this model, ϰ = JS2ZN(a1a2)/2; n0 is the refractive index; and φ1,k(t) = \(\frac{{{{{\dot {\varrho }}}_{{0,{\mathbf{k}}}}} \mp 2T_{{\mathbf{k}}}^{\varphi }}}{{{{\omega }_{\varphi }}}}\) is a solution of Eqs. (12), where ‒ and + signs refer to domains with φ = π/2 and 3π/2. It is noteworthy that the ellipticity of light, which is proportional to the birefringence \(\Delta n\) and the thickness of the medium, rather than the birefringence is directly measured in the magneto-optical Cotton–Mouton effect [28]. Thus, we discuss below the ellipticity induced by the pump pulse in the antiferromagnet with unit thickness.

7 SIMULATION OF SPIN DYNAMICS

The solution of the system of differential equations (12) describing the spin dynamics of the antiferromagnet induced by the pump pulse is very cumbersome in the general case at hφ ≠ 0 and hρ ≠ 0. Figure 1d presents the color (hφ, hρ) map of the difference between the amplitude of oscillations of the ellipticity induced by pump pulses parametrically (ωp = 2ωM) and resonantly (ωp = ωM). Regions beyond the black area correspond to the coefficients hφ and hρ at which the parametric instability is developed and the amplitude of parametrically induced oscillations is larger than that of resonantly induced oscillations and their d-ifference is positive. The boundaries of the black region correspond to the threshold of the parametric instability.

To reveal the main features of parametric resonance and to enhance the Cotton–Mouton effect due to the parametric instability, we simulated only a particular case where \({{a}_{1}} = - {{a}_{2}}\) (hφ = 0 and h = \({{h}_{\varrho }}\)). In this case, the system of Eqs. (12) can be reduced to a second-order differential equation with constant coefficients, which was studied in detail in many theoretical researches of parametric resonance [8, 11, 12, 29]. We introduce a new variable \(\varrho \) = ρk + ρk, which describes the coherent paired state with zero quasimomentum mentioned above. Then, the system of Eqs. (12), where the variable \({{\varphi }_{{0,{\mathbf{k}}}}}\) is excluded, can be transformed to the second-order equation for \(\varrho \)

$$\ddot {\varrho } + \frac{2}{\tau }\dot {\varrho } + {{\omega }_{\varphi }}{{\omega }_{\rho }}[1 + h\psi (t)]\varrho = d\dot {\psi }(t)\cos 2\alpha ,$$
(17)

where the damping constant τ, the spin dynamics eigenfrequency \(\sqrt {{{\omega }_{\varphi }}{{\omega }_{\rho }}} \), the parameter h proportional to the laser pump intensity, the parameter d determining the magnitude Cotton–Mouton effect and the pump intensity, and the time t in units of the period of pump pulse oscillations are the dimensionless parameters used for the simulation. The function \(\psi (t)\) describes the duration and shape of the pump pulse according to the expression

$$\psi (t) = \left[ {\frac{1}{{{{e}^{{n(t - {{t}_{2}})}}} + 1}} - \frac{1}{{{{e}^{{n(t - {{t}_{1}})}}} + 1}}} \right]\cos ({{\omega }_{{\text{p}}}}t),$$
(18)

where \({{t}_{1}}\) and \({{t}_{2}}\) are the start and end times of the pump pulse (\(\Delta \tau = {{t}_{2}} - {{t}_{1}}\)) and n is the steepness of the pulse edges. It is worth noting that the time derivative of the canting angle of the sublattices \(\dot {\rho }\), which is a solution of Eq. (17), determines the time dependence for the induced ellipticity according to Eq. (16).

The 45 fs pump pulse with a central energy of 1.03 eV was used in a laser-induced magnetic dynamic experiment in the Heisenberg antiferromagnet RbMnF3 [30]. The parameters in Eq. (18) were chosen such that the model pump pulse includes about ten cycles (\(\Delta \tau = {{t}_{2}} - {{t}_{1}} = 10\)) as in the experiment, as shown by the green line in Figs. 4a and 4b. Solving Eq. (17) with the pump frequency ωp = 2ωM and the pump intensity h above the threshold of the parametric instability, we obtained time dependences of the induced ellipticity, which were calculated using Eq. (16) and include oscillations at the frequency ωM, as shown in Fig. 4a. In addition, it is seen that the phase of the induced ellipticity for two polarizations of the pump pulse E along the x (α = 0°, blue line) and y (α = 90°, red line) axes differs by π. Furthermore, as seen in Figs. 4c and 4d, at the rotation of the plane of polarization of the pump pulse (α = 0°–360°), the amplitude of the induced ellipticity behaves as cos2α, where the sign change corresponds to a change in the phase of oscillations by π, which corresponds to the experimental result obtained in [30]. In the case of Eq. (17) with zero right-hand side (d = 0), pump pulses with the same parameters will excite oscillations of the induced ellipticity at the same frequency ωM (see Fig. 4b), but their amplitude is independent of the angle of polarization, as shown in Figs. 4e and 4f.

Fig. 4.
figure 4

(Color online) (a, b) Time dependences of the ellipticity induced by laser pulses with the polarizations \({\mathbf{E}}\) along the x (α = 0°, blue line) and y (α = 90°, red line) axes due to parametric resonance and the inverse Cotton–Mouton effect in the two-sublattice antiferromagnet calculated with (a) nonzero (d ≠ 0) and (b) zero (d = 0) right-hand side of Eq. (17). (c–f) Amplitude of oscillations of the induced ellipticity versus the polarization angle α of the electric field \({\mathbf{E}}\) of the laser pulse calculated with (с, d) nonzero (d ≠ 0) and (e, f) zero (d = 0) right-hand side of Eq. (17). The sign change in panel (c) under the variation of the angle α corresponds to change in the phase of oscillations of the ellipticity by π. Polar diagrams (d, f) present the amplitude of oscillations of the induced ellipticity.

The analysis of solutions of Eq. (17) with a fixed spin dynamics eigenfrequency ωM at the varying pump frequency ωp and several intensity values h showed that the ellipticity at coinciding frequencies ωp = ωM is induced at any positive intensities h, as shown in Fig. 2a. On the contrary, the amplitude of oscillations of the induced ellipticity under the parametric excitation at the pump frequency ωp = 2ωM has a pronounced threshold in the pump intensity h. The spin dynamics at h > 0.57 is more efficiently excited under parametric pumping at the frequency ωp = 2ωM than under resonant pumping at the frequency ωp = ωM, as seen in Figs. 2a and 2b. The threshold of the parametric instability generally depends on the duration of the pump pulse: the longer the duration, the lower the threshold. In addition, the threshold depends on the shape of pump pulses, which requires further investigation. The amplitude of induced oscillations of the ellipticity under the parametric excitation at the frequency ωp = 2ωM exponentially depends on the pump intensity h above the threshold, as shown in Fig. 2b. It is noteworthy that the amplitude of oscillations of the induced ellipticity in the experiment linearly depends on the pump intensity [30]. This occurs likely because nonlinear effects limiting parametric excitation were ignored in our consideration. The amplitude of oscillations of the induced ellipticity linearly depends on the parameter d on the right-hand side of Eq. (17), which indicates that this coefficient determines the magnitude of the Cotton–Mouton effect in this pr-ocess.

Figure 2с presents the color (ωp, ωM) map of the amplitude of oscillations of the induced ellipticity in the logarithmic scale, where ωp and ωM are the pump and spin dynamics frequencies, respectively. The minimum and maximum amplitudes are marked in green and brown, respectively. It is seen that pumping at the frequency ωp excites the spin dynamics not only at the half frequency ωM = ωp/2 but also at the resonant frequency ωM = ωp and at other frequencies.

Further, we considered solutions of Eq. (17) for cases where h = 0, i.e., without the parametric instability, and h ≠ 0 at various durations of pump pulses Δτ with the time profile \(\psi (t)\) and the frequency ωp = 2ωM, as shown in Fig. 3. Below, pump pulses are assumed to be polarized along the x axis, unless indicated otherwise. When h = 0, oscillations of the ellipticity are induced at any duration Δτ of pump pulses, and their amplitude is independent of the duration, as shown in Figs. 3b and 3e. On the contrary, according to Eq. (17), the amplitude of oscillations of the induced ellipticity at the parameter h ≠ 0 increases exponentially with the duration of pulses Δτ, as seen in Figs. 3c and 3e. Time dependences of the induced ellipticity are plotted with the scaling factors indicated in Figs. 3c and 3d. The sharp increase in the amplitude of oscillations of the induced ellipticity under the parametric excitation is observed immediately after the pump pulse, i.e., in the time Δτ, as seen in Fig. 3d. As mentioned above, the threshold of the parametric instability depends on the duration of the pulse Δτ and the fixed parameter h for pulses with Δτ < 6 and Δτ > 6 is below and above the threshold, respectively, as shown in Fig. 3e. As a result, the giant parametric amplification of the inverse Cotton–Mouton effect with increasing duration Δτ of the pump pulse is observed in the antiferromagnet.

The mechanism of the parametric instability considered in this work can be implemented due to the following processes of the interaction of light with spins. First, the electric field of the laser pulse \({\mathbf{E}}\), modulating the exchange interaction between spins \(J{{{\mathbf{S}}}_{i}}{{{\mathbf{S}}}_{j}}\), creates a condition for the parametric instability. This can occur through, e.g., the inverse magnetorefractive effect [31, 32]. The inverse Cotton–Mouton effect produces the external torque acting on spins, results in the ρ canting of neighbor spins, and ensures the “azimuthal” displacement of spins and, thereby, the antiferromagnetic vector \({\mathbf{L}}\) from the equilibrium position. In this case, the parametric instability is developed even under zero initial conditions. In turn, the dependence of generation of magnons on the polarization of the incident light is specified by the polarization dependence of the inverse Cotton–Mouton effect and is proportional to cos2α.

8 CONCLUSIONS

To summarize, the effect of the canting of magnetic sublattices under the action of the short laser pump pulse on the parametric excitation of the inverse Cotton–Mouton effect in a Heisenberg antiferromagnet has been theoretically studied. It has been shown that the magnitude of the inverse Cotton–Mouton effect increases very strongly near parametric resonance, i.e., in the parametric instability region. It has been found that the pump-induced ellipticity is proportional to the time derivative of the canting angle of the magnetic sublattices. It has been established that the canting angle of the magnetic sublattices is described by a Hill differential equation of parametric resonance with a driving force. It has been revealed that the parametrically induced ellipticity in the presence of the driving force in the differential equation depends on the angle between the polarization of pump pulses and the principal crystallographic axes as cos2α, and the rotation of the plane of polarization by 90° changes the phase of oscillations by π. When the driving force is zero, the amplitude of oscillations is independent of the angle of polarization. It has been found that the parametric excitation of the induced ellipticity occurs when the pump intensity exceeds the threshold value. It has been shown that the amplitude of the induced ellipticity under the parametric excitation exponentially depends on the duration and intensity of the pump pulse.