Topological frequency conversion in Weyl semimetals

We theoretically predict a new working principle for optical amplification, based on Weyl semimetals: when a Weyl semimetal is suitably irradiated at two frequencies, electrons close to the Weyl points convert energy between the frequencies through the mechanism of topological frequency conversion from [Martin et al, PRX 7 041008 (2017)]. Each electron converts energy at a quantized rate given by an integer multiple of Planck's constant multiplied by the product of the two frequencies. In simulations, we show that optimal, but feasible band structures can support topological frequency conversion in the"THz gap"at intensities down to $ 2{\rm W}/{\rm mm^2}$; the gain from the effect can exceed the dissipative loss when the frequencies are larger than the relaxation time of the system. Topological frequency conversion provides a new paradigm for optical amplification, and further extends Weyl semimetals' promise for technological applications.

We theoretically predict a new working principle for optical amplification, based on Weyl semimetals: when a Weyl semimetal is suitably irradiated at two frequencies, electrons close to the Weyl points convert energy between the frequencies through the mechanism of topological frequency conversion from [Martin et al, PRX 7 041008 (2017)]. Each electron converts energy at a quantized rate given by an integer multiple of Planck's constant multiplied by the product of the two frequencies.
In simulations, we show that optimal, but feasible band structures can support topological frequency conversion in the "THz gap" at intensities down to 2W/mm 2 ; the gain from the effect can exceed the dissipative loss when the frequencies are larger than the relaxation time of the system. Topological frequency conversion provides a new paradigm for optical amplification, and further extends Weyl semimetals' promise for technological applications.
While the model from Ref. [52] has been experimentally implemented and studied [55,56], actual observation of topological frequency conversion is still lacking. The reasons are two fold: First, in the magnetic realm, topological frequency conversion in the desirable frequency regime of THz and above requires extremely high amplitudes of the oscillating magnetic field (of about 1 Tesla and above, corresponding to radiation intensities of more than 240 MW/mm 2 ). Even then, measurable -and especially useful -conversion rates would require many spins acting synchronously.
In this work we propose a Weyl semimetal as the medium of choice for realizing topological frequency con- Near εF = 4 meV, these values correspond to radiation intensities of 0.5 W/mm 2 , and 2 W/mm 2 , respectively. Blue, orange and green data result from relaxation times 200, 400, and 600 ps, respectively. (c) Energy transfer per unit volume, as a function of relaxation time τ , for an isolated Weyl node with Fermi velocity 3.87 · 10 5 m/s. Modes 1 and 2 have amplitudes 900 kV/m and 1800 kV/m inside the material, respectively. Mode 2 has frequency f2 = 1.23 THz, and mode 1 has frequency f1 = ( √ 5 − 1)f2/2, (blue), f1 = 2f2/3 (green), and f1 = 2f2/(3 + 0.001π) (orange). See Secs. VI D and IV for further details of the calculation used for panels (b) and (c), respectively.
version at high frequencies and with large conversion rates. For that we consider a Weyl semimetal, subjected to incoming radiation at two incommensurate fre-arXiv:2201.07804v3 [cond-mat.mes-hall] 15 Jul 2022 quencies, as depicted in Fig. 1(a). Under appropriate driving, individual electrons near the Weyl nodes act as an ensemble of topological frequency converters (as in Ref. 52), with the (pseudo-)spin of each electron playing the role of the spin, and the vector potential potential inside the material playing the role of the magnetic field (the "transduction" being provided by the Fermi velocity of the Weyl point). As a result, the system hosts an ensemble of electrons that each convert energy from mode 2 to mode 1 at the quantized rate ±hf 1 f 2 per electron; the number of active frequency converters is controlled by the magnitude of the vector potential. Importantly, topological frequency conversion can be realized in Weyl semimetals at relatively modest radiation intensities. This is because the effective spins interact directly with the (strongly coupled) electric field of the radiation rather the than the magnetic field. As another benefit, Weyl semimetals host a macroscopic number of active frequency converters, giving rise to very large conversion rates.
As a bulk response, topological frequency conversion is unique to Weyl semimetals, and constitutes a fundamentally new mechanism for optical amplification. The phenomenon has novel features of intrinsic interest: first, it is a 2-wave mixing effect that does not require an idler beam or phase matching. Secondly, it is in essence a nonperturbative effect, beyond the regime of standard "χ n " responses: in the ideal, fully adiabatic, limit, we show that the rate of topological frequency conversion is non-analytic as a function of the driving amplitude, and hence cannot be captured through a standard Taylor expansion. Away from this limit (i.e., in the presence of finite driving frequency and relaxation), the nonperturbativeness persists in the form of a highly nonlinear amplitude-dependence.
The novel features above, along with the modest radiation intensities required and the macroscopic number of active frequency converters give Weyl semimetals a significant potential for optical amplification. This is demonstrated in Fig. 1(b): here we plot the gain coefficient (i.e., the exponential rate at which the intensity of the amplified mode increases inside the material), obtained from simulations with a somewhat optimized, but feasible band structure of a Weyl semimetal. The material is irradiated at frequencies in the "THz gap," where new effective amplifiers are in high demand, due to a lack of powerful coherent radiation sources. Assuming sufficiently slow relaxation, our simulations indicate gain coefficients of order 100 cm −1 can be achieved at intensities of order 1 W/m 2 . This value is comparable with current methods such quantum cascade lasers [57][58][59][60][61], which report gain coefficients, 20−50cm −1 range [57,58]. We emphasize it may be possible to realize significantly larger gain coefficients than O(100 cm −1 ) in other parameter ranges; e.g., with stronger intensities.
There still are challenges that need to be overcome before optical amplification can become reality: being a conductor, Weyl semimetal respond with plasma oscilla-tions to radiation which renormalize the vector potential inside the material. It is therefore necessary to drive the system above its plasma frequency to allow the vector potential enter the material. The plasma oscillations on the other hand provides an opportunity: driving the material close to resonance with the plasma frequency amplifies the internal vector potential, thus significantly enhancing the rate of energy conversion. Indeed, we exploit this resonance effect to achieve the simulated gain coefficients of ∼ 100 cm −1 for the data depicted in Fig. 1(b).
Another, more serious, challenge is electronic relaxation processes. These counteract the frequency conversion by providing a channel for trivial energy dissipation -material heating.
For the parameters considered in Fig. 1(b), net energy gain of the pumped mode becomes possible for a characteristic relaxation time of order 300 picosecond at THz frequencies. Such relaxation times are longer than the relaxation times that have been mostly reported experimentally to date, which range from 0.25 ps−3 ps [62][63][64][65] to 40 ps [66]. The nature and timescales for scattering processes in Weyl semimetals is an interesting subject on its own which is still being explored, however: some experiments report signatures with much longer lifetimes [17,67,68] that can even exceed 1000 ps [69]. In addition experiments and theoretical studies indicate regimes dominated by non-standard, momentum-conserving channels of dissipation, resulting in hydrodynamical behavior [70,71].
We speculate that slower relaxation rates can be achieved, e.g. through improvement of materials quality and bath/substrate engineering. As another example, we show that dissipation is significantly reduced at commensurate frequencies, without affecting the energy transfer from topological frequency conversion [see Fig. 1(c)]. Excessive heating can be countered through pulsed driving, by allowing the system to dissipate away heat between the pulses. If sufficiently slow relaxation can be reached through such or similar incremental improvements, there is a potential for significant benefits in the form of a new and powerful mechanism for optical amplifcation.
The rest of this paper is structured as follows: in Sec. I we review the characteristic properties of Weyl semimetals, which forms the basis for our discussion. In Sec. II, we present the mechanism for frequency conversion from a single-particle perspective. Sec. III shows how topological frequency conversion arises in a realistic many-body system, taking into account the effects of finite frequency and dissipation. In Sec. IV, we support our conclusions with numerical simulations. In Sec. V, we summarize the conditions that a Weyl semimetal and driving modes must satisfy to allow for topological frequency conversion. In Sec. VI, we incorporate the effects of plasmons on the single-grain frequency converter, calculate the work in the context of Maxwell equations for the problem, and propose a practical implementation of an amplifier based on this effect using a "phase array" of Weyl grains. We conclude with a general discussion in Sec. VII. Details of derivations are provided in Appendices.

I. REVIEW OF WEYL SEMIMETALS
We begin by reviewing the characteristic properties of Weyl semimetals. This review forms the basis for our subsequent discussion.
Weyl semimetals are 3-dimensional materials in which two adjacent energy bands touch at isolated points in the Brillouin zone [11,12], as depicted in Fig. 2(a). These band-touching points are known as Weyl points. To understand Weyl points better, we consider the Bloch Hamiltonian of the system, H(k), near one such Weyl point, which we (without loss of generality) take to be located at wave vector k = 0. When restricted to the subspace spanned by the two touching bands, and linearized in k around k = 0, H(k) takes the following characteristic form: (1) where σ x , σ y and σ z denote the Pauli matrices acting on the subspace spanned by the two touching bands in some given basis, R is a real-valued symmetric full-rank 3 × 3 matrix, while V = (V 1 , V 2 , V 3 ) and ε 0 is a realvalued velocity and energy, respectively. Evidently, the two energy bands of H(k) included above touch at the Weyl point (k = 0). When the touching energy bands are plotted in the plane k i = 0 (for i = x, y, or z), the bands form a characteristic "touching cones" structure, as for example in Fig. 2(a). ε 0 determines the location of the touching point on the energy axis, while V determines the "tilt" of the cones. The eigenvectors and spectrum of R determines the anisotropy (or "squeezing") of the band gap around the Weyl point.
Once it is present, a Weyl point is a very robust feature: as long R remains full-rank, any infinitesimal perturbation to the system can only shift the location of the band-touching point, but not eliminate it. This is straightforward to verify through direct calculation. Hence, a smooth change of system parameters can only cause Weyl points to continuously move around in the Brillouin zone [72]. As a result of this robustness, Weyl semimetals are a generic class of materials. Indeed, many materials have recently been shown to be Weyl semimetals [10][11][12][73][74][75][76][77][78][79].
Another novel feature of Weyl semimetals is the nontrivial band topology associated with the eigenstates of the Bloch Hamiltonian, {|ψ α (k) }. These topological properties are captured by the Berry curvature with ijk denoting the Levi-Civita tensor and ∂ i the partial derivative with respect to the ith component of crystal wave vector, k i (we discuss the physical significance of the Berry curvature below). Weyl points act as point sources for Berry curvature: for two bands, 1 and 2, touching at an isolated Weyl node at k = k i , the Berry curvature of the upper band, 2, satisfies were | · | denotes the determinant, and ∇ the nabla operator in k-space. The sign is reversed for the the lower band. The relationship between Weyl points and Berry curvature is in exact analogy to point charges and the electric field. In this analogy, the index q ≡ sgn|R| determines the "charge", or chirality, of the Weyl point [80].
The net charge of all Weyl points that appear within a given gap is zero [1]; thus any gap must hold an even number of Weyl points.
For a system with many bands and multiple Weyl points, Eq. (3) generalizes to where the sum runs over all Weyl points in the system, q i denotes the chirality of Weyl point i, and s i,α indicates how the Weyl points of the system connect the bands: specifically s i,α = 1 if Weyl point i connects band α with the adjacent band above, s i,α = −1 if it connects band α with the band below, and s i,α = 0 if band α is not involved at Weyl point i.
Eq. (4) can equivalently be expressed using the divergence theorem: for a closed surface in the Brillouin zone, C, the total Berry flux of band α, C d 2 S · Ω α (k) (which is identical to the Chern number of band α when constricted to the 2-dimensional closed surface C), is given by ki∈C q i s α,i where the sum runs over all Weyl points contained within C.
Berry curvature acts as a magnetic field in reciprocal space: an electron in band α with a relatively well-defined position and wavevector, r and k, acquires a transverse velocity proportional to Ω α (k) when subject to a weak external force [81],k: This second term above is known as anomalous velocity, and can be seen as a canonically-conjugate analog to the Lorentz force: whereas a magnetic field B generates a velocity in reciprocal space perpendicular to the real-space velocity,k = − e B ×ṙ (the Lorentz force), Berry curvature generates a real-space velocity perpendicular to the reciprocal space velocity,ṙ =k × Ω n (the anomalous velocity); here e denotes the elementary charge. Eq. (3) implies that the Berry curvature diverges near Weyl points. Hence electrons with wavevectors near a Weyl point experience a divergent anomalous velocity [82]. When subject to an applied electric field, E, such thatk = −eE/ , Weyl semimetals can thus produce a large current response which may be nonlinear as a function of E. This significant nonlinearity makes Weyl materials particularly attractive as nonlinear optical media,  ), and with E1 = 720 kV/m, f1 = 2 3 f2 (resulting in a commenusrate frequency ratio). The different value of E1 is chosen to ensure that the vector potential of mode 1 has the same amplitude in panels (b) and (d), such that the topological phase boundary B0 is the same for panels (b-d).
In principle, any material the band geometry that has large local Berry curvature near the Fermi level is prone to having strong nonlinear response; Weyl semimetals are just a prominent example of those thanks to the divergent Berry curvature near the Weyl points. However, this is not the full story: the exotic band topology of Weyl semimetals (i.e. the nontrivial winding of the Berry curvature around Weyl points) in itself gives rise to unique nonlinear response phenomena. The effect we explore in this paper -topological frequency conversion -is an example of such an inherently topological response phenomenon.

II. FREQUENCY CONVERSION FROM A SINGLE ELECTRON
Here we show how the nontrivial band topology of Weyl semimetals allows electrons to act as topological frequency converters [52,83]. We consider a Weyl semimetal irradiated by two electromagnetic waves, or "modes", with distinct propagation angles and frequencies, and with elliptical or circular polarization. Fig. 1(a) depicts a concrete example in which the two waves are circularly polarized in the xz and yz planes. We let E 1 (t) and E 2 (t) denote the electric fields resulting from mode 1 and 2, respectively, such that the net electric field in the Weyl semimetal at time t is given by E(t) = E 1 (t)+E 2 (t). We assume the wavelengths of the incoming waves to be much longer than the relevant length scales we consider, and hence take E i (t) to be spatially uniform. The two modes are oscillating with frequencies f 1 and f 2 , such that, For simplicity, we first consider the case where f 1 and f 2 are incommensurate; we consider the case of commensurate frequencies in Sec. II A.
The coupling between the Weyl semimetal and the electromagnetic radiation is captured by the Peierls substitution [84], which causes the driven system to be governed by the time-dependent Bloch Hamiltonian where A(t) = − t 0 ds E(s) denotes vector potential induced by E(t) [85]. In the following H(k) denotes the Hamiltonian of the system in the absence of the driving, while H(k, t) denotes the Hamiltonian in the presence of the driving.
It is useful to decompose the vector potential as is generated by electromagnetic radiation, its time-average vanishes; hence A i (t) is also T i -periodic with respect to t. Without loss of generality we take both A 1 (t) and A 2 (t) to have time-average zero (recall that constant shifts in A(t) correspond to benign gauge transformations). It is convenient to represent the vector potentials A and A i as explicit functions of the phases of the two modes , α and α i (rather than single time variable).
where ω i ≡ 2πf i denotes the angular frequency of mode i. We similarly let (φ 1 , φ 2 ) and i (φ i ) denote the electric fields E and E i as functions of the individual phases the two modes.
To reveal how topological frequency conversion emerges, we consider the dynamics of a single electron in band α, in a wavepacket with some relatively well-defined position, r, and wavevector, k. The rate of energy transferred to mode 1 by the wavepacket, P α (k, t), is given by Ohm's law [86], whereṙ α (k, t) denotes the velocity of the wavepacket in band α at wavevector k, given the Hamiltonian H(k, t).
When ω 1 and ω 2 are small enough so that the timedependence of H(k, t) is adiabatic [87],ṙ α (k, t) is given by Eq. (5), with the instantaneous reciprocal space velocity given byk(t) = −eE(t)/ : Our goal is to compute the time-averaged rate of energy transfer into mode 1, Here and throughout this work, we use the· accent to indicate time-averaging, such that, for any function of time and, possibly, other parameters Here v α (k; φ 1 , φ 2 ) is obtained from the expression forṙ α (k, t) in Eq. (8) after replacing A(t) and E(t) with α(φ 1 , φ 2 ) and (φ 1 , φ 2 ), respectively. Since we assume ω 1 and ω 2 to be incommensurate, the time-averaged value of E 1 (t) ·ṙ α (k, t) is identical to the phase-averaged value of 1 (φ 1 ) · v α (k; φ 1 , φ 2 ). Hence, Using the expression for v we described above, along with The integral above has a direct geometrical interpretation: gives the differential area element of the closed surface defined by eα(φ 1 , φ 2 )/ in reciprocal space, The direction of the differential area element (∂ φ1 α × ∂ φ2 α) defines the orientation of B 0 . In Fig. 2(b) we depict B 0 for the case where modes 1 and 2 are circularly polarized in the xz and yz planes respectively, and have electric field amplitudes E 2 = 1000 kV/m, E 1 = 740kV/m, and frequencies f 2 = 1 THz, f 1 = √ 5−1 2 f 2 . For incommensurate frequencies, the trajectory of eA(t)/ fills out B 0 completely at long times, as also illustrated in Fig. 2 With the above geometric interpretation, we find where B0 d 2 k denotes the surface integral of k over the surface B 0 . From Sec. I we recall that this integral is quantized as 2π times the net charge of Weyl points of band α enclosed within the surface B 0 after displacing it by k from the origin in reciprocal space, Q α [k] (here the enclosed charge is weighted by the orientation of B 0 with respect to the volume in which the Weyl point is enclosed):P where we used h = 2π .
For an isolated Weyl point with charge +1 located at k = 0 in a two-band system, Q α [k] is given by the following for the upper band (α = 2): where the function W (k) is integer-valued and denotes the net winding number of B 0 around k as a function of φ 1 and φ 2 . In Fig. 2(c) we plot W (k) for the configuration of two circularly polarized modes also considered in Fig. 2(b). The sign is reversed for mode 2 [i.e., when replacing E 1 with E 2 in Eq. (7)]. Hence the electron acts as a conversion medium that transfers energy between mode 2 and 1.
For a system with multiple bands, Q α [k] = i W (k − k i )q i s i,α , where the index s i,α encodes how Weyl point i connects the bands of the system (see Sec. I). We hence arrive atP This constitutes one of our main results.
Eq. (16) shows that each electron in the Weyl semimetal transfers energy from mode 2 to mode 1 at a rate which is quantized, as an integer multiple of hf 1 f 2 . The value of the integer depends on the location of the electron in the Brillouin zone, k. Specifically, the conversion rateP α (k) is nonzero for electrons whose wavevectors k are located within the surface B 0 relative to a Weyl point. Thus, a nonzero conversion power can be realized for electrons near Weyl points.
The energy conversion predicted in Eq. (16) can be seen as a realization of the topological frequency conversion that was discovered in Ref. [52]. Ref. [52] showed that a 2-level system (such as a spin-1/2) initialized in its lower band and adiabatically driven by two modes with frequencies f 1 and f 2 can transfer energy between the modes at an average rate quantized as hf 1 f 2 z, where z is an integer. Ref. [52] explained this conversion as an anomalous velocity along the synthetic dimensions that correspond to the photon numbers of the two modes. To understand the relationship between our result and Ref. [52], note that for fixed k, H(k, t) is a Hamiltonian of a 2-level system of the exact same form as considered Ref. [52], with the pseudospin of the electron playing the role of the physical spin in Ref. [52]. Indeed, the arguments of Ref. [52] show that for the two-level system described by H(k, t), z = W (k). In this way, each electron in a Weyl semimetal can be seen as a topological frequency converter from Ref. [52], with the quantized rate of conversion controlled by its location in the Brillouin zone.

A. Commensurate frequencies
The discussion above for simplicity assumed the frequencies f 1 and f 2 incommensurate. Here we consider the case where the frequencies of the modes are commensurate such that f 1 /f 2 = p/q for some integers p and q. In this case, E(t) and A(t) thus are time-periodic with the extended period T ext = pT 1 = qT 2 . This time-periodicity significantly affects the electron's trajectory in the BZ (relative to its equilibrium wavevector), eA(t)/ . For incommensurate frequencies, the trajectory fills a closed surface, namely B 0 , as illustrated in Fig. 2(b). In contrast, commensurate frequencies causes the trajectory to form a closed curve, C 0 , as in Fig. 2 For commensurate frequencies, the driving experienced by the electron depends on the initial phase difference between the modes, ∆φ; here nonzero ∆φ corresponds to a shift of the phase of mode 2 such that E(t) = E 1 (t) + E 2 (t + ∆φ/ω 2 ), resulting in E(t) = (ω 1 t, ω 2 t + ∆φ) and A(t) = α(ω 1 t, ω 2 t + ∆φ). For incommensurate frequencies, different values of ∆φ are equivalent to shifts in the time origin and hence do not affect the long-term dynamics of the electron. In contrast, for commensurate frequencies, each distinct value of ∆φ results in a different closed trajectory of A(t), C 0 . The surface B 0 is recovered by combining the curves C 0 for all possible values of ∆φ.
For commensurate frequencies, the quantization ofP breaks down. The breakdown of quantization arises because the trajectory of the modes' phases (φ 1 (t), φ 2 (t)) = (ω 1 t, ω 2 t + ∆φ) does not cover the whole 2d phase Brillouin zone over time, φ 1 , φ 2 = [0, 2π), thus invalidating the step leading to Eq. (10). However, quantization is recovered when averagingP over all possible values of ∆φ: for commensurate frequencies, Eq. (10) remains valid for the average value ofP with respect to ∆φ. Thus, for commensurate frequenices it is possible to enhance conversion rates relative to the quantized value by tuning the phase difference to a value where the conversion rate exceeds its average value. For uncontrolled (random) phase differences, the conversion rate remains quantized on average.

III. FREQUENCY CONVERSION IN MANY-BODY SYSTEMS
Our next goal is to show how topological frequency conversion emerges in a realistic Weyl semimetal where electrons are affected by interactions, impurities and phonons. We focus on the rate of energy transfer to mode 1 per unit volume for a Weyl semimetal driven by two modes, η(t). If η(t) is positive, there is a net flow of energy into mode 1, implying amplification of this mode. This energy must originate from mode 2. The conversion rate η(t) can be computed from the current density, j(t), using Ohms law: To obtain the current density j(t) we characterize the many-body state of the Weyl semimetal in terms of the momentum resolved density matrix, whereρ F (t) denotes the full density matrix of the Weyl semimetal at time t, which is subject to interactions, impurities, and phonons, while Tr k =k [·] denotes the trace over all possible occupations of electronic states with crystal momentum other than k.ρ(k, t) is a matrix in the 2 d dimensional Fock space associated with the d orbitals (or bands) accessible by the electrons at wavevector k [88] Below, the "hat" accent· indicates operators that act on many-body orbital Fock states. Operators without the accent, such as the Bloch Hamiltonian from Secs. I-II, H(k, t), are single-particle operators.ρ(k, t) encodes the band occupancies alongside with inter-band coherences and all multi-particle correlations of electrons with the same wavevector k. The inter-band coherences are crucial for capturing topological energy conversion, since they give rise to the anomalous velocity in our formalism. ρ(k, t) determines the current density in the system, j(t), through where momentum integrals are taken over the full Brillouin zone, andĤ(k, t) denotes the second-quantized Bloch Hamiltonian of the system: Here H ij (k, t) ≡ i|H(k, t)|j , and |i denotes the ith orbital state in the standard Bloch space.
In the presence of drivingρ(k, t) approaches to a timedependent steady state. We obtain this steady state by solving a master equation forρ(k, t), in which the effects of interactions, impurities and disorder are included as a dissipative term. The master equation and steady-state solution are summarized in Sec. III A below. The calculation of the steady-state is straightforward, but involved, and is detailed in Appendix B. A key feature of the steady-state solution is that the current response can be split into an energy-conserving "adiabatic component", j 0 (t), and a dissipative correction due to nonadiabaticity and scattering, δj(t): This decomposition allows us to identify an energyconserving and dissipative component of η(t): The component j 0 (t) is responsible for topological frequency conversion, and we find that this term dominates in the limit of adiabatic driving and slow relaxation. As a central result, we find that In what follows, we first discuss the steady state solution of the density matrix (Sec. III A). We then compute the time-averaged energy pumping resulting from the non-dissipative component of the current response, η 0 . We finally consider the dissipative component of η, η dis in Sec. III C. It is crucial to estimate η dis , since amplification is only achieved when η 0 exceeds η dis .

A. Steady state solution
We now discuss how we obtain the steady-state of ρ(k, t). Details of this discussion are provided in Appendix B.
For a clean and non-interacting system,ρ(k, t) evolves according to the von Neumann equation, Interactions, phonons, and impurities cause a dissipative correction to this equation. For sufficiently weak dissipation, this correction can be derived approximately from first principles and takes the form of a trace-and positivity-preserving linear operator acting onρ(k, t), D(k, t) [89]. Thus,ρ(k, t) is governed by the following Lindblad-type quantum master equation In Appendix. B, we obtain a solution to Eq. (24). The solutionρ(k, t) is accurate as long as the driving is adiabatic with respect to the energy gap of H(k, t), δε, and much faster than the the magnitude of the dissipator D [90]: We term this limit, as the coherent adiabatic regime.
When the above conditions are satisfied, we find thesteady state value ofρ(k, t) is diagonal in the eigenbasis of the Hamiltonian,Ĥ(k, t), up to minor nonadiabatic corrections. The corresponding eigenvalues (which determine the the occupations of the instantaneous bands ofĤ(k, t)) are nearly stationary, except for minor fluctuations of order D /O(ω 1 , ω 2 ). These fluctuations, along with the subleading (second-order) nonadiabatic corrections toρ(k, t) give rise to the dissipative current δj(t). The term j 0 (t) results from just keeping the (dominating) time-independent component of the eigenvalues ofρ(k, t) and including leading-order nonadiabatic correction to its eigenbasis. Here the leading-order non-adiabatic correction to the eigenbasis is responsible for the anomalous velocity which enters in j 0 (t).
Our solution to Eq. (24) applies to any dissipator D, and this dissipator can be derived from first principles [89]. However, for illustrative purposes, we now demonstrate our solution for the concrete example where D takes a particular phenomenological form: the "Boltzmann" form. In this approximation, the dissipator uniformly relaxes electrons towards their instantaneous equilibrium state at some given ambient temperature 1/β and chemical potential µ, and with some rate 1/τ : Hereρ eq (k, t) is the instantaneous equilibrium state described above, and is given by We also use this dissipator in our numerical simulations (see Sec. IV).
The Boltzmann-form dissipator [Eq. (26)] leads to the following steady state density matrix: where f β (E) denotes the Fermi-Dirac distribution at temperature 1/β. This result (see also Appendix B) indicates a steady state occupation which is the average band-population on the trajectory k + eA(t)/ , as if the equilibrium distribution is "smeared" over a characteristic wavevector scale eA/ , where A is the drive vector potential magnitude. This smearing is confirmed in our numerical simulations (see Sec. IV and Fig. 4 in particular).

B. Non-dissipative frequency conversion
We first compute the average rate of energy transfer in the limit of adiabatic driving and zero dissipation. I.e., we compute the the time-average of the component η 0 (t), η 0 . We find that η 0 (t) can have nonzero time-average because of the mechanism of topological frequency conversion that we discovered in the last section.
To computeη 0 we first note η 0 can be written where P α (k, t) ≡ eE(t) ·ṙ α (k, t) [see also Eq. (7)]. We find the time-average of the above using the main result from Sec. (16)]. Here q i , s i,α , and k i denote the charge, band connectivity, and wave vector of Weyl point i in the system, respectively, while W (k) measures the net winding of the surface B 0 around wavevector k (see Sec. III C for further details). with this, Eq. (28) becomes: Thus, each Weyl point is surrounded by a region of reciprocal space (namely the region where W (k − k i ) = 0), in which electrons act as topological frequency converters.
In this region, each transfers energy to mode 1 at the quantized rate ±hf 1 f 2 . This is a many-electron generalization of Eq. (16) and constitutes another main results of this paper. In the following we thus refer toη 0 as the topological frequency conversion rate of the system, to distinguish it from the dissipation rate, which is given by the time-average of η dis (t).
While the conversion rate from each electron is quantized, the net number of electrons with nonzero conversion rate is not fixed, but depends on the amplitude and configuration of the driving field (through the function W (k)) and the steady-state distribution surrounding each Weyl point,ρ α (k). This steady-state distribution is in turn controlled by the band structure of the system, as well as the configuration and intensity of the external driving.
To explore how the band structure and driving configuration controls the conversion rate, we first estimate the "gross" rate of topological frequency conversion from a Weyl point (i.e., not taking into account cancellation between electrons that transfer energy at opposite rates). Note that W (k) is positive within volume of As an example, for ω i ∼ 2π THz and E i ∼ 1500 kV/m. the above estimate yields η gross ∼ 500 kW/mm 3 . The actual, net, topological conversion power,η 0 is significantly smaller than the gross rate we estimated above, due to cancellation between electrons that convert energy at rates hf 1 f 2 and −hf 1 f 2 . Specifically, when modes 1 and 2 only contain a single harmonic each, the driving [91] [this symmetry is clearly evident in Fig. 2(c)]. Hence the regions of reciprocal space characterized by conversion rates hf 1 f 2 and −hf 1 f 2 have equal net volumes. In realistic situations, both volumes will be occupied by electrons, implyingη 0 η gross . However, because η gross can be quite large, even a small imbalance in the filling of the two regions can lead to significant net frequency conversion.
To achieve a nonzeroη 0 , the steady-state occupation of the bands,ρ α (k), must be anisotropic around the Weyl point to counteract the antisymmetry W (k) = −W (−k). Such an anisotropy is generally achieved when the "Weyl cone tilt", V, is nonzero, since we expect the steady-state inherits the same symmetry properties as the equilibrium state [see Eq. (1)]. Additionally, the Weyl point must be within a distance of order ∼ eA i / from the Fermi surface to ensure thatρ α (k) does not take constant value (1 or 0) within B 0 . Indeed, our numerical simulations demonstrate that nonzeroη 0 can arise when V = 0 and the Fermi surface lies close to the Weyl point.
Topological frequency conversion is in essence a nonperturbative effect: it is controlled by the overlap of the quantized (i.e., nonanalytic) function W (k) with the steady-state distribution. Hence topological frequency conversion does not have a simple power-law dependence on A in the limit of small A, and is therefore beyond standard nonlinear response theory. In Sec. IV [ Fig. 6(b)] we provide data from numerical simulations indicating this highly nonlinear nature of the phenomenon.

C. Dissipative energy loss
For topological frequency conversion to cause a net amplification of mode 1,η 0 must exceed the rate of energy loss due to dissipation, η dis . It is therefore crucial to estimate this dissipation rate. This is the goal of this subsection.
Our solution of the master equation in Appendix B shows that the dissipative current response, δj(t) [see Eq. (21)], contains two components: which we interpret as arising from momentum-relaxation (δj mr ) and nonadiabaticity-induced particle-hole pair creation (δj na ); see discussion below. Consequently, η dis (t) can also be separated into these two components: where η mr (t) ≡ E 1 (t) · j mr (t), and η na is defined likewise. While η mr and η na are given in Appendix B, here we discuss their origin and estimate their magnitudes based on a phenomenological discussion.
Energy loss due to momentum relaxation, η mr , arises when perturbed electrons in the close vicinity the Fermi surface relax due to their displacement from instantaneous equilibrium, as schematically indicated in Fig. 3(a). In contrast, η na (t) arises from the particle-hole pair creation that results because the driving-induced electric field inevitably overpowers the gap sufficiently close to Weyl points. Equivalently, η na arises because effective gap closing ofĤ(k, t) near a Weyl point gives rise to Landau-Zener tunneling from the conduction to the valence band upon driving. These excited electrons dissipate energy as they relax back to the conduction band, as schematically illustrated in Fig. 3 Below we estimate η mr and η na based on phenomenological arguments. For simplicity, we consider system with two bands and an isolated Weyl point at k = 0 (which is easily generalized to multiple Weyl points). Moreover, we do not distinguish between the characteristic relaxation rates associated with momentum relaxation and particle-hole pair creation (which may be different in real materials), but use as an estimate for both characteristic relaxation rates.

Momentum relaxation
We first estimate the rate of energy loss arising from momentum relaxation, η mr (t). For convenience, in the following we let A = E/Ω denote the characterstic magnitude of the driving-induced vector potential, E the characteristic magnitude of E(t), and Ω the characteristic scale of ω 1 and ω 2 .
Effects of momentum relaxation can only emerge within a distance ∼ eA/ from the Fermi surface, where electronic occupation fluctuates. Therefore, only a density of eS F A/(2π) 3 contributes to δj mr (t), where S F is the area of the Fermi surface. Electrons near the Fermi surface on average gain an energy of order eAv F /2 due to the driving (with v F the characteristic Fermi velocity). Assuming their relaxation rate is of order Γ, the average rate of energy loss in the system due to momentum relaxation thus is given by . We estimate that half of this comes from mode 1. This estimate agrees well with our predictions based on the definition of δj mr (t) in Appendix B. Using that A ∼ E/Ω, we find It is interesting to note that η mr is proportional to the size of the Fermi surface, S F . Therefore type-I Weyl points are more suitable for frequency conversion than type-II Weyl semimetals, since the surrounding Fermi surface forms a compact ellipsoid for the former case, in contrast to an extended hyperboloid in the latter.
As an example, we estimate η mr for the same parameters we used to estimate η gross in Sec. III B i.e., E ∼ 1500 kV/m, v F ∼ 5 · 10 5 m/s, Ω = 2π THz. For an ellipsoid Fermi surface with principal semi axes (1.5, 1.5, 2.4)eA/ (yielding S F ≈ 0.06Å −2 ), our estimate then results in η mr ∼ 5 · 10 −5 J/mm 3 τ , where τ = 1/Γ. Note that our estimated value of η mr is proportional to the Fermi surface area and thus can be easily adjusted to other values of this quantity. Recalling that η gross ∼ 500 kW/mm 3 for the same parameters, we expect net frequency conversion to only exceed momentum relaxation when τ 100 ps. This expectation is confirmed in our numerical simulations [see Fig. 1(b)].

Non-adiabatic heating
Non-adiabatic heating arises from electrons at k-points where the time-dependence of H(k, t) = H(k + eA(t)/ ) overwhelms the band gap. In Weyl semimetals, such kpoints inevitably exist (even for arbitrarily slow driving), because the band gap of H(k) closes at Weyl points. For each such k-point, the band gap of H(k, t) effectively closes at certain times t, namely when k + eA(t)/ is sufficiently close to the Weyl point at k = 0 (see below for more detailed conditions). At each such gap-closing event, electrons at wavevector k will undergo partial or complete Landau-Zener transition from the conduction to the valence band. This mechanism effectively heats the electrons and eventually results in dissipation once the excited electrons relax. The dissipation induced by the mechanism above is captured by η na (t).
To estimate η na we first identify the set of k-points for which the time-dependence of H(k, t) is non-adiabatic; we term this region of reciprocal space as the "nonadiabatic" region and denote it by V na . The Landau-Zener formula [92,93] states that time-dependence of H(k, t) is non-adiabatic if, for some t, where δε(k) ≡ ε 2 (k) − ε 1 (k), and ε α (k) denotes the αth energy band of H(k). Using the linearized form of H(k) in Eq. (1), a straightforward derivation (see Appendix C) shows that this condition is satisfied at k-points for which while R and v 0 denotes the largest and smallest eigenvalue of the velocity matrix R, respectively [see Eq.
(1)] [94]. For incommensurate frequencies, V na thus consists of all k-points within a distance d 0 from the topological phase boundary B 0 by our estimate. For commensurate frequencies V na consists all k-points within a distance d 0 from C 0 , which forms a closed curve on B 0 , as in Fig. 2(d).
Electrons with wave vectors k within V na encounter a vanishing gap ofĤ(k, t) at times t where |k + eA(t)/ | d 0 . These electrons then undergo Landau-Zener tunneling, which effectively heats them to a high-temperature state, as explained in the beginning of this subsection. These high-temperature electrons then relax back to equilibrium after a characteristic time 1/Γ. η na (t) is then the rate of energy loss, or heating, (per unit volume) arising from this relaxation. We estimate where n na is the concentration of excited electrons within V na , and ∆ε na denotes the characteristic average value of δε(k+eA(t)) for k within V na . Here the factor of 2 comes because we estimate that the other half of the dissipated energy comes from mode 2. We obtain ∆ε na using that V na is located a distance ∼ eA/ from the Weyl node, such that ∆ε na ∼ eA R , and R is the largest velocity implied by the velocity tensor R. Using A ∼ E/Ω, we obtain To estimate n na , it is crucial to know the characteristic time interval between successive gap-closing events experienced by electrons with a given wavevector within V na , ∆t. To build intuition, let us first consider what happens when ∆t 1/Γ, i.e., when electrons have time to fully relax between successive gap-closing events [95]. Electrons at wavevector k are taken to a high-temperature state whenever k comes within a sphere of radius ∼ d 0 from eA(t)/ . Electrons are in equilibrium as they "enter" the sphere (due to our assumption ∆t 1/Γ), and we estimate that half of them are excited to the conduction band as they "leave" the sphere. The concentration of electrons per unit time that are heated by this process is hence given by the cross-section of this sphere times 0.5 (2π) 3 e|∂ t A|/ . Therefore, we expect the concentration of electrons heated per unit time to be given by ed 2 0 |∂ t A(t)|/16π 2 . Assuming the electrons relax with characteristic rate Γ, we estimate n na as the fixed point of ∂ t n na = eπd 2 0 |∂ t A(t)|/16π 2 − Γn na . Using ∂ t A = E, we thus find n na ∼ ed 2 0 E/16π 2 Γ. Next, we consider the case where ∆t 1/Γ. In this case, a significant fraction of electrons are already in a high-temperature state when they experience a gap-closing event (i.e., when they "enter" the sphere with radius d 0 centered at eA(t)/ ). Assuming that the gap-closing event effectively randomizes the state of the electrons (i.e., the electrons are in a infinitetemperature state right after "leaving" the sphere, regardless of their initial state), a subsequent gap-closing event only re-heats a reduced number of electrons to a high-temperature state. We estimate the fraction of preexcited electrons to be of order 0.5 e −Γ∆t right before the gap closing and 0.5 right after; thus the heating rate is reduced by a factor O(1 − e −Γ∆t ), resulting in Combining this result with Eqs. (36)- (38), we obtain Below we estimate ∆t (i.e., the characteristic time between gap closing events) for the two cases incommensurate and commensurate frequencies; as we find these two situations lead to significantly different ∆t, and hence also different values of η na . Evidently, the bound above is controlled by the ratio between the largest and smallest eigenvalues of the matrix R, R /v 0 . As we argued in Sec. I, this number quantifies the anisotropy of the band gap around the Weyl point.
Note that the first factor in Eq. (40) is larger than the gross rate of topological frequency conversion in Eq. (30) (this follows from R /v 0 ≥ 1, and π 2 > 8). Thus, ∆t needs to be much shorter than Γ −1 for nonadiabatic heating not to overwhelm the net rate of topological frequency conversion. In particular, since ∆t is at least 2π/Ω, we expect Γ Ω to be a necessary condition for topological frequency conversion.

Nonadiabatic heating at incommensurate frequencies
For incommensurate frequencies, we estimate ∆t as the time-window over which trajectory of eA(t)/ has length |B 0 |/d 0 , where |B 0 | denotes the area of the surface on which eA(t)/ is confined to at all times, B 0 ≡ {eα(φ 1 , φ 2 ), 0 < φ i < 2π} [see Sec. II and Fig. 2(b)] Since ∂ t A(t) = E(t), the trajectory of eA(t)/ over the time-window ∆t has length eE∆t/ . Estimating |B 0 | ∼ 4π 2 e 2 A 2 / 2 and using d 0 = eE R / v 2 0 along with A ∼ E/Ω, we hence obtain For the parameters we used to estimate η na above [E ∼ 1500 kV/m, v F ∼ 5 · 10 5 m/s, and Ω = 2π THz] this estimate yields ∆t ∼ 30 ps. To achieve topological frequency conversion at incommensurate frequencies with these parameters, the characteristic relaxation time τ = Γ −1 must be much longer than this timescale. In this limit (i.e., τ ∆t), we find η na ∼ 5 · 10 −8 kJ/mm 3 τ , which is smaller than η gross when τ 100 ps.

D. Lissajous Conversion
We now consider the case of commensurate frequencies, which can give a marked reduction of the non-adiabatic losses.
Eq. (40) shows that small time-intervals between subsequent gap-closing events, ∆t, leads to suppression of nonadiabatic heating, η na . An important consequence of this is that η na is strongly suppressed for commensurate frequencies, i.e., when for some integers p and q (which we, without loss of generality, take to have no common divisor). In this case, eA(t)/ forms a 3-dimensional Lissajous figure in reciprocal space, as in Fig. 2(d), and is time-periodic with period qT 1 = pT 2 . Thus, a given k-point within V na experiences gap-closing events with periodicity [96] ∆t = qT 1 = pT 2 .
As a consequence, η na is strongly suppressed for highly rational frequency ratios, i.e., when p and q are small. The suppression of η na means that net amplification from topological frequency conversion is significantly enhanced at highly commensurate frequencies. We term this mechanism of topological frequency conversion at commensurate frequencies Lissajous conversion. The dramatic suppression of nonadiabatic heating results in an enhanced net frequency conversion rate in the Lissajous conversion regime, as is evident in our numerical simulations (see Fig. 1(b), and Sec. IV). : net rate of energy transfer to mode 1 from electrons with a given wavevector k,P (k), as a function of k in the plane ky = 0, for same parameters as in Fig. 4(a). (c) occupation number in equilibrium of electronic modes with wavevector k, as a function of k in the plane ky = 0, for the same system as in panels (a-b). (d): net time-averaged occupation number of electronic modes with wavevector k, n(k) as a function of k in the plane ky = 0 for the system depicted in panel (a-c).
As an example, we consider Lissajous conversion at frequencies ω 1 = 2π THz ω 2 = 3 2 ω 1 . These parameters result in ∆t ∼ 3ps. To compare, recall that ∆t was estimated to 30 ps in the same frequency range. Using Eq. (37), we estimate the nonadiabatic heating rate to be given by η na ∼ 3.5 · 10 −9 kJ/mm 3 τ , which we expect can be smaller than η gross when τ 10 ps. In contrast, recall that our estimated nonadiabatic heating rate at incommensurate frequencies in the same frequency range is given by 5·10 −8 kJ/mm 3 τ , and thus more than 10 times larger.
In the limit of large p, q, we expect that our estimate for ∆t saturates at the expression we obtain for incommensurate frequencies.

IV. NUMERICAL SIMULATIONS OF FREQUENCY CONVERSION
We now support our theoretical predictions by data from numerical simulations.
In our simulations, we consider the dynamics of electrons near a single Weyl node in a Weyl semimetal with 2 bands. The electrons are subject to the linearized Bloch Hamiltonian We also introduce two electromagnetic modes that are circularly-polarizedthat propagate in the yz-and xzplanes, respectively. For i = 1, 2, mode i has angular frequency ω i and electric field amplitude E i inside the material. It thus induces the time-dependent electric field The irradiated electrons are governed by the timedependent Bloch Hamiltonian H(k, t) = H(k+eA(t)/ ), where A(t) denotes the driving-induced vector potential and is defined through ∂ t A(t) = E 1 (t) + E 2 (t) (see also Sec. II). As in the previous sections, we work in a gauge where A(t) has vanishing time-average. We numerically obtain the evolution of the momentum-resolved density matrix of the system, ρ(k, t) (see Sec. III for definition), using the master equation in Eq. (24).
We take the dissipator D to be given by the Boltzmann form [Eq. (26) Hereρ eq (k, t) denotes the instantaneous equilibrium state of electrons with crystal momentum k at time t at some given temperature T and chemical potential µ [see text below Eq. (26) for explicit definition]. Since Eq. (24) describes evolution in the 4-dimensional second-quantized Bloch space of the system, its numerical solution is relatively inexpensive.
For each k, we numerically solve Eq. (24) to obtain the steady-state evolution ofρ(k, t). From this steady state we extract the quantitȳ P (k) ≡ lim P (k) gives the time-averaged total rate of energy transfer to mode 1 from electrons with wavevector k. The total time-averaged rate of energy transferred to mode 1 per unit volume of the whole system,η, is obtained by integrating P (k) over all wavevectors: In our simulation, we evaluate the k-integral above by samplingP (k) over a large number of uniformly distributed values of k [97]. We solve the master equation forρ(k, t) through direct integration, not making use of any of the approximations of Sec. III C. In particular, our simulation does not distinguish between coherent and incoherent dynamics, and our obtained value forη thus includes both contributions both from topological frequency conversion and dissipation. Hence our simulation can be used to test the conclusions in Sec. III.
We probe different values of f 1 and τ , while keeping all other parameters fixed at values f 2 = 1.23 THz, v = 3.87 · 10 5 m/s, V = (0, 0, 3.1 · 10 5 m/s), µ = 115 meV, T = 20 K, E 1 = 0.9 MV/m, and E 2 = 1800 kV/m. Our chosen values of v and V have magnitudes comparable to those in real materials [98,99]. The values of µ and V are chosen to maximize the imbalance between the number of electrons acting as frequency converters at rates hf 1 f 2 and −hf 1 f 2 , as discussed in Sec. III B (see also Sec. IV B below).

A. Identification of amplification regime
In Sec. III, we showed that the time-averaged rate of energy transfer to mode 1 can be decomposed as η =η 0 +η dis . Hereη 0 can be positive due to topological frequency conversion, whileη dis is negative and measures the time-averaged rate of energy dissipated from mode 1 due to heating in the system. We expect |η dis | to decrease with increasing relaxation time τ , whileη 0 remains constant. Thus,η should increase with τ . There should also exist a critical value of τ for whichη = 0. When τ is larger than this "amplification threshold", the system will amplify mode 1 (η > 0). We expect the amplification threshold to be significantly lower in the Lissajous conversion regime (i.e., at rational frequency ratios) than for irrational frequency ratios due to the suppression of nonadiabatic heating in the former case (see Sec. IV C and Sec. IV C below).
To identify the amplification threshold for the system, we computedη as a function of τ for three representative choices of f 1 /f 2 ; namely, irrational f 1 = 1 ϕ f 2 , rational f 1 = 2 3 f 2 , and nearly-rational f 1 = 2 3+ f 2 , where ϕ is given by the "golden mean", 1 2 (1 + √ 5), and = π/1000. We keep all other parameters fixed at the values we specified earlier. The two latter values of f 1 are chosen to demonstrate the mechanism of Lissajous conversion: Whereas f 1 = 2 3 f 2 is commensurate with f 2 , f 1 = 2 3+ f 2 is not, and hence the former value of f 1 is expected to yield more efficient-Lissajous-conversion.
In Fig. 1(c) we plotη as a function of τ for the three values of f 1 above. As we expect,η increases as a function of τ for all choices of f 1 , and attains positive value for sufficiently large τ . For the irrational frequency ratio f 1 = f 2 /ϕ, the amplification threshold is reached at τ ≈ 1000 ps, for f 1 = 2f 2 /3 at τ ≈ 300 ps and for f 1 = 2f 2 /(3 + ) above τ = 1200 ps.
Note that the weak detuning of f 1 from 2f 2 /3 (green curve) to 2f 2 /(3 + ) (orange curve) reducesη by more than 100 kW/mm 3 , and pushes the amplification threshold from 300 to 1200 ps. This demonstrates the strong dependence of the net conversion rate on the commensurability of f 1 and f 2 that we discussed in Sec. III D.

B. Origin of energy conversion
Next, we confirm that the amplification of mode 1 (i.e, the positive values ofη > 0) we observed is due to topological frequency conversion. To this end, we computē P (k) as a function of k around the Weyl point.
We first review the signatures of topological frequency conversion we expect to see. For k-points where H(k, t) changes adiabatically in time, electrons should act as topological frequency converters (as in Ref. [52]) that transfer energy to mode 1 at an average rate quantized as ±hf 1 f 2 W (k), where the W (k) denotes the integervalued net winding number of the surface B 0 around k [see Fig. 2(c)]. Here + and − result from electrons in band 1 and 2, respectively. We hence expect whereρ α (k) denotes the time-averaged occupancy of band α, and P dis (k) denotes the rate of energy loss from mode 1 due to dissipation. We expect the latter is always negative, but only significant around the Fermi surface (due to momentum relaxation), and within the nonadiabatic region (due to nonadiabatic heating). In Fig. 4(a) we plot W (k) in the plane k y = 0, for f 1 = f 2 /ϕ and with all other parameters specified below Eq. (48). We also indicate the Fermi surface (dashed line) and schematically indicate the nonadiabatic region (shaded region), which surrounds the topological phase boundary (solid line). Since µ > 0, band 1 is fully occupied in equilibrium. We therefore expectρ 1 (k) ≈ 1 [see Eq. (27)] for all k away from the nonadiabatic region (where Landau-Zener tunneling can induce holes). We hence expect topological frequency conversion causes P (k) to approximately take value hf 1 f 2 [1−ρ 2 (k)] within the red region of Fig. 4(a), value −hf 1 f 2 [1 −ρ 2 (k)] in the blue region, and value 0 in the white region.
In Fig. 4(b), we plotP (k) in the plane k y = 0 in for parameters f 1 = f 2 /ϕ and τ = 51.6 ps [100]. The data shows clear signatures of topological frequency conversion, in the form of two "topological plateaux" of the Brillouin zone whereP (k) takes positive and negative values, respectively. These plateaus coincide closely with the regions in Fig. 4(a) where W (k) = ±1.P (k) approximately takes value hf 1 f 2 within the red plateau (away from the Weyl point), and value between 0 and −hf 1 f 2 within the blue plateau (close to the Weyl point).
In addition to topological frequency conversion, the data in Fig. 4(b) also shows clear signatures of the two distinct mechanisms for dissipation that we identified in Sec. III C, i.e., momentum relaxation and nonadiabatic heating:P (k) takes large negative values within the nonadiabatic region, as we expect from nonadiabatic heating, and moderate negative values around the Fermi surface, as we expect from momentum relaxation.
Note also that the data in Fig. 4(d) are in good agreement with our prediction that in the regime τ 1/Ω, the steady state band populations are effectively "smeared" versions of their equilibrium counterparts [see Eq. (27)]: The distribution in Fig. 4(d) clearly resembles a "smeared" version of the ellipsoid-profile that occurs in equilibrium [plotted in Fig. 4(c)].
Finally, the data in Figs. 4(cd) demonstrate how a nonzero value of the "Weyl cone tilt" V is needed to nonzero net rate of topological frequency conversion,η 0 , as we discussed in Sec. III B. The nonzero value of V, which causes an ellipsoid-profile of n(k) in equilibrium [ Fig. 4(c)], results in a "smeared ellipsoid" profile of n(k) in the steady state. As a result of this smeared ellipsoid profile, the region characterized by W (k) = 1 has a larger volume in which [ρ 2 (k) −ρ 1 (k)] > 0 ( n(k) ≤ 2) than than the volume where W (k) = −1, allowing for a nonzero value ofη 0 .

C. Lissajous Conversion
We finally verify that the enhancement ofη in the Lissajous regime (i.e., at commensurate frequencies) is due to the suppression of nonadiabatic heating. To this end, we plot in Fig. 5P (k) for the parameters τ = 516 ps and f 1 = 2f 2 /3 (a) and f 1 = 2f 2 /(3 + ) (b). The two choices of f 1 are very close, but whereas the former choice of f 1 is commensurate with f 2 , the latter choice is not. The negative values ofP (k) within the nonadiabatic region (which we attribute to nonadiabatic heating), are much fainter in panel (a) than in panel (b). This is consistent with our expectation that nonadiabatic heating is indeed significantly suppressed for f 1 = 2 3 f 2 compared to 2 3+ f 2 . We also compute the total dissipated power in the system due to both driving modes,P dis (k) = − lim t→∞ t 0 ds t Tr[∇Ĥ(k, s)ρ(k, s)] · (E 1 (s) + E 2 (s)), for the same parameters as in panels (a) and (b).P dis (k) measures the time-averaged rate of work done on electrons with wavevector k by the two driving modes in combination; hence it measures the total rate of dissipation, and is guaranteed to be positive due to the second law of thermodynamics. In Fig. 5(c) and (d) we plot P dis (k) for the parameter sets we considered in panels (a) and (b), respectively. While outside the nonadiabatic region,P (k) andP dis (k) effectively take the same values for the two frequency ratios, nonadiabatic heating is much weaker in panel (d) than in panel (c). The very different values ofη at frequencies f 1 = 2 3 f 2 and f 1 = 2 3+ f 2 must therefore be due to this suppression of nonadiabatic heating in the commensurate case.

D. Nonanalytic amplitude dependence
We next explore the relationship between the topological frequency conversion rate, the amplitudes of the incoming modes. Fixing E 2 = 2E 1 , Fig. 6(a) plotsη as a function of E 1 for the isolated Weyl node studied in the previous subsections, with f 1 = 3 2 f 2 = 1.23 THz, τ = 516.3 ps, and µ = 115 meV. The error bars in Fig. 6(a) indicate the estimated uncertainty due to the finite number of k-points we sample [101]. The conversion rate exhibits a clear cusp when E 1 ≈ 1000 kV/m; by inspecting the k-dependent frequency conversion rate as in Fig. 4, we verified that this is the amplitude where topological frequency conversion sets in due to the surface B 0 crossing the Fermi surface. The cusp ofη reveals a non-monotonous and nonlinear dependence on driving amplitude, supporting our conclusion that topological frequency conversion is an effectively nonperturbative response phenomenon. The amplification threshold is reached at amplitude E 1 ≈ 1300 kV/m.

V. CONDITIONS FOR FREQUENCY CONVERSION
There are several conditions that a Weyl material must satisfy to realize topological frequency conversion. The conditions can be grouped into the conditions that

A. Conditions on individual Weyl nodes
The rate of topological frequency conversion from Weyl node i is given by the ith term in the sum over Weyl nodes in Eq. (29): Therefore only the Weyl nodes sufficiently near the Fermi energy can contribute to topological frequency conversion. If a Weyl node is too far from the Fermi energy, the two touching bands are either both full or empty within a distance ∼ eA/ from the Weyl node, implyinḡ For the most natural case where each of the two modes contains only a single harmonic, W (k) = −W (−k), as is also evident in Figs. 2(c) and 4(a). For Weyl point i to contribute to frequency conversion,ρ 1 (k) orρ 2 (k) must hence break inversion symmetry around k i . Specificallȳ This constitutes our second condition. This symmetry breaking can be achieved with a nonzero value of the Weyl cone tilt, V [see Eq. (1)], as we demonstrated in our numerical simulations (Sec. IV).

B. Condition on symmetry class
We now identify the symmetries a Weyl semimetal must break to support topological frequency conversion.
The two symmetries that are central to the Weyl semimetals are the inversion and the time-reversal symmetry [12]; at least one of these symmetries must be broken for the Weyl nodes to exist. Both inversion and timereversal symmetry results in inversion-symmetric energy bands, ε α (k) = ε α (−k) (with α denoting the band index after indexing them according to their energy). Thus, for both symmetries, a Weyl point at wavevector k implies the existence of a Weyl point at wavevector −k. The conjugate Weyl nodes at k and −k have equal charges for time-reversal symmetric Weyl semimetals, and opposite charges for inversion-symmetric Weyl semimetals [12]. We expect the steady-state to approximately inherit the same inversion symmetry, such thatρ α (k) =ρ α (−k). For the most natural case where modes 1 and 2 each contain a single harmonic, W (k) = −W (−k). Hence the contributions toη 0 from symmetry-conjugate nodes cancel out for Weyl semimetals with time-reversal symmetry, but not for Weyl semimetals with inversion symmetry.
We conclude that broken time-reversal symmetry is required for topological frequency conversion, while inversion symmetry does not need to be broken. Other crystal symmetries, such as reflection and discrete rotation symmetry, do also not preclude frequency conversion, since the incoming modes (and hence W (k)) can be configured in a way that breaks these symmetries. Hence magnetic Weyl semimetals, such as Co 3 Sn 2 S or Co 2 MnGa [77,79], intrinsically support topological frequency conversion, while non-magnetic Weyl semimetals (such as TaAs) require an externally-provided timereversal symmetry breaking. This external symmetry breaking is already achieved with the circularly-polarized driving; a higher degree of asymmetry can further be achieved with e.g. a current bias or externally applied magnetic field.

C. Condition on driving
Next, we identify the conditions that the driving amplitudes and frequencies must satisfy to support frequency conversion for a given Weyl semimetal.
Sec. III C 2 concluded that the dynamics of electrons is non-adiabatic within a distance d 0 from the boundary B 0 , where d 0 ∼ eE v0 , and v 0 is the smallest singular value of the matrix R in Eq. (1) [see Eq. (36)]. For a nonzero number electrons to act as frequency converters, d 0 must hence be smaller than the linear dimension of B 0 , which we estimate to be of order eA/ . These considerations imply that is required for frequency conversion. Hence, "steep" Weyl points (i.e., Weyl points with large v 0 ) are most useful for frequency conversion, as they support topological frequency conversion at lower intensities.
Weyl points in known compounds support topological frequency conversion at experimentally accessible parameters: for example, TaAs has Weyl points for which v 0 ∼ 10 5 m/s [98,99] At frequency Ω ∼ 2π THz, we hence expect these Weyl points can support topological frequency conversion at moderate intensities of order 100 W/mm 2 and above.
Finally, we require that the bandwidth of the bands containing the Weyl node be larger than the driving frequency; otherwise, driving cannot be considered adiabatic anywhere in the system. This puts an upper limit for the frequencies that could achieve frequency conversion in a given Weyl semimetal. As an example, for TaAs the characteristic band gap between Weyl points is of order 20 meV [98], corresponding to a maximum frequency limit of ∼ 5 THz.

D. Condition on relaxation
A final condition for amplification, is that the rate of topological energy conversion,η 0 , must overcome the (negative) rate of dissipation,η dis . Our analysis and numerical simulations identified two sources of dissipation: momentum relaxation (η mr ) and nonadiabatic heating (η na ): η dis = η mr + η na . Amplification of mode 1 thus requiresη In Sec. III C 2 we concluded that the gross rate of topological frequency conversion, η gross (i.e., the rate that results when not taking into account cancellation between electrons that convert energy in opposite directions), can only exceed η na if τ 1/Ω. Since the net rate of topological frequency conversion,η 0 , is just a small fraction of η gross , and since energy is also lost to momentum relaxation, we hence expect net amplification (η > 0) can only be achieved when The phenomenological discussion in Sec. III C 2 shows that η na ∝ (1 − e −Γ∆t ) where ∆t denotes the characteristic time between instances where a given wavevector k is taken to the Weyl point by the applied drive (k → k + eA(t)/ ). Here ∆t is significantly smaller for Lissajous conversion (highly rational frequency ratios) than for incommensurate frequency ratios. Thus, the threshold for net amplification is significantly lower at commensurate frequencies. Indeed, in our simulations, a small adjustment of f 1 from 2f 2 /(3 + ) to 2f 2 /3 lowered the amplification threshold from above 1200 ps to ∼ 300 ps.
Our quoted values in Sec. III provide an example of how to estimate the break-even relaxation rate. For intensity E ∼ 1600 kV/m and Ω ∼ 2πTHz, we estimated η gross ∼ 500 kW/mm 3 . The net rateη 0 will be only a fraction of this value. For the same parameters, and with isotropic band gap matrix R of order 5·10 5 m/s and Fermi surface area 0.06Å −2 , we found η mr ∼ 5 · 10 −8 kJ/mm 3 τ , and η na ∼ 5 · 10 −8 kJ/mm 3 τ (for incommensurate frequencies) or 3.5·10 −9 kJ/mm 3 τ (for Lissajous conversion at frequency ratio 2/3). Hence, we expect topological frequency conversion can exceed the rate of dissipation when τ are several times larger than 100ps. We moreover expect the threshold to be significantly lower for Lissajous conversion (i.e. commensurate frequencies), than for incommensurate frequencies. This is in good agreement with our data in Fig. 1(b) which indicate that τ must exceed 300 ps in order to achieve net frequency conversion in the Lissajous regime for the parameters above, and 1000ps for incommensurate frequencies.
The different scaling behaviors of dissipation and topological frequency conversion point toward the parameter regimes beneficial for amplification.
First, note that (at a fixed area of the Fermi surface), η mr scales linearly with electric field E, while η na and η gross scales as E 3 (specifically, η na ∼ I 3 in the Lissajous regime) [see Eqs. (34), (37), and (30)]. Thus, we expect that the relative contribution of η mr decreases at high intensity, while the ratio of η na and η gross remains fixed, implying that frequency conversion becomes more efficient at higher intensities.
Second, for a given intensity, the topological frequency conversion rate scales as 1/Ω, while η na scales as 1/Ω 3 (for incommensurate frequencies) or 1/Ω 2 (for commensurate frequencies). Similarly, momentum relaxation scales as 1/Ω 2 . Thus, we expect amplification is most easily reached at the top of the frequency range that supports topological frequency conversion, given the driving intensity and band structure of the system.
The requirements on the relaxation rate pose the biggest current challenge to realizing topological frequency conversion. Relaxation times in known Weyl semimetals have been reported to be in the range 0.25 − 3ps [62][63][64]67], although transient signatures with lifetimes above 100 [67,68] and 1000 [69] ps have also been reported in some compounds. Thus further improvements in the quality of materials are needed to fulfill the requirements of topological frequency conversation in the practically interesting THz range.

VI. THINKING OUTSIDE THE GRAIN: GLOBAL ELECTRODYNAMICS CONSIDERATIONS AND IMPLEMENTATION USING PHASE ARRAYS
The full understanding of the frequency conversion effect requires thinking about the global electromagnetic field, and the material response of the Weyl grains to an external drive. Specifically, in this section we incorporate the dielectric response to our analysis, and propose a phase-array geometry of the Weyl grains as a prototype for a Weyl topological amplifier.

A. Renormalization of electric field by plasma oscillations
Let us begin with considering the macroscopic response of a single grain to external driving. For i = 1, 2, we let E i (r, t) denote the (plane-wave) electric field from mode i as a function of position r and time t and let E 0 (r, t) = E 1 (r, t) + E 2 (r, t) denote the net "incoming" field resulting from the driving. The current and charge oscillations in the grain induced by the external driving creates an additional electric field, E ind (r, t). The total electric field inside the sample is thus given by E(r, t) ≡ E 0 (r, t) + E ind (r, t); this is the field driving the response of the material, and is the one we considered in the calculation in the previous sections. Evidently, the internal field in the sample gets renormalized by the charge and current in the material.
Our first order of business is to find the internal field E(r, t) (which we used in our analysis above) in terms of the external fields. We can find E(r, t) self-consistently by solving Maxwell's equations, taking account the current and charge dynamics in the grain induced by E(r, t).
While an exact (geometry-dependent) analysis is in principle possible, the small size of the grain allows us to make some simplifications, such as ignoring the skin effect. Thus, inside the grain E ind (r, t) is approximately given by the electrostatic field resulting from the instantaneous charge configuration in the system. These in turn arises from the driving-induced oscillations of the grain's bulk plasmonic mode [102].
For a small spherical grain, we can assume E 0 (r, t) uniform within the grain, and moreover ignore retardation effects of the electromagnetic field (this is equivalent to neglecting the skin effect). Inside the grain, E ind (r, t) is thus given by the electrostatic field arising from the instantaneous charge distribution at time t. The charge distribution is nontrivial due to the oscillating currents, which produce surface charges. Specifically, ∂ t ρ(r, t) = ∇ · j(r, t), implying ρ(r, ω) = i ω ∇ · j(r, ω). We now show that the equations of motion above have a solution in which the current density and E ind (r, t) are also uniform within the sample. To show this, note that a uniform current density in the grain, j(r, t) = j(t), implies that the charge accumulates on the surface. The surface charge density at the angle specified by unit vectorr on the sphere, is given by ∂ t λ(r, t) = j(t) ·r. Hence λ(r, t) =r · λ 0 (t), with λ 0 (t) denoting the unique zero-mean solution to ∂ t λ 0 (t) = j(t). Inside the sphere, the electrostatic field from a surface charge distribution λ(r) =r · λ 0 is uniform and given by E ind (t) = λ0(t) 3 0 . Hence, the electric field is uniform within the sample and given by E(t) = E 0 (t) + λ0(t) 3 0 . Thus, a uniform current density j(t) and E(t) solves the dynamics of the grain.
Transforming to frequency domain, and using that j(ω) = −iωλ 0 (ω), we finally arrive at This gives the frequency-dependent renormalization of the electric field inside the grain, and is an exact solution in the limit where the grain size is smaller than the wavelength of the driving modes.
In linear response theory, the time-derivative of the current response is assumed proportional to the electric field, implying j(ω) ≈ −iκE(ω)/ω for some constant κ. The resulting solution leads to a frequency dependent relative permittivity, E(ω) ≈ (ω)E 0 (ω) with (ω) = (1 − ω 2 p /ω 2 ) −1 with ω p = σ/3 0 denonting the plasma frequency of the system. The plasma frequency omega p is estimated for generic Weyl semimetals, in Ref. [103]: it is typically given by an O(1) constant times the Fermi energy.
The linear response analysis above is useful for elucidating the qualitative features of the plasmonic response. However, the regime we consider potentially supports a significant nonlinear response due to the nonquadratic dispersion and large Berry curvature surrounding Weyl nodes -indeed, topological frequency conversion is a nonlinear response phenomenon. We thus go beyond the linear response regime in our analysis below: for a given internal field configuration, E(t), the current response j(t) can be easily computed in the limit of weak relaxation and adiabatic driving without any linear response approximation, using Eq. (21), j(t) = j 0 (t) + δj(t) with j 0 (t) given in Eq. (23) and the dissipative component δj(t) is negligble in the limit we consider [104].
The driving frequency controls whether the plasma oscillations amplify or screen the electric field from the incoming radiation. This qualitative behavior is evident in the linear response result we quoted above, but also endures after taking into account the nonlinear response. To see this, consider what external electric field E 0 (ω) is needed to cause a given internal field E(ω) [which determines the current response j(ω)]. As in the linear response regime, j(ω) is controlled by vector potential, and thus scales with E(ω)/ω [see Eq. (23)]. The plasmoninduced electric field hence is negligible in the limit of large ω (but grain size still smaller than the wavelength), meaning the grain is effectively transparent to the radiation: E(ω) ≈ E 0 (ω). Conversely, for small ω, j(ω)/ω 0 will be considerably larger than E(ω), implying that E 0 in turn has to be much larger in E(ω) for Eq. (55) to hold. Thus, for small frequencies, the plasma oscillations severely screens the electric field inside the sample relative to the external field. At some intermediate fre- 3 0ω , and a very weak external field thus causes a large internal field. In this case, driving resonates with the plasma oscillations, causing significantly enhanced amplitude of the electric field. As we will see, this mechanism allows for significant enhancement of the topological frequency conversion rate.

B. Radiation output of grain
Next, we want to compute E(r, t) outside the grain, to determine the profile of the emitted radiation. Here, the grain's small size means that E ind (r, t) to a good approximation takes the form of dipole radiation generated by some nontrivial trajectory of the dipole moment. Using that the surface charge distribution we obtained above, we find the dipole moment to simply be given by p(ω) = 4π 3 ir 3 j(ω) ω . The energy converted to mode 1 leaves the grain as radiation energy at frequency ω 1 . The bulk of the emitted radiation energy comes from constructive interference between the incoming plane wave E 0 (r, ω 1 ) and the emitted dipole radiation E ind (r, ω 1 ) (i.e., the ω 1 -Fourier components of E ind (r, t) and E 0 (r, t), respectively).
To compute the frequency-resolved radiation energy emanating from the grain, we consider the total energy flux density, given by the Poynting vector field, S(r, t) = 1 µ0 E(r, t) × B(r, t). By using the Fourier decomposition E(r, t) = dωE(r, ω)e −iωt along with E(r, ω) = E * (r, −ω) [and likewise for B(r, t)], we find that the time-averaged energy flux density,S(r), is thus given byS We identify S(r, ω) ≡ 1 µ0 Re [E(r, ω) × B * (r, ω)] as the energy flux density from modes with frequency ω. The total radiation power from the grain at frequency ω is given by where C is some surface containing the grain. To compute P (ω), we use the divergence theorem to find P (ω) = 2 µ0 C dV Re [∇ · (E(ω) × B * (ω)], with C dV denoting the integral over the volume contained in C. Next, we apply the cross product identity ∇ · (E × B * ) = −E · (∇ × B * ). Using Ampere's law ∇ × B(r, ω) = −iµ 0 0 ωE(r, ω) − µ 0 j(r, ω), where j(r, ω) is the Fourier transform of the current density, yields P (ω) = 2 µ0 dV Re [iω 0 |E 2 (ω)| + E(ω) · j * (ω)]. The first term in the parenthesis evidently is fully imaginary, and thus gives a vanishing contribution to the integral. This leave us with P (ω) = 2 dV Re [E(r, ω) · j * (r, ω)]. Since E(r, ω) and j(r, ω) are uniform within the grain, we find which is exactly the quantity we calculated in Sec. III. Interestingly, the plasma-induced electric field does not directly contribute to the power output, since it is proportional to −ij(ω); rather it indirectly modifies the power output through its effect on the current response. We thus arrive at withη denoting the frequency conversion rate within the grain. This gives the output intensity of the dipole radiation emitted with frequency ω 1 .
C. Implementation using a phase array In order to produce an amplifier out of the frequencyconversion effect, we must consider combining many grains together to create a phase array. While we do not intend to analyze such a device in detail in this manuscript, we will here outline its design. The phasearray geometry we envision has Weyl grains arranged in a 3d cubic-lattice, with one axis along the propagation direction of an external plane wave with circularly polarization (mode 1), which we intend to amplify. In the direction of the propogation of mode 1 the structure will have a lattice constant of a quarter wavelength, such that the backscattering element of the amplified mode 1 is be eliminated by destructive interference.
The array will also be subject to a normal-incident radiation from mode 2, which is the amplification source beam to be converted. In order to maximize the amplification effect, we anticipate that in-phase arrangement of all layers with respect to mode 2 would be beneficial, hence the lattice constant along the mode 2 direction should be mode-2 wavelength,λ 2 (see Fig. 7).
Indeed, obtaining amplifiers from single gain elements is a common practice. Josephson traveling wave amplifiers (JTWA; see, e.g., Ref. 105) , for instance, are essentially a chain of individual Josephson parametric amplifiers. Other traveling wave optical amplifiers rely on nonlinear crystals such as LiNb or β−BaB 2 O 4 (BBO). It is only in the macroscopic constructs of parametric JTWA, and non-linear crystals which mix light modes, that issues related to phase-matching arise. Similarly, a single Weyl grain in a topological-frequency-conversion regime requires no mode-matching beyond the need to have a rational frequency ratio to be in the Lissajous regime. In fact, the grains are expected to be smaller than the wavelengths involved. It is only when these grains are combined into an array that we need to consider how the wavelengths of the amplified radiation correspond to the spatial structure of the amplifier.

D. Numerical simulations of plasmon-enhanced amplification
We now demonstrate that topological frequency conversion remains possible even after including screening effect of plasma oscillations. Furthermore, we show that tuning frequencies near the plasma resonance dramatically enhances the frequency conversion effect.
We consider an inversion-symmetric Weyl semimetal whose Fermi surface consists of two Weyl nodes that are related by inversion symmetry. The Hamiltonian of one Weyl node is given by H(k) = vk · σ + k · V with v = 3.87 · 10 5 m/s and V = (0, 0, 3.1 · 10 5 m/s) (i.e., the same dispersion as considered in Sec. IV); the Hamiltonian of the other Weyl node is given by H(−k). We consider the case where the electric field inside the grain, E(t), is fixed and given by two circularly polarized FIG. 7. Phase-array proposal for combining the gain from Weyl grains into a topological amplifier. The spacing of the grains in the direction of propagation should match a quarter of the wavelength of the amplified wave. This will suppress reflection, and will concentrate the contribution of each grain into a single forward-propagating beam. The pump beam is expected to be normal to the source beam, and the spacing along its direction should be its wavelength.
In Fig. 6(b) we plot the resulting relative increase of mode 1 inside the material, G = |E(ω 1 ) · E * 0 (ω 1 )|/|E 0 (ω 1 )| 2 , as a function of chemical potential, µ (blue). We also plot the corresponding relative increase of mode 2 (orange), as well as the net relative gain of all remaining modes (induced by the nonlinear oscillations of the plasmons), G 3 = ( dω|δE 0 (ω)| 2 ) 1/2 , with δE 0 (ω) denoting the component of external field E 0 (ω) which is orthogonal to the internal field E(ω) (for ω / ∈ {ω 1 , ω 2 }, δE(ω) = E 0 (ω)). The data in Fig. 6(b) shows that the plasma oscillations do not affect the modes when the chemical potential is smaller than 2.5 meV. Moreover, for µ 5 meV, G 3 1, implying that the external radiation field does not need lead to appreciable amplitude of any higher-harmonic or orthogonal modes to provide the bichromatic electric field inside the grain which we require. In other words, the plasma oscillations do not significantly excite any modes other than the pump and signal mode when µ 5 meV. For values of µ above 5 meV, the internal electric field gets severely suppressed by the plasmon screening, while the plasma oscillations begin to significantly excite modes other than the pump and signal modes. Here, the frequency conversion rate is significantly reduced. Furthermore, a more careful analysis is needed in this regime to account for the higher harmonics of E(t) induced by the nonlinear plasma oscillations.
In the range 2 meV µ 4 meV, the internal field is significantly enhanced by the plasma oscillations, without nonlinear corrections playing a role. This plasma resonance dramatically enhances topological frequency conversion: we first compute the frequency conversion rate for the same parameters considered for Fig. 6(b), using the approach of Sec. IV. From our obtained frequency conversion rate,η, we compute the gain coefficient of the material, g =η/I 1 with I 1 = cε 0 |E 0 (ω 1 )| 2 denoting the intensity of mode 1 outside the material; we use the data from Fig. 6(b) to compute E 0 (ω 1 ) (recall we consider a fixed value of the internal field, E(ω 1 ), but allow E 0 to vary). The gain coefficient has dimension of inverse length, and gives the characteristic rate at which mode 1 gets amplified inside the material. In Fig. 1(b), we plot the gain coefficient as a function of chemical potential, using τ = 200 ps (blue), 400 ps (orange) and 600 ps (green). When the plasmon resonance is reached at µ ≈ 4 meV, the gain coefficient increases dramatically, reaching values of order 100 cm −1 , exceeding, e.g., the THz gain coefficients of 20-50 cm −1 reported in Refs. [57,58].

VII. DISCUSSION
In this manuscript, we showed that Weyl semimetals can efficiently convert energy between two driving modes, through the mechanism of topological frequency conversion [52]. This effect makes Weyl semimetals promising media for THz and possibly even infrared amplification. Our analysis shows that Weyl semimetal with feasible band dispersions support topological frequency conversion in the "THz gap" at experimentally accessible intensities of order ∼ 50 W/mm 2 , or even less (∼ 1 W/mm 2 ) if one drives near the plasma frequency. Topological frequency conversion is supported both for incommensurate frequencies and commensurate frequencies, but is most efficient in the latter case, due to the mechanism of Lissajous conversion. Our numerics and estimates focused on topological frequency conversion in the THz regime, where there is the biggest need for new photonic control elements, but the effect may also be supported at other frequency ranges.
The primary obstacle to Weyl semimetals operating as topological frequency converters is drive-induced heating. Heating both wastes energy from the beams we would like to amplify, and may even damage the material. Phonons, interactions and impurities all lead to electron relaxation processes which cause this heating.
Through phenomenological arguments and numerical simulations, we identified two important mechanisms for dissipation: momentum relaxation and nonadiabatic heating. Momentum relaxation occurs when electrons near the Fermi surface relax their energy by changing their momentum, and which is common to all irradiated materials. Nonadiabatic heating emerges when electrons undergo Landau-Zener transitions between the valence and the conduction band. This mechanism is particularly relevant in topological semimetals, due to the existence of gap-closing points in these materials. Even so, nonadiabatic heating is strongly suppressed in the Lissajous regime, which makes it much preferred for amplification.
In our simulations and phenomenological discussion, relaxation was parameterized through a single relaxation time, τ . In particular, we took electron-hole recombination, and intra-band momentum relaxation (which is supported by phonons) to have the same characteristic rates. Needless to say this treatment could be made more realistic by considering separate relaxation rates for these processes, as suggested by experiments [64,67,68]. Nonetheless, we believe our simple dynamical model captures the conditions for amplification.
To achieve amplification, where energy gain due to topological frequency conversion exceeds the loss due to dissipation, the characteristic relaxation time τ must be sufficiently long. To limit non-adiabatic heating as well as momentum relaxation we need τ f 1. This condition was clearly evident in our simulations: even for optimal parameters, and in the Lissajous regime, breakeven was only reached when 1/τ 300f ,(for incommensurate frequencies, amplification required τ f > 1000). So far, τ 's were reported in the range 0.25 − 5 ps [62][63][64] in pump-probe experiments. This suggests net amplification of continuous-wave THz frequencies is currently beyond reach. That said, transient experimental signatures with τ > 100 ps have been seen in Weyl semimetals [67,68], emphasizing that a more discriminating analysis may reveal a broader amplification regime.
Notwithstanding, signatures of topological frequency conversion effect could be observed even if the relaxation time is too short to allow for amplification. That is because the direction of energy conversion of a Weyl semimetal driven bichromatically by two circularlypolarized modes is determined by the product of the two mode's polarizations. Hence a reversal of the circular polarization of either mode should will lead to an increase in the output intensity of one mode, and a decrease for the other mode. This topological effect could be accessed experimentally. Another group recently proposed to utilize this chirality-sensitive intensity shift to extract enantioselective information from a gas of chiral molecules [106].
Furthermore, we can suggest several strategies to approach the amplification regime. Commensurate frequency conversion, i.e. Lissajous conversion, already provides a dramatic improvement by suppressing nonadiabatic effects by an order of magnitude. Momentum relaxation is harder to control. Note, however, that momentum relaxation energy loss scales linearly with radiation intensity, I. In contrast, topological energy conversion (and nonadiabatic heating) scale as I 3/2 . Therefore the relative significance of momentum relaxation should decrease at larger intensities. Moreover, at a given intensity, topological frequency conversion scales inversely with the driving frequency, f , while dissipative energy absorption decreases as f −2 (specifically, η na ∼ f −2 in the Lissajous conversion regime). The amplification threshold of τ f will therefore be lower at higher frequencies.
If an issue, driving-induced heating can possibly be circumvented by using pulsed lasers instead of continuous wave beams: by allowing the system to dissipate heat between pulses, such a scheme would allow us to reach the high-intensity regime without causing material damage; while a detailed investigation would be an interesting topic for future studies, we expect pulses with durations more than a few periods, or randomly-timed pulses, will yield conversion rates consistent with topological frequency conversion at continuous-wave operation. This way the large-amplitude regime required for frequency conversion could be realized while allowing time for the system to dissipate absorbed heat between pulses even if relaxation times are short. In addition to these considerations, materials with a steeper velocity makes realizing the large frequency regime easier, as the velocity at the Weyl point is the "coupling constant" that converts the electric-field amplitude into an energy scale.
Weyl nodes need to be located near the Fermi surface to support topological frequency conversion, and moreover need to be surrounded by an asymmetric electron distribution, in order to ensure an imbalance in the numbers of electrons that convert energy at opposite rates. Optimal imbalance can be reached in the presence of a "Weyl cone tilt", and through appropriate tuning of the chemical potential. Additionally, our analysis indi-cates that time-reversal symmetry needs to be broken to acheive frequency conversion. Hence, we expect magnetic Weyl semimetals, such as Co 3 Sn 2 S or Co 2 MnGa [77,79] are best-suited for topological frequency conversion.
Topological frequency conversion could also be achieved in non-magnetic Weyl semimetals, or even be enhanced in magnetic Weyl semimetals, by "priming" the particle distribution into an out-of-equilibrium state. Such priming could e.g. be achived by driving the system with ultrashort laser pulses or with a DC current, and would create a transient state more suited for frequency conversion than the steady states we have considered in this work. Similarly, purification of the material, alongside bath or substrate engineering are other potentially important directions for realizing amplification by suppressing dissipation. Indeed, these research directions are also important for the general nonlinear response of Weyl semimetals (e.g., chiral photogalvanic effect) [20].
Acknowledgements -We thank N. Peter Armitage, Chris Ciccarino, Cyprian Lewandowski, Prineha Narang, and Mark Rudner for valuable discussions. where A αβ ≡ i ψ α |∇ψ β . Combining Eqs. (B67)-(B69), we hence find We identify the first term in the right-hand side above as the contribution arising from the group velocity.
To rewrite the second term, we use that |∂ t ψ α = e E · |∇ψ α , implying M αβ = e E · A αβ [114]. Thus αβ and E i denotes the ith vector component of A αβ and E, respectively. Next, we note (B73) We identify the right-hand side as −i α,k ρ α ijk Ω k α , where ijk denotes the Levi-Civita tensor and Ω k α denotes the kth component of Ω α . Thus, Hence the second term in Eq. (B70) gives the contribution to the particle velocity from the anomalous velocity. In particular, by inserting the above result into Eq. (B70), and dividing through with , we establish Eq. (B66), which was the goal of this subsection.

Derivation of auxiliary results
In this subsection we derive the auxiliary results which we quoted in the subsections above. Specifically, we derive Eqs. We first show that ρ α is nearly stationary. Our starting point is the equation of motion for the diagonal matrix elements ofρ(k, t) in the orbital basis, {f n (k, t)}, ∂ t f n (t) = − k R mm nn (t)f m (t). We note that R mm nn (k, t) is of order Γ, but oscillates with characteristic frequency Ω. As a result, we expect f n (k, t) to deviate from its time-average,f n (k), by a correction of order Γ/Ω: Inserting this into Eq. (B58), we thus find whereρ α (k) denotes the time-average of ρ α (k, t), and we neglected a correction of order λ 2 , since it is subleading relative to Γ/Ω.