Light-induced anisotropic skyrmion and stripe phases in a Rashba ferromagnet

TL;DR

This paper explores how off-resonant light can be used to control magnetic phases in Rashba ferromagnets, enabling stabilization and manipulation of skyrmions and stripe phases through induced anisotropy of the Dzyaloshinskii-Moriya interaction.

Contribution

It introduces a theoretical and numerical framework showing how linearly polarized off-resonant light induces anisotropic DMI, stabilizing novel magnetic phases and controlling skyrmion properties.

Findings

Off-resonant light induces strong anisotropy in DMI.

Stabilization of noncollinear magnetic phases.

Control over skyrmion stability, shape, and size.

Abstract

An external off-resonant pumping is proposed as a tool to control the Dzyaloshinskii-Moriya interaction (DMI) in ferromagnetic layers with strong spin-orbit coupling. Combining theoretical analysis with numerical simulations for an $s$-$d$-like model we demonstrate that linearly polarized off-resonant light may help stabilizing novel noncollinear magnetic phases by inducing a strong anisotropy of the DMI. We also investigate how with the application of electromagnetic pumping one can control the stability, shape and size of individual skyrmions to make them suitable for potential applications.

Equations

$$H=p^2/2m+\alpha \mathbf{\sigma }\cdot \left(\mathbf{p}\times \hat{\mathbf{z}}\right)+\Delta \mathbf{m}\cdot \mathbf{\sigma },$$

$$U_t=2^{-1/2}({\sigma }_x+{\sigma }_z)e^{-it\left[u_t+{\sigma }_{z}e\alpha A_0\sin (\omega t)/\omega t\right]/\hslash },$$

$$u_t=\frac{e^2{A}_0^2}{4m}\left(1+\frac{\sin (2\omega t)}{2\omega t}\right)+\frac{e{p}_y{A}_0}{m}\frac{\sin (\omega t)}{\omega t}.$$

$$H=\frac{p^2}{2m}+\left(\alpha p_y+\Delta m_x\right){\sigma }_z+\sum _{n=-\infty }^{\infty }{h}_n{e}^{in\omega t},$$

$$h_n=\left[\Delta m_z{\sigma }_x+\left(\alpha p_x-\Delta m_y\right){\sigma }_y\right]J_n(\gamma ),$$

$$H_{\mathrm{eff}}=p^2/2m+\alpha \left[p_y{\sigma }_x-p_x{\sigma }_y{J}_0(\gamma )\right]+V,$$

$$V=\Delta \{m_x{\sigma }_x+\left[m_y{\sigma }_y+m_z{\sigma }_z\right]J_0(\gamma )\},$$

$$G_0(i\omega ,\mathbf{k})=\frac{1}{2}\sum _{s=\pm }\frac{{\Lambda }_s({\theta }_{\mathbf{k}})}{i\hslash \omega -{\epsilon }_{\mathbf{k}}^s},$$

$${\Lambda }_{\pm }(\theta )=1\pm \left[{\sigma }_x\sin \theta -{\sigma }_y{J}_0(\gamma )\cos \theta \right]/g(\theta ).$$

$$\delta \mathcal{S}[\mathbf{m}]=\frac{1}{\beta }\sum _{i\omega ,\mathbf{k}}\sum _{ij}{\Pi }_{ij}(i\omega ,\mathbf{k})m_i(\mathbf{k})m_j(-\mathbf{k}),$$

$${\Pi }_{ij}(i\omega ,\mathbf{k})=\frac{1}{4}{\Delta }_i{\Delta }_j\sum _{\mathbf{q}}\sum _{s,s^{\prime }=\pm }\frac{f({\epsilon }_{\mathbf{k}+\mathbf{q}}^s)-f({\epsilon }_{\mathbf{q}}^{s^{\prime }})}{{\epsilon }_{\mathbf{k}+\mathbf{q}}^s-{\epsilon }_{\mathbf{q}}^{s^{\prime }}-i\hslash \omega }$$

$$\times \mathrm{Tr}\left[{\Lambda }_s({\theta }_{\mathbf{k}+\mathbf{q}}){\sigma }_i{\Lambda }_{s^{\prime }}({\theta }_{\mathbf{q}}){\sigma }_j\right],$$

$$\delta \mathcal{S}[\mathbf{m}]=\int d^{2}r\left(D_x{L}_{xz}^{(x)}+D_y{L}_{yz}^{(y)}\right),$$

$$D_x=\frac{{\Delta }^2}{\pi \alpha \hslash }\frac{\left|J_0(\gamma )\right|}{1+\left|J_0(\gamma )\right|},D_y=\frac{{\Delta }^2}{\pi \alpha \hslash }\frac{J_0^2(\gamma )}{1+\left|J_0(\gamma )\right|},$$

$$E=-J\sum _{\mathbf{r}}{\mathbf{S}}_{\mathbf{r}}\cdot \left[{\mathbf{S}}_{\mathbf{r}+\hat{\mathbf{x}}}+{\mathbf{S}}_{\mathbf{r}+\hat{\mathbf{y}}}\right]-H_{\text{ext}}\sum _{\mathbf{r}}{\left({\mathbf{S}}_{\mathbf{r}}\right)}_z$$

$$-D_x\sum _{\mathbf{r}}{\left({\mathbf{S}}_{\mathbf{r}}\times {\mathbf{S}}_{\mathbf{r}+\hat{\mathbf{x}}}\right)}_y+D_y\sum _{\mathbf{r}}{\left({\mathbf{S}}_{\mathbf{r}}\times {\mathbf{S}}_{\mathbf{r}+\hat{\mathbf{y}}}\right)}_x,$$

$$H=\sum _{n=-\infty }^{\infty }{h}_n{e}^{in\omega t},\mathrm{with}{h}_n=\left[\frac{p^2}{2m}+\left(\alpha p_y+\Delta m_x\right){\sigma }_z\right]{\delta }_{n,0}+\left[\Delta m_z{\sigma }_x+\left(\alpha p_x-\Delta m_y\right){\sigma }_y\right]J_n(\gamma ).$$

$$H_{\mathrm{eff}}=H^{(0)}+H^{(1)}+H^{(2)}+\ldots=h_{0}+\frac{1}{\hbar\omega}\sum\limits_{m\neq 0}\frac{h_{m}h_{-m}}{m}+\frac{1}{(\hbar\omega)^{2}}\bigg{(}\sum\limits_{m,n\neq 0}\frac{h_{m}h_{n-m}h_{-n}}{mn}-\sum\limits_{m\neq 0}\frac{h_{m}h_{-m}h_{0}}{m^{2}}\bigg{)},$$

$$\displaystyle H_{\mathrm{eff}}=\frac{p^{2}}{2m}+\bigg{(}1-2\beta\sum\limits_{n\neq 0}\frac{J_{n}(\gamma)J_{-n}(\gamma)}{n^{2}}\bigg{)}\left(\alpha p_{y}+\Delta m_{x}\right)\sigma_{z}$$

$$\displaystyle+\bigg{(}J_{0}(\gamma)+\beta\sum_{\begin{subarray}{c}m\neq n\\ m,n\neq 0\end{subarray}}\frac{J_{m}(\gamma)J_{n-m}(\gamma)J_{-n}(\gamma)}{mn}\bigg{)}\left(\Delta m_{z}\sigma_{x}+[\alpha p_{x}-\Delta m_{y}]\sigma_{y}\right),$$

$$H_{\mathrm{eff}}=\frac{p^2}{2m}+({\alpha }_y{p}_y+{\Delta }_x{m}_x){\sigma }_x-({\alpha }_x{p}_x-{\Delta }_y{m}_y){\sigma }_y+{\Delta }_z{m}_z{\sigma }_z,$$

$$C=\Big{|}\sum_{\begin{subarray}{c}m\neq n\\ m,n\neq 0\end{subarray}}\frac{J_{m}(\gamma_{0})J_{n-m}(\gamma_{0})J_{-n}(\gamma_{0})}{mn}\Big{|}=6\cdot 10^{-6},$$

$$N_i(\mathbf{k})=\frac{1}{2}\int \frac{d^{2}q}{(2\pi {)}^2}\sum _{s,s^{\prime }=\pm }\frac{f({\epsilon }_{\mathbf{k}+\mathbf{q}}^s)-f({\epsilon }_{\mathbf{q}}^{s^{\prime }})}{{\epsilon }_{\mathbf{k}+\mathbf{q}}^s-{\epsilon }_{\mathbf{q}}^{s^{\prime }}}\left(s\frac{{\hat{\mathbf{e}}}_i\cdot {\mathbf{n}}_{\mathbf{k}+\mathbf{q}}}{g({\theta }_{\mathbf{k}+\mathbf{q}})}-s^{\prime }\frac{{\hat{\mathbf{e}}}_i\cdot {\mathbf{n}}_{\mathbf{q}}}{g({\theta }_{\mathbf{q}})}\right),$$

$${\Pi }_{xz}(\mathbf{k})=-{\Pi }_{zx}(\mathbf{k})=J_0(\gamma ){\Delta }_x{\Delta }_z{N}_x(\mathbf{k}),\mathrm{and}{\Pi }_{yz}(\mathbf{k})=-{\Pi }_{zy}(\mathbf{k})={\Delta }_y{\Delta }_z{N}_y(\mathbf{k}),$$

$$N_x^{(1)}(\mathbf{k})=-\frac{1}{\pi \alpha \hslash }\underset{0}{\overset{2\pi }{\int }}\frac{d\theta }{2\pi }\frac{k_x{\cos }^2\theta +k_y\sin \theta \cos \theta }{{\sin }^2\theta +J_0^2(\gamma ){\cos }^2\theta }=-\frac{k_x}{\pi \alpha \hslash }\frac{1}{|J_0(\gamma )|}\frac{1}{1+|J_0(\gamma )|},$$

$$N_y^{(1)}(\mathbf{k})=-\frac{1}{\pi \alpha \hslash }\underset{0}{\overset{2\pi }{\int }}\frac{d\theta }{2\pi }\frac{k_y{\sin }^2\theta +k_x\sin \theta \cos \theta }{{\sin }^2\theta +J_0^2(\gamma ){\cos }^2\theta }=-\frac{k_y}{\pi \alpha \hslash }\frac{1}{1+|J_0(\gamma )|}.$$

$$\mathcal{S}[\mathbf{m}]={\mathcal{S}}_0[\mathbf{m}]+\delta \mathcal{S}[\mathbf{m}]={\mathcal{S}}_0[\mathbf{m}]+\int d^{2}r\left(D_x{\Lambda }_{xz}^{(x)}+D_y{\Lambda }_{yz}^{(y)}\right),$$

$$(1+{\alpha }^2)\frac{d{\mathbf{S}}_{\mathbf{r}}}{dt}=-\gamma {\mathbf{S}}_{\mathbf{r}}\times {\mathbf{H}}_{eff}-\alpha \gamma {\mathbf{S}}_{\mathbf{r}}\times {\mathbf{S}}_{\mathbf{r}}\times {\mathbf{H}}_{eff},$$

$${\mathbf{H}}_{eff}=-\frac{1}{\hslash \gamma }\frac{\delta E}{\delta {\mathbf{S}}_{\mathbf{r}}}=H_{ext}\hat{\mathbf{z}}+J\sum _{\mathbf{d}=\pm \hat{\mathbf{x}},\hat{\mathbf{y}}}{\mathbf{S}}_{\mathbf{r}+\mathbf{d}}+D_x({\mathbf{S}}_{\mathbf{r}+\hat{\mathbf{x}}}-{\mathbf{S}}_{\mathbf{r}-\hat{\mathbf{x}}})\times \hat{\mathbf{y}}-D_y({\mathbf{S}}_{\mathbf{r}+\hat{\mathbf{y}}}-{\mathbf{S}}_{\mathbf{r}-\hat{\mathbf{y}}})\times \hat{\mathbf{x}}.$$

Full text

Low-dimensional magnetic structures provide an exciting playground for condensed matter physics and technology applications. Some of them, e. g., helical magnets Ishikawa1976 ; Moriya1982 ; Pfleiderer2004 ; Li2012 , are known to support topologically nontrivial magnetic textures Kishine2011 ; Togawa2012 ; Koumpouras2016 . Such noncollinear states emerge as a result of the competition between Heisenberg exchange, antisymmetric Dzyaloshinskii-Moriya exchange, and magnetocrystalline anisotropy yielding magnetic ground states that are far more intricate than those in homogeneous ferromagnets Bak1980 ; Maleyev2006 ; Belitz2006 ; Vedmedenko2007 ; Grigoriev2009 . Recent progress in the fabrication of magnetic materials motivated an interest in particlelike domains, such as magnetic Skyrmions Muhlbauer2009 ; Yu2010 ; Heinze2011 ; Yu2011 ; Kiselev2011 ; Tonomura2012 ; Nagaosa2013 ; Rybakov2013 ; Pereiro2014 ; Gayles2015 ; Keesman2015 ; Gungordu2016 ; Leonov2016b , that are typical for nonhomogeneous ferromagnets. Skyrmions and other topologically protected magnetic textures have been proposed as building blocks for logical operations and information storage in the rapidly advancing fields of magnon spintronics Chumak2015 and Skyrmionics Krause2016 ; Leonov2016a .

Spintronics, as a branch of applied science, is traced back to the pioneering work on giant magnetoresistance by Grünberg Saurenbach1988 ; Binasch1989 and Fert Baibich1988 . The subsequent discovery of spin-transfer torque Slonczewski1989 ; Berger1996 provoked the idea to exploit the spin degree of freedom rather than the charge for processing and transferring information. This idea is currently being extended to magnetic materials supporting localized magnetic excitations such as Skyrmions.

Ways to create and control Skyrmions and other noncollinear magnetic textures are essential for the practical implementation of this emerging technology. In this Letter, we investigate microscopically how the noncollinear magnetic domains in thin ferromagnet layers with strong spin-orbit coupling may be controlled by linearly polarized light. The effects predicted may be observed in thin films such as Co/Pt heterostructures subject to short light pulses Huisman2016 .

The exchange interaction alone may lead only to a collinear orientation of magnetic moments in a cubic crystal. Spatially inhomogeneous magnets are usually associated with a lack of lattice inversion symmetry. In his seminal work on noncentrosymmetric magnets Dzyaloshinskii1964 , Dzyaloshinskii identified the one-dimensional magnetic spiral states stabilized by the Dzyaloshinskii-Moriya interaction (DMI) Dzyaloshinskii1958 that favors the noncollinear orientation of neighboring spins. A nontrivial ground state arises in helical magnets as a consequence of the competition between the Heisenberg exchange and the DMI. In two-dimensional structures, this competition leads to a helical spin-spiral ground state configuration that becomes unstable in the presence of a magnetic field with the tendency towards the formation of Skyrmions Rosler2006 . Quite generally, Skyrmions correspond to the solitonlike solutions of the field equations of Dzyaloshinskii’s theory that destroy the homogeneity of magnetic order Bogdanov1989 . The existence of such localized states and the mechanism of their nucleation as mesoscopic objects are rather common for continuum systems described by the free energy functional with Lifshitz invariants Rosler2011 . The strength of DMI can be rigorously approached by the microscopic theory, as an indirect exchange interaction between two neighboring spins facilitated by itinerant (conduction) electrons Imamura2004 ; Kundu2015 ; Kikuchi2016 , as well as adopted from the first-principles simulations Nossa2012 ; Dmitrienko2014 ; Koretsune2015 . In a thin film, the DMI can be induced by the Rashba spin-orbit coupling.

An external electromagnetic field may strongly modify the properties of an electronic system providing an important tool for manipulating materials in a controllable fashion Mikhaylovskiy2015 ; Fujita2017a ; Fujita2017b . It has been recently shown that the effect of off-resonant electromagnetic radiation (with the frequency exceeding the bandwidth of the system) may be described by effective time-independent models with strongly renormalized parameters Oka2009 ; Kitagawa2011 ; Lindner2011 ; Wang2013 ; Cayssol2013 ; Yudin2016a ; Yudin2016b ; Stepanov2016 ; Yudin2017 ; Sato2017 .

In this Letter, we derive the effective $s$-$d$ exchange model of a Rashba ferromagnet in the regime of strong coupling to external radiation. We find that the DMI strength can be effectively controlled by the application of off-resonant pumping that opens up exciting opportunities for controlling the stability, size, and shape of individual metastable Skyrmions. Also, we show that the application of linearly polarized radiation induces anisotropy of the DMI that not only provides a finer control over the individual DMI strengths in two orthogonal directions but also leads to the appearance of novel anisotropic phases.

From the theory point of view, the effect of time-periodic fields may be described, with some reservations, by the so-called Floquet theory Grifoni1998 ; Kohler2005 . Periodicity of the driving field enables one to map the original time-dependent problem to the eigenvalue problem of Floquet states. Off-resonant pumping takes place if the frequency of the driving field is so high that electrons are not able to follow field oscillations. In this case, real absorption or emission of light quanta cannot happen due to restrictions imposed by energy conservation for radiation frequency exceeding the electron bandwidth. Still, such off-resonant radiation affects the system via virtual processes leading to a significant renormalization of the parameters of the initial Hamiltonian of an electron subsystem. Below, we restrict our attention to the effects of linearly polarized light, since it has a greater impact on the noncollinear magnetic textures. For the sake of microscopic treatment, we rely upon the Floquet-Magnus expansion and its generalizations that have been developed in Refs. Blanes2009 ; Goldman2014 ; Eckardt2015 ; Bukov2015 .

For microscopic analysis, we consider a weak two-dimensional ferromagnet that yields an $s$-$d$-like Rashba model for conduction electrons:

[TABLE]

where $\bm{m}$ is the unit ($|\bm{m}|=1$) local magnetization vector due to, e. g., localized $d$ electrons, $\Delta$ is the $s$-$d$-like exchange energy, $\bm{p}$ stands for the momentum operator for conduction electrons with an effective mass $m$, $\alpha$ is the Rashba spin-orbit interaction strength, $\hat{\bm{z}}$ is the unit vector in the direction perpendicular to the two-dimensional electron gas, and $\bm{\sigma}=\left(\sigma_{x},\sigma_{y},\sigma_{z}\right)$ denotes the vector of Pauli matrices. Models of the type of Eq. (1) were originally proposed to explain the physics of ferromagnetic metals beyond the Heisenberg exchange picture Vonsovsky1946 ; Zener1951a ; Zener1951b ; Zener1951c . This approach relies upon a formal distinction between a localized (classical) magnetic subsystem (e. g., $d$ or $f$ electrons, that are described by an $\bm{m}$ field which is governed by a classical Heisenberg model) an itinerant subsystem [e. g., $s$ electrons described by Eq. (1)] that are coupled to each other by means of exchange interaction.

An external ac electromagnetic field of frequency $\omega$ footnote is introduced in the model of Eq. (1) by means of the Peierls substitution $\bm{p}\to\bm{p}+e\bm{A}_{0}\cos\omega t$, where $\bm{A}_{0}=\bm{E}_{0}/\omega$ and $\bm{E}_{0}$ is the electric field component of the field. In what follows, we restrict our analysis to the case of linearly polarized light by choosing $\bm{E}_{0}=E_{0}\hat{\bm{y}}$, where $\hat{\bm{y}}$ is the in-plane unit vector (in the $y$ direction).

The Hamiltonian of Eq. (1) is conveniently rotated as $H\to U_{t}^{\dagger}HU_{t}$ with the time-dependent unitary operator

[TABLE]

The transformed Hamiltonian yields the matrix Floquet model of the form

[TABLE]

with the coefficients $h_{n}$ defined by

[TABLE]

where the parameter $\gamma=2e\alpha E_{0}/\hbar\omega^{2}$ describes the effective light-matter coupling and $J_{n}(\gamma)$ stands for the $n$th order Bessel function of the first kind.

The high-frequency expansion in the form of the Brillouin-Wigner perturbation theory recently developed for this class of problems Mikami2016 maps Eq. (3) onto an effective time-independent Hamiltonian:

[TABLE]

that is valid away from resonance frequencies. The effective model fails only in a tiny vicinity $\delta\gamma$ of the zeros of the Bessel function, $\delta\gamma\approx 10^{-5}\Delta^{2}/(\hbar\omega)^{2}$, which is well beyond our numerical resolution sup . The model is equivalent to that of Eq. (1) with anisotropic renormalization of coupling constants: Rashba spin-orbit interaction strength and $s$-$d$ exchange coupling. In what follows, we assume a weak ferromagnet and treat the exchange interaction term $V$ perturbatively.

Based on the symmetry analysis, Dzyaloshinskii discovered that the effective Ginzburg-Landau free energy functional may allow for terms linear in magnetization gradients provided the absence of lattice inversion symmetry Dzyaloshinskii1958 . Later, Moriya argued, on the basis of Anderson’s theory of superexchange, that the microscopic mechanism of spin-orbit coupling is responsible for such an interaction Moriya1960a ; Moriya1960b . The latter can also be thought as a coupling between an excited state of a magnetic ion and the ground state of the neighboring ion. Such a coupling can be derived microscopically from the correction to the bare action $\mathcal{S}_{0}[\bm{m}]$ (that collects all terms corresponding to magnetic subsystem) computed to the second order with respect to the perturbation $V$ sup .

To construct the perturbation theory, we take advantage of the bare Matsubara Green’s function for the Hamiltonian of Eq. (5a):

[TABLE]

where $g(\theta)=\left[\sin^{2}\theta+J_{0}^{2}(\gamma)\cos^{2}\theta\right]^{1/2}$ and the spectrum $\varepsilon_{\bm{k}}^{\pm}=\hbar^{2}k^{2}/2m\pm\alpha\hbar k\,g(\theta_{\bm{k}})$ acquires the dependence on the direction of the wave vector $\bm{k}=\bm{p}/\hbar=k(\cos\theta_{\bm{k}},\sin\theta_{\bm{k}})$ due to the linear polarization of light.

The second-order contribution to the effective action in the imaginary time representation is given by

[TABLE]

where $\beta$ stands for the inverse temperature while the indexes $i$ and $j$ denote the Cartesian vector components. The polarization operator is expressed as

[TABLE]

where $f(\varepsilon)$ is the Fermi-Dirac distribution function and $\Delta_{x}=\Delta$, $\Delta_{y}=\Delta_{z}=J_{0}(\gamma)\Delta$.

The straightforward expansion of $\Pi_{ij}(\omega=0,\bm{k})$ around $\bm{k}=0$ up to the terms linear in $\bm{k}$ yields a fully antisymmetric contribution to the effective action sup :

[TABLE]

where we introduce the DMI couplings

[TABLE]

and Lifshitz invariants $L_{ij}^{(l)}=m_{i}\partial_{l}m_{j}-m_{j}\partial_{l}m_{i}$.

In the absence of an electromagnetic field, i. e., for $\gamma=0$, we obtain an isotropic DMI with $D_{x}=D_{y}=\Delta^{2}/(2\pi\alpha\hbar)$. In the presence of linearly polarized light, the DMI coupling becomes essentially anisotropic as given by Eq. (10). We stress that the employed high-frequency expansion is legitimate only as far as there are no resonant transitions and the parameter $\gamma$ is away from zeros of the Bessel function $J_{0}(\gamma)$ sup .

To illustrate our results, we consider a classical two-dimensional Heisenberg exchange model on a square lattice, that is given by the total energy

[TABLE]

where $\bm{S}_{\bm{r}}$ is the spin on a lattice site $\bm{r}$, $H_{\textrm{ext}}$ is an external magnetic field (in energy units) perpendicular to the two-dimensional plane, $\hat{\bm{x}}$ and $\hat{\bm{y}}$ stand for the unit vectors in $x$ and $y$ direction, correspondingly, and the lattice constant is set to unity.

With the help of the numerical approach described in Ref. sup , we analyze the influence of the anisotropic DMI on the Skyrmion profile. The profile is obtained by relaxing a trial Skyrmion ansatz using the dynamical Landau-Lifshitz-Gilbert equation until the stationary state is reached. To avoid nonuniversal effects of boundary conditions, the numerical simulation is performed in a box of a large size that exceeds characteristic size of a Skyrmion by a large factor (only the central part of the box is shown in Fig. 1).

We find that the anisotropic renormalization of the DMI strength of Eq. (10) results in the anisotropic squeezing of a Skyrmion. The Skyrmion becomes elongated along the light polarization direction and develops an elliptic profile as shown in Fig. 1. In the model with $\gamma\lesssim 1$ and positive $D_{x}$ and $D_{y}$, the in-plane spin projections are directed towards the Skyrmion center. Such a configuration may be referred to as the inverted hedgehog Skyrmion that is distinguished from the hedgehog Skyrmion in which the in-plane spin projections point outwards. The Skyrmion type is, therefore, defined by the overall sign of the DMI coupling.

It is also worth stressing that the model (11) supports only Neél-type Skyrmions which are the Skyrmions with a radial orientation of spins. Such Neél-type Skyrmions were observed recently in GaV4S8 Kezsmarki2015 . Those should be contrasted with Bloch-type Skyrmions that are characterized by spin orientations perpendicular to the radial direction. The Bloch-type Skyrmions are thought to be characteristic for materials like FeGe Wilhelm2011 and MnSi Yu2015 .

It can be shown using the methodology of Ref. sup that individual Skyrmions are metastable only in certain areas of the parameter space as illustrated in Fig. 2(a). The metastability regions of individual Skyrmions in the model of Eq. (11) must be distinguished from the phases that characterize the absolute minimum of the energy functional. By extending our numerical analysis to search for a ground state sup , we obtain the phase diagram depicted in Fig. 2(b). The diagram consists of three phases: (i) the homogeneous ferromagnetic order phase denoted by points, (ii) the Skyrmion crystal phase (a crystal of elliptic Skyrmons) denoted by circles, and (iii) the stripe crystal phase (periodic stripes in the direction of light polarization) denoted by vertical lines.

Anisotropy induced by pumping distorts the symmetry of the Skyrmion crystal from the equilateral triangular at $\gamma=0$ to an isosceles triangular at nonzero $\gamma$ sup . The stripe crystal phase is analogous to the helical phase discussed, e. g., in Ref. Iwasaki-NComm-2013 , although, in contrast to the conventional helical phases, the stripe phase arising for the model (11) is the Néel type, with spins rotating in the radial direction, as opposed to a helix. The orientation of the stripe phase depends on the direction of the induced anisotropy of the DMI. This provides a control over the orientation of the stripe phase by changing the polarization of the applied radiation — a property which may be employed in future light-controlled magnetic logic gates.

Interestingly, the range of metastability of individual Skyrmions does not generally coincide with the phase boundaries. However, we find that metastable Skyrmions generally do not exist in a stripe crystal phase that is dominant at low values of the magnetic field. In this region, individual Skyrmions quickly become unstable with respect to stretching in the $y$ direction to form a stripe. For intense light with the parameter $\gamma$ exceeding the first zero of $J_{0}(\gamma)$, i. e., for $\gamma\gtrsim 2.4$, the phase diagram is dominated by the stripe phase at small magnetic fields. The Skyrmion crystal phase is limited to a moderate light intensity as shown in Fig. 2(b).

The numerical studies of Skyrmion dynamics in the absence of the field $\gamma=0$ were performed in Refs. Iwasaki-NN-2013 ; Iwasaki-NComm-2013 ; Iwasaki-Nano-2014 ; Koshibae-2016 ; ZZ-2016 . In the absence of light, i. e., for $\gamma=0$, the phase diagram of Fig. 2(b) reproduces these known results. Indeed, the obtained values of the critical fields for the transition between the stripe and Skyrmion-crystal phases ($H_{c1}=0.0072J$) and between the Skyrmion-crystal and ferromagnetic phases ($H_{c2}=0.026J$) are very close to those given in Ref. Iwasaki-NN-2013 . To simplify the comparison with Ref. Iwasaki-NN-2013 , we have used the same parameter values of the DMI strength $D_{x}=D_{y}=0.18J$ at $\gamma=0$.

In conclusion, the field of magnetic Skyrmions has attracted considerable attention due to the potential applications of Skyrmions in information processing. The major advantage of such noncollinear spin configurations as compared to domain walls is the possibility to make the Skyrmion size as small as a few nanometers without losing its stability. In this Letter, we employ the $s$-$d$-like exchange model for a weak two-dimensional ferromagnet with strong spin-orbit coupling to show that the off-resonant linearly polarized light can be used to tune the strength of the DMI and induce a large DMI anisotropy in the two orthogonal directions. This effect leads to the appearance of novel anisotropic phases — an elliptic Skyrmion crystal phase and a stripe phase — and can provide a new tool to control the stability, size, and shape of individual Skyrmions, as well as a control over the stripe phases by changing the light polarization direction.

The predicted effects may be observed in thin films such as Co/Pt using femtosecond laser pulses Huisman2016 . The light pulses must be sufficiently long to drive the structure into a nonequilibrium state that can be considered quasistationary. The typical experimental facility with the pump fluence 2-3 mJ/cm2 should be sufficient to test the theoretical results. We also expect that qualitatively the same physics persists at room temperature, helping to create a controlled set of Skyrmions that can be used to make the concept of Skyrmion racetrack memory viable Monchesky2015 .

We thank Alexey Kimel for helpful discussions. The support from the Russian Science Foundation under Project No. 17-12-01359, from the Dutch Science Foundation NWO/FOM 13PR3118, and the EU Network FP7-PEOPLE-2013-IRSES Grant No. 612624 “InterNoM” is gratefully acknowledged.

I Derivation of the effective Hamiltonian

The model of Eq. (1) of the main text is transformed by the unitary transformation (2a) to the model (3) which can also be written as

[TABLE]

where $J_{n}(\gamma)$ is the $n$-th order Bessel function of the first kind and $\gamma=2e\alpha A_{0}/\hbar\omega$.

If the frequency $\omega$ of a driving field is the largest energy scale in the system (the frequency $\omega$ is much larger than the bandwidth so that optical excitations of electrons are forbidden) one may replace the time-dependent problem with the time-independent Hamiltonian of an effective stationary model. The corresponding formalism of high-frequency expansion in the form of Brillouin-Wigner perturbation theory Mikami2016 has recently been developed for this class of problems. This formalism applied to the model of Eq. (s1) gives the following effective model

[TABLE]

where we kept terms up to the second order in $1/\omega$. Substituting $h_{n}$ from Eq. (s1) into the general expression of Eq. (s2) we obtain the effective Hamiltonian

[TABLE]

where $\beta=(\Delta^{2}m_{z}^{2}+[\alpha\hbar k_{y}-\Delta m_{y}]^{2})/(\hbar\omega)^{2}\ll 1$ under off-resonant radiation conditions. It is also clear that we can neglect the terms proportional to the sum of the Bessel functions in Eq. (s3). Moreover, if the parameter $\gamma$ does not lay in the immediate vicinity of zeros of the Bessel function $J_{0}(\gamma)$ (the second line in Eq. (s3)) the system in question can be effectively described (in the original basis) by the Hamiltonian

[TABLE]

which is the effective model of Eq. (5a) considered in the main text. Similarly, one can show that all high order terms arising in $1/\omega$ expansion are negligible away from zeros of $J_{0}(\gamma)$.

In a close vicinity $\delta\gamma$ of the zeros of $J_{0}(\gamma)$ the high order terms are formally important. However, it can be shown that the actual value of $\delta\gamma$ is negligibly small. To provide a quantitative estimate of $\delta\gamma$ we evaluate the quantity

[TABLE]

at the first zero $\gamma_{0}=2.4048$ of the function $J_{0}(\gamma)$. Thus, the effective model of Eq. (5) formally breaks down only in a tiny region around $\gamma_{0}$ that is determined by $\delta\gamma=|\gamma-\gamma_{0}|\sim 10^{-5}\beta$, which is well beyond our numerical resolution.

II Derivation of the DMI strength

In this section we present the derivation of antisymmetric exchange interaction which is linear in magnetization gradients and can be identified as the DMI. To analyze the polarization operator (8) we compute the quantity

[TABLE]

where $\varepsilon_{\bm{q}}^{\pm}=\hbar^{2}q^{2}/(2m)\pm\alpha\hbar qg(\theta_{\bm{q}})$, $g(\theta)=[\sin^{2}\theta+J_{0}^{2}(\gamma)\,\cos^{2}\theta]^{1/2}$, $f(\varepsilon)$ is the Fermi-Dirac distribution, $\hat{\bm{e}}_{x}=\hat{\bm{x}}$, $\hat{\bm{e}}_{y}=\hat{\bm{y}}$, and $\bm{n}_{\bm{q}}=\bm{q}/|\bm{q}|$. Furthermore, one can show that

[TABLE]

where we kept the notations of the main text. Taking the integral in Eq. (s6) we conclude that the quantities $N_{i}^{(1)}(\bm{k})$ in the linear order with respect to the momentum $\bm{k}$ are given by

[TABLE]

and

[TABLE]

Thus, the second order correction $\delta\mathcal{S}[\bm{m}]$ to the bare action $\mathcal{S}_{0}[\bm{m}]$ reads

[TABLE]

which coincides with Eq. (9) of the main text with renormalized DMI strength $D_{x}$, $D_{y}$ and Lifshitz invariants $\Lambda_{ij}^{(l)}$ defined by Eqs. (10) of the main text.

III Details of the numerical simulations

In our numerical calculations we use the lattice model given by the Eq. (11). To find the stationary states minimizing the total energy (11) we evolve the overdamped Landau-Lifshitz-Gilbert (LLG) equation

[TABLE]

until the stationary state is reached. Here, the effective field $\mathbf{H}_{eff}$ is defined by the functional derivative of the energy (11) over the local magnetic moment,

[TABLE]

The stationary state corresponds to the minimum of the total energy (11). To ensure that a particular state is a global minimum rather than the local minimum we compare the energies of essentially different solutions obtained for each known nontrivial phase (skyrmion crystal, stripe phase or a ferromagnetic phase). To solve numerically the LLG equation (s11) we use an explicit finite element method implemented in C and parallelized with OpenMP. Below we discuss how each of the two figures presented in Fig. 2 of the main text was obtained.

The Fig. 2a representing stability of a skyrmion is obtained as follows. We start from a seed approximate solution in the form of a skyrmion and evolve the LLG equation (s11) for sufficiently long time to obtain the numerically exact stationary solution for an isolated skyrmion. On the second stage, the stationary skyrmion solution is distorted by adding a random noise and by varying the parameters $\gamma$ and $H_{ext}$. The solution has been accepted as a stable if, upon addition of the noise and variation of the parameters it settles down and does not decay during sufficiently long evolution time of about $\sim 3000\,\hbar/J$.

To calculate the phase diagram in Fig. 2b we proceed as follows. Configurations minimizing (11) in a periodic rectangular domain of $N_{x}\times N_{y}$ cites are found as stationary solutions of the Eq. (s11) obtained by evolving (s11) from an initial seed solution for all known nontrivial phases (skyrmion crystal or stripe phase). In case of the modulated phases (skyrmion and stripe phases) the periodicity of the system for an infinite system is not known beforehand and should be found by minimizing the lattice parameters. To perform this task numerically we analyze a rectangular region of the skyrmion lattice made of $N_{x}\times N_{y}$ spins, which contains two skyrmions. The energy configuration minimizing the energy density $E/N_{x}N_{y}$ is then found by the coordinate descent method in the configuration space ($N_{x}$,$N_{y}$). The minimal energy configurations obtained for different values of the parameter $\gamma$ and for a fixed magnetic field $H_{ext}=0.015J$ are shown in Fig. s1 where a part of the skyrmion lattice with $210\times 210$ spins is shown. As seen from the Figure, the increase of $\gamma$ deforms the skyrmion configuration from the equilateral triangular lattice (at $\gamma=0.0$) to a isosceles triangular lattice (at $\gamma=0.5$ and $\gamma=1.0$). Similar procedure is implemented for the stripe crystal. In this case the problem is simplified since the optimization is required only along one direction (due to the tendency of stripes to align along the “easy” direction defined by the polarization of the applied radiation). The energies of the stationary solutions obtained are, then, compared to the minimal energy identified as a ground state.