Archimedean Screw in Driven Chiral Magnets

In chiral magnets a magnetic helix forms where the magnetization winds around a propagation vector $\mathbf{q}$. We show theoretically that a magnetic field $\mathbf{B}_{\perp}(t) \perp \mathbf{q}$, which is spatially homogeneous but oscillating in time, induces a net rotation of the texture around $\mathbf{q}$. This rotation is reminiscent of the motion of an Archimedean screw and is equivalent to a translation with velocity $v_{\text{screw}}$ parallel to $\mathbf{q}$. Due to the coupling to a Goldstone mode, this non-linear effect arises for arbitrarily weak $\mathbf{B}_{\perp}(t) $ with $v_{\text{screw}} \propto |\mathbf{B}_{\perp}|^2$ as long as pinning by disorder is absent. The effect is resonantly enhanced when internal modes of the helix are excited and the sign of $v_{\text{screw}}$ can be controlled either by changing the frequency or the polarization of $\mathbf{B}_{\perp}(t)$. The Archimedean screw can be used to transport spin and charge and thus the screwing motion is predicted to induce a voltage parallel to $\mathbf{q}$. Using a combination of numerics and Floquet spin wave theory, we show that the helix becomes unstable upon increasing $\mathbf{B}_{\perp}$ forming a `time quasicrystal' which oscillates in space and time for moderately strong drive.


I. INTRODUCTION
The Archimedean screw has benefited humanity as a mechanical tool since antiquity. There is evidence that it was already used in ancient Egypt to pump water, but even to this day, it is still used extensively, for example to transport materials such as powders and grains in factories. In addition, some bacteria use helical screws, so-called flagella, to propel themselves through liquids. Usually an Archimedean screw consists of a helical surface encased in a tilted tube, a simplified version of which is shown in Fig. 1(a). By rotating the screw on its axis as shown, the helical surface can be made to push material inside upwards, as indicated by the blue spheres and vertical arrows.
Helical surfaces analogous to the Archimedean screw have been predicted and observed in chiral magnets [1][2][3]. There the helical surface is spanned be spins winding around the corresponding pitch vector q, see Fig. 1(b). These structures form naturally in chiral magnets at low temperatures [4,5]. Chiral magnets are dominantly ferromagnetic materials in which inversion symmetry is broken by the crystal lattice, allowing weak spin-orbit interactions to induce a so-called Dzyaloshinskii-Moriya interaction. It is the competition between these two interactions that favors the formation of long-wavelength helical structures [3]. In addition to the helical phase, chiral magnets also host other phases. The conical phase can simply be viewed as a helical phase oriented parallel to an external magnetic field where spins uniformly tilt towards the magnetic field. In a small phase pocket close to the critical temperature T c , a skyrmion phase -a lattice of topologically quantized magnetic whirlscan form. Skyrmion phases can be manipulated by ultrasmall external forces created, e.g., by electric currents [6,7]. The coupling to currents is directly proportional to the winding number of skyrmions. This mechanism is absent for the topologically trivial helical and conical phases, which are therefore more difficult to control. It has also been suggested to use oscillating fields to move a single skyrmion [8][9][10], to create skyrmions [11] or to melt skyrmion crystals [12]. Similarly, the motion of domain walls by oscillating fields has been studied in simulations [13]. When a weak oscillating magnetic field is applied to a magnet, to linear order in perturbation theory, spin waves are excited. Early experiments by Onose et al. [14] showed that the helical, conical, skyrmion lattice and ferromagnetic phases exhibit a characteristic pattern of collective spin wave excitations. These excitations have been quantitatively described by linear spin wave theory in a range of different materials [15][16][17][18]. In the case of the helical and conical states at k = 0, the oscillating external field couples to two modes, often referred to as ±Q modes, for a review see [17]. They can be viewed as (spin-compression) waves traveling up or down the helix.
To second order in perturbation theory, a magnetic field oscillating with the frequency Ω is expected to generate a response at frequencies 0 and 2Ω. We will argue that the zero-frequency response couples to the Goldstone modes of the helical and conical phases. Here, naïve perturbation theory breaks down and a slow precessional motion with frequency Ω screw of all spins is induced as sketched in Fig. 1(b). This type of motion is precisely of the type characteristic of an Archimedean screw. Equivalently, the net rotation can also be interpreted as a translation with velocity V screw = λ Ω screw /(2π), where λ is the pitch of the helix. In the absence of pinning by disorder, this screw-like motion is induced for arbitrarily weak oscillating fields.
Upon increasing the strength of the driving field, the Archimedean screw solution ultimately becomes unstable. The discrete time-translational invariance of the driven system is spontaneously broken and an incommensurate spin wave oscillating in space and time is macroscopically occupied. Such a state can be viewed as a time crystal, or, more precisely, as a time quasicrystal as it is an incommensurate state [19,20]. In magnets such states are also referred to as magnon Bose-Einstein condensates (BECs). Such magnon BECs have, for example, been observed in YIG samples driven by GHz frequencies [21,22].
In the following, we will first analyze the equations of motion of spins in a helical or conical state driven by a perpendicular magnetic field B 1 (t) to second order in the amplitude O B 2 1 . We will show that a screw-like motion is induced and compare our analytical results to numerical micromagnetic simulations. To investigate the stability of the perturbative solution, we calculate the Floquet spin wave spectrum of the system and identify leading instabilities. Micromagnetic simulations show that these instabilities lead to the formation of a time quasicrystal at intermediate driving strength while chaotic behavior sets in at stronger driving. Finally, we show how the helical or conical state can be used as an Archimedean screw to transport electrons.

II. MODEL
We consider a chiral magnet in the presence of Dzyaloshinskii-Moriya interactions described by the free energy where F demag [M] encodes the dipole-dipole interactions and we use Heisenberg spins of fixed length, |M| = M 0 withM = M/M 0 . We consider an external magnetic field We will consider both linearly polarized fields, B y ⊥ = 0, and circular polarization, B x ⊥ = ±B y ⊥ . Throughout the paper we consider small oscillating fields and use 1 for bookkeeping purposes in perturbation theory.
We have performed all analytical and numerical calculations in the presence of dipolar interactions. To avoid overly long formulas, the analytical formulas presented in the main text are given in the absence of dipolar interactions. The effects of dipolar interactions are discussed in App. C.
In the absence of oscillating fields, = 0, the free energy for B 0 < M0D 2 J is minimized by the conical state described by where the helical pitch vector and the conical angle are given by q = D Jẑ and cos(θ 0 ) = B0M0J D 2 , respectively. When B 0 = 0 the free energy is isotropic and there is no preferred direction for q, but we are still free to choose the z-axis as the direction of spontaneous symmetry breaking, q z. In this case θ 0 = π/2, corresponding to a helical state where magnetization and q are perpendicular to each other everywhere.
The magnetic texture, Eq. (3), is translationally invariant in the xy plane. It is also invariant under a combined spin-rotation and translation along theẑ direction.

III. ARCHIMEDEAN SCREW
For an oscillating field, we calculate the time evolution of M(r, t) using the Landau-Lifshitz-Gilbert (LLG) equationṀ Here δM is functional derivative of the free energy Eq. (1), α is a phenomenological damping term and γ is the gyromagnetic ratio. Note that we use a convention where γ is negative, with γ = − |e|g 2me for an electron with charge −|e|, mass m e and g-factor g. The prefactor in front of the damping term ensures that all formulas remain valid independent of the sign of γ.
The goal is to calculate the response to the oscillating magnetic field, Eq. (2), by doing a Taylor expansion in . To this end we now update the parametrisation of M given in Eq. (3) to allow for some small dynamical excitations, replacing This parametrization assumes that the system remains translationally invariant in the xy plane. Importantly, we assume here that the q vector does not tilt in the presence of the oscillating magnetic field. This is justified as such a tilt would nominally lead to frictional forces diverging in the thermodynamic limit. Experimentally, such a tilting occurs only for extremely slow changes of the field direction [23].
Substituting Eq. (5) into Eq. (4) and dotting with ∂M ∂θ , ∂M ∂φ gives two sets of coupled differential equations for θ 1,2 , φ 1,2 . To first order we get The equations to second order in take the form with c = cos(θ 0 ) and s = sin(θ 0 ). Note that we switched to dimensionless units b and also use dimensionless space and time units: z → q −1 z and t → JM0 D 2 |γ| t, respectively. The latter also motivates the definition of a dimensionless driving frequency ω = JM0 D 2 |γ| Ω. Eq. (6) and (7) are to be solved consecutively, as the first order solutions θ 1 , φ 1 enter in the second order equations.
To linear order, O( 1 ), the driving terms proportional to b x,y on the right hand side (RHS) of Eq. (6) have (dimensionless) Fourier momentum and frequency components ±1, ±ω. The steady state solutions of θ 1 , φ 1 are composed of these Fourier components only and therefore have the form which translate physically to two traveling waves running up or down the helix, depending on the relative sign in ωt ± z. The analytical forms of the pre-factors θ For circular polarized driving, b x = ±b y , only one of the traveling wave modes gets excited. When dipolar interactions are switched on this only remains true if the demagnetization factors N x , N y (see App. C for definition) in the plane perpendicular to q are identical. Right and left polarized circular driving are defined as the magnetic field rotating anticlockwise and clockwise in time, respectively, when we position ourselves at the origin and look in the positiveẑ direction. If we drive with a right polarized magnetic field b x = b y only the down-traveling (ωt + z) wave will be excited, and vice versa for left polarized driving.
To second order, O( 2 ), the coupled equations in Eq. (7) have driving terms with Fourier components k = 0, ±2, ω = 0, ±2ω. Here the mode k = 0, ω = 0 is special as it couples to the Goldstone mode of the system, arising from the spontaneously broken translational symmetry of the conical state. Therefore, we can expect a diverging response. If we substitute the naïve choice of the (k = 0, ω = 0)-Fourier modes θ 0,0 2 , φ 0,0 2 , independent of t, z into the first equation of (7), we quickly run into trouble, as the left and right sides of the equation do not balance each other. Mathematically, this conundrum can be solved by assuming that φ 2 obtains a correction linear in t φ 2 (z, t) = φ osc 2 (z, t) + ω screw t.
Physically, this term does not describe an instability of the system but a rotation of the helix with angular velocity ω screw which induces a screw-like motion. As shown in Fig. 1(b), individual spins precess rapidly with the driving frequency Ω (small circles in Fig. 1(b)). In analogy to the physics of a spinning top, this rapid local motion triggers a slow net precession of all spins around the q axis giving rise to a rotation of the helix with frequency ω screw . Equivalently, this screw-like rotation can also be interpreted as a translation of the helix in space parallel to q with constant velocity, Within our perturbation theory ω screw is quadratic in the oscillating fields. An analytic formula for ω screw is given in App. B. In the absence of a constant magnetic field, B 0 = 0 and θ 0 = π/2, we obtain where b R/L = b x ±b y are the amplitudes of the right-and left polarized oscillating magnetic field. In the second line of Eq. (11) we expanded around the resonance frequency ω res = √ 2 + O(α 2 ) in the limit of small damping α. Translating this back to physical units, Eq. (11) reads Ω screw ≈ Ω res 3 32 where Ω res = √ 2 D 2 |γ| JM0 . In Fig. 2(a) we show ω screw as function of the (dimensionless) driving frequency ω for a vanishing external field. Switching from left-to right-polarized oscillating B-fields changes the sign of ω screw . At the resonance frequency of the the helix ω screw is strongly enhanced in the limit of weak damping by the factor 1/α 2 . For linear polarization b y = 0 or, equivalently, b R = b L , there is no rotation of the helix, ω screw = 0, as predicted by Eq. (11). This changes when a static magnetic field parallel to q is switched on, see Fig. 2(b). In the resulting conical state the response to right-and left-polarized fields become different, see App. B, and one also obtains a finite result for linearly polarized fields oscillating only in the x direction. In this case one can control the sign of ω screw by changing the direction of the field B 0 .
All formulas above are given in the absence of dipolar interactions. If they are included, an analytical calculation is still possible but the resulting formulas are too long to be displayed. The analytical result is plotted in Fig. 2(c) for vanishing external field (helical state) and in occurs when both a static external field B 0 and dipolar interactions are considered together. In this case the resonance splits into a right-handed and a left-handed mode which selectively couple to the right-and left-polarized oscillating fields. If in this situation a linearly polarized oscillating field is considered, b y = 0, one can control the sign of ω screw by changing the frequency of the applied field, see Fig. 2

(d).
To confirm our results, we employ micromagnetic simulations. Using mumax3 [24,25], we solve the LLG Eq. (4) numerically for a conical state driven by an oscillating magnetic field. Parameters are chosen to describe Cu 2 OSeO 3 , where we choose the damping parameter to be α = 0.03. For a quantitative comparison between simulations and analytical calculations (both including the effects of demagnetization fields), we determine Ω screw as a function of driving frequency Ω. From a set of simulations with different excitation frequencies Ω, we extract Ω screw as the linear slope of the azimuthal angle φ(t) of a single spin. For the chosen parameters rotation frequencies Ω screw are in the MHz range, for driving frequencies in the GHz range. In Fig. 3 we compare the numerical result to the analytical formula and find an excellent agreement.

IV. FLOQUET SPIN WAVE THEORY
As we will discuss below, the Archimedean screw solution becomes unstable when the driving fields get too large. This motivates us to investigate the stability of our solution using spin wave theory, or, more precisely, the "Floquet" variant of spin wave theory, which can be used to describe periodically driven systems. For this we have to expand the magnetization around the (perturbative) solution (5), M = M screw + δM to derive an equation for δM. In the following we use a notation similar (but not identical) to the one which is familiar from the Holstein-Primakoff treatment of quantum spins in the large S limit [26]. Importantly, we will also include the effects of the phenomenological damping α, which cannot easily be described by a quantum Hamiltonian. The magnetization is parametrized by where a(r, t) and a * (r, t) are complex space-and timedependent expansion coefficients. We use a coordinate system where e 3 points parallel to the local magnetization of the Archimedean screw solution while e ∓ are perpendicular, with where the angles θ(r, t) and φ(r, t) are given by the solutions (5) discussed in Sec. III. The expansion coefficients a and a * have Poisson brackets to leading order in a Taylor expansion in a. Using the notation of classical Hamiltonian dynamics, the LLG equation (4) takes the form We will only be interested in terms linear in a and a * and up to quadratic order in the oscillating external fields. More precisely, we consider to quadratic order only the contributions giving rise to the screw-like motion, omitting tiny oscillating terms at frequencies 2ω. A useful check of the expansion (and the Archimedean screw solution of Sec. III) is that all constant terms O(a 0 ) drop out. Projecting the resulting equation onto the directions e ∓ , see App. D, gives rise tȯ where F (2) is the contribution to F quadratic in a and a * , and the factor of |γ| has been absorbed into F (2) . The fact that we have periodic driving and are expanding around a state that is -in a frame of reference co-moving with our Archimedean screw -periodic in space and time makes Eq. (16) an ideal candidate for a Floquet treatment. We begin by defining the space and time Fourier transformed fieldsã m k ,ã m * k as Eigenvalues λ k of M F as a function of k , the component of k q, for a system with parameters α = 0.03, c = 0.71, δ = 1.76, Nx = Ny = Nz = 1/3, k ⊥ = 0. All graphs are in the first Brillouin zone −q/2 < k < q/2. The real parts of the eigenfrequencies Re[λ] are plotted in the first Floquet zone, between −ω/2 < Re[λ] < ω/2 with ω = 2. The two graphs in the left column are the real and imaginary parts of λ for an undriven system. All imaginary parts are negative, indicating that the conical static state is stable, as expected. The two graphs in the right column show the band energies for a driven system where the driving magnetic field is left circular polarized: bL = 0.01, bR = 0, with driving frequency ω = 2 very close to the resonance frequency ωres,+. The parts of the imaginary spectrum highlighted in red are unstable, and occur for k ∼ 0.13q, which corresponds to a Re[λ] = ±0. 16. The crosses denote the spectrum at k = 0, which differs from the spectrum for k → 0 due to the long-ranged dipolar interactions.
Note the factor r + v screw t arising from our comoving coordinate system. Within our perturbative scheme, only the fieldsã m k+nq ,ã m * k+nq with indices m = −1, 0, 1 and n = −1, 0, 1 couple to each other. We collect those in a 18-component vector Ψ F k . The restriction of the Floquet space is formally justified because we investigate the system in the limit of small B ⊥ and our results for eigenenergies and decay rates are formally exact to quadratic order in B ⊥ . Here it is important to realize that to conserve the Poisson brackets of the fields, one has to perform Bogoliubov transformations to diagonalize the dynamical matrix describing our system. Taking this into account, it is possible to recast Eq. (16) as a 18 × 18 matrix equation (see App. D for details) with Ψ F k (t) = e −iλt Ψ F k . Importantly, the Floquet-Bogoliubov matrix M F is not a Hermitian matrix, both because of the underlying Bogoliubov transformation and the damping terms. Its eigenvalues are therefore complex in general, There is a clear physical interpretation for the real and imaginary parts of λ: the real part gives the temporal frequency of oscillation of the spin wave, whereas the imaginary part determines how fast it grows or decays in time. Importantly, the sign of the imaginary part determines whether the spin wave decays (negative imaginary part) or grows (positive imaginary part) exponentially in time, signaling an instability. As we show below, such instabilities are quite common in driven bosonic systems.
To understand how the oscillating fields affect the spin wave spectrum it is useful to consider first the case without oscillating fields, B ⊥ = 0, shown in the left panels of Fig. 4. The upper left panel shows the real parts of the eigenmodes, Re[λ k ], as function of the momentum k parallel to the q direction. They always come in pairs i.e., they are calculated modulo the driving frequency ω. The lower left panel displays the imaginary parts, Im[λ k ], which in the absence of driving are always negative and describe the decay of modes due to the damping term α. Note that for k → 0 the Goldstone mode becomes overdamped and purely diffusive: the real part vanishes and λ k ∼ −iαk 2 . This is the behavior expected for Goldstone modes in systems where translational symmetry is spontaneously broken but where at the same time the underlying model lacks momentum conservation [27][28][29].
The panels on the right of Fig. 4 show how a finite oscillating field modifies the spin wave spectrum. Here the most dramatic effect occurs for the imaginary parts in the lower right panel: when the oscillating field is sufficiently large, they change sign and become positive. Thus the system becomes unstable when the oscillating field increases. The physics of the instability can be traced back to a resonance described by a simple 2 × 2 matrix Here 0 i,k > 0 denotes the energies of spin waves with band index i of the unperturbed system and αΓ i are the corresponding lifetimes. The frequency-dependent prefactors µ (i) ω describe how the oscillating fields couple the energy level on the diagonal of the matrix. The coupling is most efficient when the driving frequency hits a k = 0 resonance of the helix. Schematically, we find The instability is most pronounced when In this case the oscillating field can resonantly create a pair of spin waves out of the vacuum. In contrast, we do not find instabilities at energies 0 i,k − 0 j,k = ±ω when spin waves are resonantly coupled. At this spin wavecreation resonance, the eigenvalues of M res are given by Importantly, the sign of Im[λ + res ] changes when b ⊥ grows, signaling an instability. Assuming Γ 1 ∼ Γ 2 ∼ Γ the system is only stable if (up to numerical prefactors) More precisely, this formula is only valid for ω ≈ ω res . If one stays away from this point, then µ In the limit α → 0 our calculation predicts that an arbitrarily weak oscillating field induces an instability. This is, however, an artifact of our approximation which ignores that the modes with finite energy and momentum can also decay via scattering processes. In this case an extra calculation of these lifetimes would be necessary to estimate when the instability occurs. As a function of momentum, the resonance condition Eq. (21) is met along planes in momentum space. We therefore have to find the leading instability, i.e., the one where upon increasing b ⊥ the instability occurs first. In Fig. 5, Im[λ k ] is shown as a function of k ⊥ and k . At least for the parameter regime investigated by us, we find that the dominant instability occurs for k ⊥ = 0.
To track the leading instability as function of frequency, we plot in Fig. 6 the largest Im[λ] as a function of driving frequency ω, for a range of amplitudes of linearly polarized driving b x . To produce this figure, we diagonalized M F at k ⊥ = 0, choosing k to fulfil the resonance condition Eq. (21). As expected from Eq. (23) and (24), the system first becomes unstable at the resonance frequencies ω ± res of the underlying conical state.

V. FORMATION OF A TIME QUASICRYSTAL
The spin wave calculation rigorously shows that for α > 0 the Archimedean screw solution is stable for small amplitudes of the oscillating field but becomes unstable upon increasing the field strength. However it cannot predict the fate of the unstable system. Therefore we again used numerical solutions of the LLG equation to analyze this regime. As the instability is expected to occur at finite momentum k , it is essential to make the system sufficiently large in this direction. We therefore simulated a system with a length of up to 15 times the pitch of the helical state. In the perpendicular direction we use periodic boundary conditions using the fact that the instability occurs at k ⊥ = 0, see Fig. 5. In very good agreement with our analytical solution, we find the stable Archimedean screw solution for small amplitudes of the driving field as discussed above in Fig. 3. We can interpret this new mode as a kind of laser-type instability (or, equivalently, as a Bose-Einstein condensate) of the resonantly driven magnons. As the mode oscillates in time and space it defines a "time crystal", or more precisely, a "time quasicrystal", as the frequency and momentum of oscillation determined by Eq. incommensurate with the driving frequency ω and the pitch vector q of the underlying conical state. From the viewpoint of symmetry, due to the presence of the oscillating field, time-translation invariance is only discrete. This discrete symmetry is then spontaneously broken by the time quasicrystal.
In Fig. 8 we show the screwing frequency as a function of the amplitude of the oscillating magnetic field. For small amplitudes, Ω screw grows quadratically in B x ⊥ , following exactly the prediction of perturbation theory. The screwing frequency continues to grow in the regime where the time quasicrystal forms but the rate of growth is strongly reduced.
When we increase the driving further, the time quasicrystal also becomes unstable, see Fig. 8 and Fig. 7. We enter a chaotic regime discussed in more detail in App. E. Note that our simulations are not reliable in this regime as they assume translational invariance in the direction perpendicular to the q vector, which is valid both for the Archimedean screw solution and the time quasicrystal, but not in the chaotic regime.

VI. TRANSPORT
Archimedean screws have been widely used for technological applications since antiquity, for example to transport water in irrigation systems, dehumidify low lying mines, or more recently even to deliver fish safely from  Fig. 3, with driving frequency f = 4.15 GHz) for two different sizes of the simulated system (7 and 15 times the pitch of the helix). For small fields Ωscrew, the Archimedean screw solution is found, following the analytic prediction which high accuracy. Similarly, an instability resulting in the formation of a time quasicrystal occurs as predicted. The onset of the instability is slightly delayed for the smaller system, as the predicted wavelength of the instability λ ≈ 7.7 2π q does not match the boundary conditions in this case. For B x ⊥ 3.8 mT we obtain chaotic solutions discussed in more detail in Appendix E).
one tank to another in so-called "pescalators" on fish farms. But could they also be used for transport in our system? In this section we want to show how coupling electrons to our rotating helical magnet gives rise to a finite DC current parallel to the q vector of the magnet. We model the electronic system by the following Hamiltonian where C(r) = (c ↑ (r), c ↓ (r)) T is a spinor containing the up and down components of the electron annihilation operators. In addition to the free energy termp 2 /2m we have a spin-orbit coupling term λ sop ·σ and the exchange coupling n · σ of the electrons' spins to the local magnetization n. For a static helix, the spin-orbit term induces the formation of exponentially flat mini-bands of periodicity q in the k direction [30]. As we want to study the transport of electrons, it is essential to include the effects of disorder, which we model by a spin-independent random potential V (r). In the following, we will model the effect of scattering from disorder by a scattering rate 1 τ . We assume the following hierarchy of energy scales, F > J H τ , λ so k F , typical for magnets with weak spin-orbit coupling, where F and k F are the Fermi energy and Fermi momentum, respectively.
We are interested now in a moving helix. We use a simplified Archimedean screw ansatz for n(r, t) where we have suppressed all the oscillations which are multiples of the driving frequency ω and kept only the ω screw time dependence.
In the absence of disorder (and also in the absence of Umklapp scattering due to electron-electron interactions), the problem can be solved by moving to a frame of reference comoving with the helix using the transformation C † (r) → C † (r − v screw t). The current in the comoving frame vanishes and therefore the electronic current density j in the lab frame is simply given by where n ↑/↓ are the electron densities of majority and minority electrons, respectively. More realistically, one has to take into account the effects of disorder (or Umklapp scattering) which is expected to dominate transport properties. Here it is useful to consider a transformation where (i) impurities do not move, and (ii) the Hamiltonian is diagonal in the dominant term J H , i.e. the spins of the electrons are aligned and anti-aligned with the time-dependent local magnetization n(r, t). This can be achieved by rotating the spin-quantization axis using the unitary matrix U U (r, t) = cos(θ 0 /2) sin(θ 0 /2)e −iφ sin(θ 0 /2)e iφ − cos(θ 0 /2) , where φ = qz − ω screw t. We can then define C(r) = U (r)D(r), such that d † ↑ , d † ↓ now create electrons with spins parallel and anti-parallel to the local time-dependent magnetization n, respectively. Rewriting Eq. (25) in terms of D, D † and switching to Fourier space we obtain approximatelỹ Here we ignored some small static correction terms to σ,k as well as spin-mixing terms of type d † ↑ d ↓ , which can be ignored because of the large splitting between the minority and majority-spin Fermi surfaces due to J H . Importantly, the unitary transformation does not affect the potential scattering term H dis .
Following the rotation by U , the only time-dependent term H 1 (t) in the Hamiltonian comes from spin orbit interactions. For H 1 = 0, we obtain electrons with σk describing majority and minority electrons. Their Fermi surfaces are shifted by k = ±k 0 both due to the spinorbit interactions and the rotation of the spins by the matrix U .
We would like to evaluate the expectation value of the parallel component of the current operator treating H 1 (t) as a small time dependent perturbation. We can formulate this as a Keldysh problem 1(t )dt and we denote asÕ(t) = e iH0t Oe −iH0t the operators in the interaction picture. Ultimately we are interested in the DC component of J , which to lowest order comes in at second order in the perturbation H 1 (t) ∼ e iωscrewt . After some algebra (see App. F for technical details) we arrive at where we have used that ω screw = qv screw . Here, n σ,k is the Fermi distribution function (1 + e β( σ,k − σ,k F ) ) −1 .
Performing the integral in k-space at T = 0 amounts to integrating over the two Fermi spheres located at ±k 0 discussed earlier (see App. F for details). We obtain In the limit when the mean free path of the electrons v F τ is smaller than the wavelength of the helix 2π q , the current is quadratically dependent on the electron's lifetime τ . In the opposite limit, v F τ 2π/q, in contrast, the current is linear in τ and thus proportional to the conductivity of the system. Eq. (34) has been derived in perturbation theory in λ so and thus cannot describe the formation of band-gaps and minibands triggered by λ so . These minibands have a band splitting of the order of ∆ ∼ q √ v F λ so [30] and thus perturbation theory is only reliable for τ ∆/ 1 which sets an upper limit for the regime of validity of the second line in Eq. (34). These results are summarized in Fig. 9.

VII. EXPERIMENTAL SIGNATURES AND CONCLUSIONS
Within our numerical and analytical calculations, we found that even for a weak oscillating magnetic field, the electron lifetime τ electronic current j /envscrew FIG. 9. Schematic plot of the electronic current density j as a function of electron lifetime τ . Note that we have suppressed the spin indices, which is justified for a strongly spin polarized system N ↑ N ↓ . For a strongly disordered system, when is quadratic in τ . In the range (qvF ) −1 τ (q √ vF ) −1 , j grows linearly with τ . Our perturbative assumptions break down in the dashed region, but we know that for a very clean system with no disorder (τ 1) the current must plateau at j = envscrew.
magnetic helix starts to rotate in a screw-like motion. This means that naïve perturbation theory breaks down as the difference, M(r, t) − M 0 (r), of the magnetization of the perturbed system, M(r, t), and of the unperturbed state, M 0 (r), grows linearly in time. Physically this arises because the system couples to a Goldstone mode and technically it can be described by using the moving helix as a starting point of perturbation theory. Similar effects also arise in many other systems. For example, one can move skyrmions by oscillating fields [8,10] and ratchets also work by a similar mechanism, see [31] for a theoretical review and [32,33] for experiments. Friction plays a decisive role for this phenomenon. Both the force which induces the rotation of the helix and the counter force arising from the motion of the helix are proportional to the friction coefficient. As a result, the frequency Ω screw describing the screw-like rotation obtains a finite value in the limit of vanishing friction constant, α → 0. A second important effect is that friction is needed to stabilize the state and to avoid instabilities and the onset of chaos in this driven nonlinear system. The net effect is that one can reach larger values of Ω screw in systems with stronger friction. Here both extrinsic friction (parametrized by α) arising from coupling to phonons or electrons and intrinsic friction arising from magnon-magnon scattering play a role but only the first effect was included in our Floquet spin wave theory of the instabilities.
In experimental systems the role of pinning by disorder has to be considered. The following order-of-magnitude estimates are motivated by the parameters in MnSi, arguably the best investigated chiral magnet. In the presence of pinning, we expect that a critical strength of the oscillating field is required before the helix starts to move.
To obtain a rough estimate, we assume a screw frequency Ω screw ∼10 MHz (obtained using Eq. (10), for micromagnetic parameters J = 7.05×10 −13 J m −1 , D = 2.46×10 −4 J m −2 , M 0 = 1.52 × 10 5 A m −1 corresponding to MnSi [15,34,35], as well as ω screw = 0.0007 for oscillating fields of the order of 0.5 mT and α ∼ 0.01). For a helix with a pitch of 200 A, this corresponds to a speed of v screw ≈ 200 mm s −1 . We can compare this speed to the velocity of skyrmions driven by a current j, in MnSi [6]. Skyrmions are expected to have a very similar friction and pinning compared to the helical and states as the magnetization is modulated on the same length scale. They start to move above a critical current density, j c , and their speed can be estimated from measurements of the Hall effect [6]. For example, at a current density of 2j c , the skyrmion velocity has been estimated to be about 0.2 mm/s, which is three orders of magnitude smaller than our estimate for v screw . We conclude that at least for resonant driving one can likely induce the screw-like motion of the helix in materials with low pinning as realized in MnSi and similar materials.
To detect the rotation of the helix, one can try to pick up a signal from the rotating magnetization using, e.g., a detector on the surface of the crystal. A more intriguing approach would be to observe the Archimedean screw "in action". For example, in a metallic system we have shown that it can transport charge. We therefore expect that a voltage will build up parallel to the orientation of the helix. The current and voltage will depend sensitively on the amount of disorder and the strength of spin-orbit coupling in the system. For example, in the chiral magnet CoGe spin-orbit interaction lead to a band-splitting of almost 10% of the bandwidth [36] and similar values are expected for MnSi. Thus we estimate λ so /v F ∼ 10 −2 − 10 −1 . Furthermore, MnSi can be grown with exceptional crystal quality and residual resistivities well below 1 µΩ cm, giving rise to a mean free path larger than 1000 A at low T [37]. Assuming a mean free path of the order of the pitch of the helix and using n ∼ 4·10 22 cm −3 [38] our calculation yields current densi-ties of order 10 4 -10 7 A m −2 . Even for the smallest values in this range, the corresponding voltage building up in such a system will be very easy to detect. In good metals, however, the skin depth (the length scale on which electromagnetic fields penetrate the sample) is only of the order of 1 µm at microwave frequencies. Therefore one should either use thin samples or bad metals. An interesting alternative is to try to detect thermal gradients or gradients in the magnetization arising from the transport of heat and spin, respectively.
For stronger driving, we predict the formation of a 'time quasicrystal' arising in the driven system. This can probably be detected most easily by picking up the radiation arising from the oscillating magnetization which is expected to be in the 100 MHz -1 GHz range. The detection of any monochromatic emission with a frequency smaller than the driving frequency is a unique signature of such a state.
In conclusion, we have shown that using the helical and conical states of chiral magnets and weak oscillating fields one can realize an Archimedean screw on the nanoscale. As one of the archetypal machines known to mankind, it can be used to explore the transport of charge, spin or heat in a novel setup.

ACKNOWLEDGMENTS
We thank Joachim Hemberger, Christian Pfleiderer, Andreas Bauer, Markus Garst and, especially, Yuriy Mokrousov for useful discussions. NdS also thanks S. Mathey and V. Lohani for helpful discussions. We acknowledge the financial support of the DFG via SPP 2137 (project number 403505545) and CRC 1238 (project number 277146847, subproject C04). We furthermore thank the Regional Computing Center of the University of Cologne (RRZK) for providing computing time on the DFG-funded (Funding number: INST 216/512/1FUGG) High Performance Computing (HPC) system CHEOPS as well as support.
with s = sin θ 0 and c = cos θ 0 . Note that in the special cases b x = ±b y , one of each pair of pre-factors θ vanishes. Physically, b x = ±b y correspond to right and left circular polarized driving, respectively. Circular polarized driving couples only to one of the two modes. This motivates us to define In these new variables, we get left circular driving by setting b R = 0, right circular driving by setting b L = 0, and linearly polarized driving in the x, y directions by choosing b L = ±b R . Any other choice of b L , b R corresponds to the general elliptical drive.
Without dipolar interactions, the two modes are degenerate with resonant frequency ω res = √ 2 − c 2 . To evaluate the screwing frequency we substitute Eq. (B1) into the first equation of Eq. (7). Here ∂ t φ 2 = ω screw balances all the other DC components in the equation, which can be computed from the first order solutions (∂ t θ 2 and φ 2 do not contribute as they are both oscillating in time and/or space). We obtain When dipolar interactions are switched on (see App. C for details and definitions), the single resonance frequency ω res gets shifted and split into two different resonance frequencies ω ± res proportionally to δ, a dimensionless measure of the strength of the dipolar interactions The prefactors θ describing the linear-response solution with dipolar interactions are lengthy and therefore not listed here. If the shape of the crystal is cylindrically symmetric around the axis of the helix, N x = N y , circular polarized light couples only to a single mode. One can analytically calculate the screwing frequency to second order in the oscillating fields. Instead of showing the exact result of this lengthy calculation, which is too long, we display below the most singular contribution obtained from a Taylor expansion around the two resonance frequencies ω ± res ω screw ≈ ∓ sgn(γ)b 2 R/L A sgn(γ)± ω − ω sgn(γ)± res 2 + ∆ω 2 , ∆ω 2 = α 2 ω sgn(γ)± res 2 c 2 (5δ + 6) − 7δ − 18) 2 4(δ + 6) (c 2 (5δ + 6) − 6(δ + 2)) . (B6) The prefactor A ± turns out to be finite in the limit α → 0 and therefore ω screw ∝ 1/α 2 at resonance. If the system is driven with |ω − ω ± res | ∆ω, in contrast, ω screw ∼ (ω − ω ± res ) −2 remains finite in the limit α → 0.

Appendix C: Calculation of dipolar interactions
In contrast to the Heisenberg and DMI energy terms which are local, dipolar interactions are long-ranged. They are not only much more computationally costly to calculate but require different treatment for the case k = 0 (the so-called demagnetization field limit) and the limit k → 0 in the thermodynamic limit [39]. Here the calculation for k = 0 has to take into account the energy stored in the magnetic fields outside of the sample.

Demagnetization fields
Let us denote the k = 0 or DC component of the magnetization as M i with where V is the total volume of the sample. For our helical or conical texture, this integral only needs to be done over the z-axis due to the translational invariance in the x, y directions. Applying Eq. (C1) to the static helical or conical ansatz Eq. (3), we find that M x = M y = 0 and the only non-zero DC component is M z = M 0 cos(θ 0 ). In general, the magnetization of the system will set up internal demagnetization fields in directions opposite to the applied external magnetic fields. For a sample with an ellipsoidal shape, the mathematical expression for these internal demagnetization fields is particularly simple and given by where N i are the demagnetization factors which solely dependent on the shape of the sample and obey the identity Tr(N ) = 1. The corresponding contribution to the free energy Eq. (1) is Applying this to the static conical helix, we obtain an additional contribution F dip,k=0 = 1 2 µ 0 N z M 2 0 cos 2 (θ 0 ). For the static conical state, q is unaffected but cos(θ 0 ) = b0 1+δNz changes. As δ, N z > 0, this means that θ 0 increases. This is a consequence of the internal demagnetization field B demag opposing and reducing the applied field B 0 .
For the dynamical calculation, we need to add Eq. (C2) to B eff in the RHS of Eq. (4), but now we need to substitute the dynamic ansatz Eq. (5) with Eq. (8) into M. This gives many new terms on the RHS of the first and second order equations Eq. (6) and (7). Here are the first order in contributions (C4) Using this result, one can calculate the corresponding magnetic fields using Eq. (C2) which contribute to the effective magnetic field in the LLG equation, Eq. (4). Mathematically, the oscillating finite-k contributions of θ 1 (z, t), φ 1 (z, t) modify the k = 0 magnetization because we are Taylor expanding around a spatially modulated static helix. At second order 2 we have many more terms because there are more combinations between the perturbing terms and static solution which modify the k = 0 magnetization. We do not list them here as they are lengthy, but the procedure to obtain them follows exactly from that used for the first order terms Eq. (C4).

Finite k contributions
At finite momentum, for k larger than the inverse system size, the contributions of the dipolar interactions to the free energy take the form where M k = 1 V d 3 rM(r)e −ik·r . We will first analyze how this term affects the static conical helix. The Fourier transform of Eq. (3) is therefore, since the only non-zero Fourier components of magnetization M ±q ⊥ q, F dip,k =0 = 0 for the static conical helix. Thus only the DC k = 0 components of magnetization play a role in the determination of q, θ 0 for the static conical helix.
Moving on to dynamics, we need to extract a magnetic field  (15). To calculate the excitation spectrum of the Archimedean screw solution, we first project Eq. (15) onto e ± , defined in Eq. (14). We use the Holstein-Primakoff expansion Eq. (13) keeping only terms linear in a, a * on both sides of the equation. Note that the terms which are zeroth order in a, a * simply correspond to the LLG equation, which we solved correctly up to second order in amplitude of driving b x , b y in Sec. III. The following identities for the three basis vectors e 3 , e ± turn out to be useful e 3 · e 3 = 1 e ± · e ± = 0 e ± · e ∓ = 1 e ± ·ė ∓ = ±i cos(θ)φ e ± ·ė ± = 0 e ± · (ė 3 × e ∓ ) = 0 e ± · (ė ± × e 3 ) = 0 e ± · (ė ∓ × e 3 ) = cos(θ)φ.
On the RHS of the projected Eq. (15) we have the Poisson bracket {F,M}. As we only consider F (2) -the contribution to F quadratic in a, a * -the only possibility for linear terms coming from the overall Poisson bracket is when we take the Poisson bracket with the linear in a, a * components S, i.e. e − a + e + a * . After projecting onto e + we obtain e + · i{F (2) , e − a + e + a * } = i{F (2) , a} Next we consider the damping term −αŜ×Ṡ. Here there are two possibilities for obtaining linear terms in a, a * : either we take a linear term from the firstŜ and zeroth term fromṠ, or vice versa giving e + ·(e − a ×ė 3 + e 3 × (ė − a + e −ȧ +ė + a * )) = iȧ−cos(θ)φa Setting the LHS equal to the RHS we obtain sgn(γ)(ȧ + i cos(θ)φa) = i{F (2) , a} − iα(ȧ + i cos(θ)φa).
Multiplying both sides of the equation by sgn(γ)−iα 1+α 2 we obtain the equation of motion for a Eq. (16) given in the main text. The equation of motion for a * can be obtained by following the same procedure as above, but projecting onto e − rather than e + , or simply by noticing that it should be the complex conjugate of Eq. (16).

Derivation of Floquet matrix M F
The goal of this subsection is to explain how we obtain the Floquet equation Eq. (18) from the equation of motion for the a, a * Eq. (16). The first step is to substitute the Fourier space and time expansions Eq. (17) into a, a * Eq. (16). The back Fourier transform lets us express a, a * as where k = k + k ⊥ and we are working in cylindrical coordinates r = (ρ, z) T , with ρ = (x, y) T . Due to the cylindrical symmetry of the problem, the azimuthal angle between k x , k y makes no difference and can be set to 0. Also note that k is only defined in the first Brillouin zone, −q/2 < k < q/2. It is also useful to define the column vector Ψ from which we the build the Floquet vector Ψ F , Eq. (D3).
where we definedz = z + v screw t and f (l) = 1 2 (l − lmax 2 ) , where l runs between l = 1 and l = l max , and l max stands for the maximal index of Ψ. It is sufficient to choose l max = 6 to obtain equations which are accurate to second order in the oscillating fields. l max is always even because we always include the same number ofã m k andã m * k operators. For theȧ expression we sum over only odd l = 1, 3, . . . l max − 1, whereas for theȧ * expression we sum over even l = 2, 4, . . . l max .
Let's now look at the RHS of Eq. (16). First we have to compute the Poisson bracket {F (2) , a/a * }. As previously mentioned, F (2) is obtained by inserting Eq. (13) into Eq. (1) and keeping only the terms quadratic in a, a * . By using the Fourier convention Eq. (17) we obtain F (2) in terms of theã m jq+k ,ã m * jq+k operators. F (2) contains both number conserving operatorsã m * jq+kã n lq+k and non-number conserving operatorsã m jq+kã n jq+k ,ã m * jq+kã n * jq+k . In general, this type of Hamiltonian can be diagonalized by Bogoliubov transformations, and the method we will use implicitly accomplishes the same thing. We denote the Fourier components of F (2) byF n (k). With this convention and the vectors Ψ m (k) defined in Eq. (D2) we obtain where j is odd if we are evaluating for the Poisson bracket with a and even if it is the Poisson bracket with a * . Above we used Ψ m i (k), Ψ n * j (k ) = (−1) i−1 δ i,j δ m,n δ k,k Ψ m i (k), Ψ n j (k ) = (δ i∈odd:i,j−1 − δ i∈even:i,j+1 )δ m,−n δ k,−k Ψ m * i (k), Ψ n * j (k ) = (δ i∈even:i,j+1 − δ i∈odd:i,j−1 )δ m,−n δ k,−k . The final remaining term on the RHS of Eq. (16) is theφ cos(θ) term. This term, which is built from the steady state solutions θ(z, t), φ(z, t), oscillates in both space and time with components e imωt , e inqz , m, n ∈ Z, and can be written asφ cos(θ) = m,n∈Z e −i(mωt+nqz) g m n . (D11) Finally, putting together Eq. (D4), (D7) and (D11) and setting the coefficients of the terms which oscillate at the same spatial and temporal frequencies equal to each other, we obtain an equation of motion for Ψ m j Ψ m j = i (mω + (f (l)q + k )v screw )δ m,l δ j,j + 2 sgn(γ)(−1) j + iα 1 + α 2F The Floquet matrix M F is non-Hermitian. Its complex eigenvalues describe the energy and decay rate of the the spin wave excitations on top of the Archimedean screw solution.

Appendix E: Onset of Chaos
As already mentioned in the main text and visualized in Fig. 7, we find a transition to chaotic behavior at some critical strength of the driving field in our simulations, as can be expected from a driven nonlinear system with many degrees of freedom. In order to further examine this transition, we analyze the dynamics of a single spin in more detail. First, we determine Ω screw from a linear fit to its azimuthal angle φ(t), see also Fig. 7. Then, we evaluate the orientation of the spin stroboscopically at times t n = 2πn/(ω − Ω screw ), and rotate the result by an angle −Ω screw t n around q e z , to eliminate the screw rotation. In Fig. 10 the projection of this spin onto the xy plane is shown for a range of driving field amplitudes. For weak driving, B ⊥ 0.6 mT (panels (a) and (b)), we obtain (within the numerical accuracy) a single point, which is the signature of the Archimedean screw phase. For stronger driving, the time quasicrystal forms as discussed in Sec. V. In this case the spin obtains an extra periodic oscillation, see also Fig. 7. Within our stroboscopic projection, this manifests in closed orbits visible in panels (c)-(o). For B ⊥ 4 mT, panel (p), chaos sets in. In Fig. 10, this manifests in aperiodic trajectories that fill certain regions of the plot. For stronger driving a larger area is filled, see panels (q)-(r).
We would like to emphasize that both the onset of chaos and the nature of the chaotic trajectories depends on the size of the unit cell used in our simulations (in Fig. 10 we use 15 2π q ). Smaller unit cells suppress chaos as they contain fewer degrees of freedom. Also in the chaotic regime, we expect translational invariance in the direction perpendicular to the helix not to be valid anymore.