Are Microwave Induced Zero Resistance States Necessarily Static?

We study the effect of inhomogeneities in Hall conductivity on the nature of the Zero Resistance States seen in the microwave irradiated two-dimensional electron systems in weak perpendicular magnetic fields, and we show that time-dependent domain patterns may emerge in some situations. For an annular Corbino geometry, with an equilibrium charge density that varies linearly with radius, we find a time-periodic non-equilibrium solution, which might be detected by a charge sensor, such as an SET. For a model on a torus, in addition to static domain patterns seen at high and low values of the equilibrium charge inhomogeneity, we find that, in the intermediate regime, a variety of nonstationary states can also exist. We catalog the possibilities we have seen in our simulations. Within a particular phenomenological model, we show that linearizing the nonlinear charge continuity equation about a particularly simple domain wall configuration and analyzing the eigenmodes allows us to estimate the periods of the solutions to the full nonlinear equation.


I. INTRODUCTION
A novel zero-resistance state (ZRS) was observed a few years ago, when a two-dimensional electron gas (2DEG), subjected to a weak magnetic eld, was also irradiated with microwaves. 1,2,3,4In the ZRS state, the dc resistance measured in a Hall bar experiment appeared to vaniish, while correspondingly, the longitudinal conductivity σ xx , measured in a Corbino geometry appeared to be zero.A related phenomenon of microwave-induced resistance oscillations had been observed earlier in samples of lower quality, in which there are changes in the dc resistance produced by the microwave irradiation, which depend in an oscillatory manner on the ratio of the microwave frequency to the cyclotron frequency of the 2DEG, but where the resistance was not driven down to zero.

5,6
Following the discovery of the ZRS, a phenomenological explanation was put forward. 7It was assumed that for microwave power above a certain threshold, for an appropriate sample, in certain ranges of the magnetic eld and microwave frequency, the microwave-induced resistance oscillations would increase to the point where the dierential longitudinal conductivity would become negative at zero dc electric eld, and it was noted that in this case that the zero-electric-eld solution would necessarily be unstable.It was argued that the sample would spontaneously break up into a pattern of domains, separated by domain walls.Within each domain, the magnitude of the electric eld would be a constant, E c , dependent on the strength of the microwave eld and other parameters, which satises a condition that the longitudinal current generated by the eld is zero.The direction of the electric eld would change discontinuously at a domain wall, and might also vary continuously within a domain.It was assumed that the domain walls could move rather easily, so that, e.g., if a nite voltage is applied between the inner and outer edges of a Corbino sample, the system would respond by a displacement of domain walls, without causing the magnitude of the eld to deviate from the critical value E c , and without inducing a longitudinal current in the sample.
If one accepts this general picture, an obvious question is what determines the domain pattern in any given sample?In fact, one may question whether there is necessarily a static domain pattern at all.Since the system is driven out of equilibrium by the constant absorption of microwave radiation, it is possible in principle that the system will enter a time-dependent state, where domain walls may move about in a chaotic, or possibly quasiperiodic, manner.It is our purpose here to address this question.We shall argue that time-dependent states are likely to occur in at least some situations, and we shall propose one experimental geometry where it should be possible to observe this time-dependence.The dc properties of a time-dependent state should be similar to those of a state with static domain walls; in particular the observed value of σ xx should be zero, or at least very small compared to the conductivity in the absence of microwave radiation.
At least two dierent microscopic mechanisms have been proposed to explain the microwave-induced resistance oscillations at lower microwave intensities, and to produce the negative zero-eld conductance assumed in the macroscopic model of ZRS: a displacement mechanism 8,9,10,11,12,13,14,15,16 , in which the absorption of a microwave photon by an electron leads (for appropriate values of parameters) to a disorder-assisted displacement of the cyclotron orbit in an up-hill direction with regard to the local dc electric eld; and a population mechanism 17,18 , in which absorption of the microwave photons leads to a population inversion in the partially lled Landau levels close to the Fermi level.Although these mechanisms may lead to quite dierent estimates for such parameters as the threshold microwave power or for the critical eld E c at a given level of microwave irradiation, it appears that they give rise to qualitatively similar phenomenological models.In any case, we shall not make any assumptions about the particular microscopic mechanism in this work.
In formulating the equations of phenomenological model, it is convenient to divide the local electric current into a longitudinal or dissipative part, and a Hall current which is always perpendicular to the local eld E.
arXiv:0806.1562v1 [cond-mat.mes-hall]10 Jun 2008 In previous work, 19,20 we considered a model where the functional form of the dissipative current could vary from one place to another in the sample, due to small inhomogeneities in the doping density, or due to other sources of disorder.However, we assumed the Hall current to be governed by a linear (in eld E) Hall conductivity, which we took to be uniform throughout the sample.In this case, we were able to derive a Lyapunov functional, whose value can only decrease in time, for a system with specied electrochemical potential at the boundaries.Then the system must eventually reach a static state, whose potential conguration is at least a local minimum of the Lyapunov functional, so that oscillatory or chaotic timedependent solutions are not possible as a steady state.
If we allow the Hall conductivity to vary from one place to another, however, it is generally not possible to nd a Lyapunov functional for the system.In this case, there is no minimization principle to determine the long-time behavior, and there is no guarantee that a stable timeindependent solution exists.In fact, we shall explore explicitly some simple geometries where, for certain ranges of parameters, the long-time behavior is time-dependent.
The rest of the paper is organized as follows.In Sec.

II. MODEL
In the presence of an external microwave eld, we assume the following relation between the dc current j(r) and the eld E(r): where E ≡ −∇V (r) ≡ ∇µ(r) e and ẑ is a unit vector normal to the plane.Here, µ(r) is the electrochemical potential and −e is the electron charge.For the sake of convenience, we will from here on refer to E as the electric eld and V (r) as the electrochemical potential.In Eq. 1, we will explicitly allow the Hall conductivity σ H to be position-dependent.The dissipative current j d satises the condition that the dierential dissipative conductivity, σ d αβ (E) ≡ ∂j d α /∂E β , is symmetric: The nonlocal (third) term in Eq. ( 1) implements an ultraviolet cuto, which will lead to a nite domain wall thickness, proportional to the square-root of the parameter λ, which will be taken to be small compared to the size of typical domain (indeed, this is almost a matter of denition of domains phase).In practice, we expect that the domain wall thickness will be of the order of the cyclotron radius l c , which is of order 1 µm in typical samples where ZRS is observed.
The vector function j d may also depend explicitly on the position r, due to inhomogeneities in the 2DEG, and its direction may not be perfectly aligned with E. In previous work, we considered explicitly the eects of inhomogeneities in j d arising from gradients in the equilibrium electrostatic potential φ d due to disorder.In the present paper, however, we shall ignore this complication, and shall assume that the function j d has no explicit dependence on r, except for a brief discussion at the end.
Eq. (1) will be supplemented by the continuity equation, where ρ is the charge density.Writing E ≡ −∇V (r), we may relate changes in the electrochemical potential V (r) to changes in ρ through the inverse capacitance matrix W : If a time-independent steady state is reached, then we have simply ∇ • j = 0, and the precise form of W is unimportant, but the form of W will be relevant for timedependent solutions.
In a Corbino geometry, one species the potential on the inner and outer boundaries of the sample, and one looks for a solution for V (r) consistent with these boundary conditions.If one assumes σ H to be a constant, the Hall current cannot contribute to ∇ • j in the interior of the sample, so it does not appear in Kircho 's equations.
Consequently, the solution V (r) is independent of σ H and we may, for simplicity set σ H = 0. To recover the Hall current, one simply inserts the resulting solution for E into the second term in (1), at the end of the calculation.Condition (2) on j d allows us to dene a scalar Lyapunov functional as g ≡ A variation of ( 5) is given by where j l = j d − λ∇ 2 E. The second integral vanishes on equipotential boundaries, or in the absence of external currents.Then, if σ H is independent of position, the extrema of G are found to be steady states, with ∇ • j = 0.
Using the positivity of the inverse capacitance matrix W , one may show that G[V (t)] is indeed a Lyapunov functional, i.e. a non-increasing function of time, so that its minima are stable steady states.In general, G may have multiple minima.Any initial choice of V (r) will relax to some local minimum of G, though not necessarily the ground state with lowest G. Nevertheless, one might expect that in the presence of noise, the system might tend to escape from high-lying minima and wind up in a state with G close to the absolute minimum.We shall assume here that the function g depends only on the norm squared E 2 of the electric eld, as would be appropriate in the limit of a uniform isotropic electron system, assuming that there is no unique axis picked out by external factors such as the polarization of the microwave eld.By hypothesis, under conditions where ZRS occurs, the function g must have its absolute minimum at a non-zero value of the electric eld, E = E c .
Expanding about this minimum, and keeping only the nontrivial terms of lowest order, we may write where the coecient σ c has the dimensions of a conductivity.In the absence of other information, we take σ c to be of the order of the dark conductivity (see Sec. VIII for estimates).When E = E c , the dissipative current, j d = ∂g/∂E, will vanish, and In the following calculations, for reasons of simplicity, and also to avoid introduction of additional parameters, we shall assume that Eq. ( 8) is exact for all values of E.
(Actually, only the range 0 < E < E c is important for our calculations.)We shall also assume that E c and σ c are independent of position.
The charge density ρ(r) is naturally broken into two parts, ρ(r) = ρ 0 (r) + δρ(r), where ρ 0 (r) is the density in thermal equilibrium with no incident microwave radiation (by denition, the electrochemical potential is constant everywhere under such conditions).In a Corbino setup, for example, ρ 0 (r) is the equilibrium density when the contacts at both edges are set to V = 0.In addition to having variations due to a non-uniform local dopant density, ρ 0 (r) could be tuned by a voltage on an external gate that is displaced from the 2DEG by a distance which varies from one point to another.The remaining term, δρ(r), is the non-equilibrium charge density produced by microwave-induced domain structure, and it is this contribution which is responsible for the position dependent electrochemical potential, according to Eq. ( 4).
In our calculations, we shall assume a local form for the capacitance matrix W , so that the electrochemical potential V (r) and the nonequilbrium charge density δρ(r) at a given point are simply proportional to each other: This form will be correct if there is a parallel conducting gate, set back from the 2DEG by a distance d small compared to the typical domain size, in which case we have C ≈ 0 /d, where is the dielectric constant of the material between the 2DEG and the gate.If a nearby conducting plane is absent, the local capacitance model will not be strictly correct, but we would expect to obtain qualitatively correct results by using (10), replacing the setback distance d by a characteristic domain size.We shall assume that the Hall conductance at point r is determined in the usual way by the charge density ρ at that point : Variations in the Hall conductance due to variations in ρ 0 (r) will be crucial for the eects we investigate below.On the other hand, variations in σ H due to variations in δρ(r) play no role in the dynamics when one assumes a local capacitance, as in Eq. (10).The contribution of this term to ∇σ H is parallel to ∇V , and, therefore, gives no contribution to the divergence of the Hall current.
Hence, we shall neglect this contribution and write the Hall conductance as which is a quantity xed at the outset, independent of the non-equilibrium charge-density induced by the microwave radiation.
We nally note that the Lyapunov functional of Eq. 5 misses one important piece of physicswhile it incorporates the fact that there is a high penalty for the eld to change on length scales shorter than the domain wall thickness λ/σ c , it neglects the fact that, at least in the local capacitance model, the domains cannot be arbitrarily big.More precisely, if we imagine drawing lines through the domain that are parallel to the eld inside of it, then these line segments inside the domain cannot get too long.Otherwise, the induced voltage dierence between points at opposite ends of these lines would result in an induced charge density δρ that alters the total electronic charge density by a large amount, which should be taken into account in the functional form of the Lyapunov function g(E).For example, one might need to take into account variations in the parameters E c and σ c in Eq. ( 8) due to variations in δρ.In a local capacitance model, the quantity where ρ0 is the mean electronic charge density in the 2DEG, will clearly provide an upper cuto to the size of any domain.Rather than modifying the functional, we shall assume here that the domains we consider are below such cuto.Note that the maximum domain size will tend to be largest if the setback distance d of the screening gate is large, so that the capacitance per unit area is small.

Dimensionless parameters
In our calculations below, we shall assume a rectangular sample with linear dimensions L x = L and L y = κL.
As postulated above, in the domains phase, the magnitude of electric eld away from the domain walls be equal to E c .We will then work with dimensionless units, with distances measured in units of L, electric eld E(r) = −∇V (r) measured in units of E c , electrochemical potential V (r) measured in units of E c L, conductivities measured in terms of the coecient σ c , and, nally, time measured in units of CL 2 /σ c .
The dimensionless Lyapunov density, c.f. Eq. ( 5), we ).When it is minimized, the electric eld is xed at 1, that is E c , away from the domain walls.Including the nonlocal contribution giving nonzero thickness to domain walls, the dimensionless longitudinal current density implied by the Lyapunov form we use is where l dw is related to parameter λ from Eq. 5 via l dw = 2 L λ σc .
In the reduced units introduced above, we will be solv- We have looked at two dierent sets of boundary conditions.In the rst one, corresponding to Corbino geometry, we assume periodic boundary along x, but the potential is being xed to be zero at the top (y = κ) and the bottom (y = 0) boundaries.In the second case, corresponding to a torus, the periodic boundary conditions along both directions are assumed.For simplicity, will only consider unidirectional Hall conductivity σ H (y).

III. UNIFORM HALL CONDUCTANCE
We rst consider the case where the equilibrium charge density, and, therefore, Hall conductivity σ H is uniform.
Then the domain patterns observed should follow directly from Lyapunov energetics.If γ is the angle that the electric eld forms with the domain wall, then, for the Lyapunov density that we have assumed, the cost of domain wall per unit length is proportional to sin 3 γ (See the Appendix for the proof ).If we neglect the energy cost associated with wall crossings and exponentially small interactions between well-separated walls, then we can describe the ground states, that is the lowest Lyapunov energy states, of the system with the two sets of boundary conditions introduced above.They are illustrated in the scope of this paper.
It should be emphasized that the existence of nonequivalent degenerate solutions in the case of κ = 0.5, and the multiple degenerate solutions for κ < 0.
Suppose that V 0 (r) Of course, this time dependence is not observable for a translationally invariant state, as in the case of a single horizontal wall in the center of the sample.However, if diagonal domain walls are present in the sample, then, if one were to detect electron density at a single point with, say, an SET, 21,22 then one would nd it changing with time.
The question of whether the states other than the one with a single horizontal wall would be present in an actual sample is hard to answer denitively.Given a random initial conguration for V (r), whether time evolution would take it to a state having diagonal walls is a matter of the relatives size of the domains of attraction for the such states compared with the domain of attraction for the state with one horizontal wall.Time-evolving the continuity equation above numerically (more details on the numerics are provided later in the paper), we found that, for κ ≤ 0.5, starting with a random initial guess for the electrostatic potential, it is very likely that we end up in the state with diagonal walls.Moreover, as mentioned above, if we were to assume a more realistic form for the Lyapunov function g than the quartic form (8), we believe that solutions with diagonal domain walls would most likely have lower Lyapunov energies than the solution with a single straight wall, in which case their basins of attractions would presumably be further enlarged.
We may also consider the situation where there is a voltage dierence V y between the top and bottom edges of the system.As long as V y is small compared to E c L y , the voltage can be accommodated by small displacements of the domain walls, giving a non-zero average of E y , without producing a net current in the y-direction.This can occur for the moving domains that we nd when σ H depends linearly on y in the same manner as for the stationary solutions appropriate to constant σ H . Thus, the property of zero dc conductance, as measured in the Corbino geometry, is preserved for the time-dependent solutions.

V. DYNAMICS ON A TORUS: NUMERICAL RESULTS
While not as readily experimentally relevant as nonstationary states in a Corbino setup, it is nonetheless instructive to consider the kinds of nonstationary states we can get on a torus.We believe that the results might be applicable to cases where equilibrium density has cor- What happens for intermediate values of α depends on n H and l dw .We nd that if l dw is big enough, then one can transition from the four-wall solution to the two-wall one without ever encountering the nonstationary solutions.As seen in Fig. 2, the potential contour lines that have four-fold symmetry at α = 0 stretch and rotate as α is tuned up, until they become horizontal.
We nd that increasing n H with everything else held xed similarly suppresses the appearance of nonstationary states.At l dw .11,however, nonstationary states do exist.For example, when l dw = .028and n H = 1, we nd nonstationary states for α ∈ [.07, .7].Oftentimes, stationary and nonstationary states (perhaps even of different types, see below) can exist at the same value of α.
Which one we nd in our simulations is both a matter of the initial conditions, as well as the relative size of the basin of attraction for that particular solution.
The nonstationary states we have seen can be broken into three groups.The rst is intimately connected to the solution one gets at small α.The periodic solutions of the rst kind can be thought as the system tunneling from one small α state to the one connected to it by translation by 1 2 x. Figure shows snapshots of such a solution, plotting V (r, t i ) for t i = i T 8 , with i = {0, 1, ..., 7} and T being the period, so that V (r, t + T ) = V (r, t).
In the periodic solutions of the second kind, the timedependent potential is conned to only half of the sample, and has the form V (r− xvt), representing a structure which drifts to the right or to the left with a velocity v.
This behavior is, of course, similar to what we found in the Corbino geometry with a uniform gradient of σ H . Indeed, the observed values of v are, to a good approximation, related to the average value of ∇σ H across a domain in the same way as the drift velocity was related to the gradient of σ H in Eq. ( 14) of Sec.IV.

VI. ANALYTIC RESULTS ON A TORUS
We would like to understand what controls the periods mentioned towards the end of the previous section.The quartic Lyapunov density of Eq. 8 is assumed throughout this section.It turns out to be useful to analyze the stability of the simplest solution, the one with two domain walls perpendicular to y-axis.With the walls xed at y = 1 4 , 3 4 (we keep the aspect ratio κ at 1), and the eld varying along y-axis only, the continuity equation is satised for all values of nonuniformity α, as well as all values of n H and domain-wall thickness l dw .However, the numerical evidence cited above at the very least seems to suggest that, as we decrease α, the basin of attraction for the two-wall solution rapidly shrinks.We will demonstrate below that, below some α c (l dw ), the two wall solution becomes unstable.As we have no analytic solution to the full-blown nonlinear partial-dierential equation, the technique of choice will be to linearize the dierential equation satised by the potential around the two-wall equilibrium and see whether the eigenmodes grow or decay with time for various values of α and domain wall thickness l dw .After linearizing the equation, we will separately consider the eects that nonzero l dw and nonzero α have on the nature of eigenmodes and the corresponding eigenvalues.We will then try to understand what happens when both of the parameters are nonzero.

A. Linearized Continuity Equation
Linearizing around the two-wall solution, we see that changing the eld by δE changes the longitudinal current by Here, E 0 (y) is the eld corresponding to the solution with two horizontal walls.For the quartic Lyapunov density, the form of E 0 (y) is derved in the Appendix for the case of a single domain wall.That result should be applicable in the current case of two domain walls, so long as the separation between them greatly exceed domain wall thickness, l dw .5.The change in the Hall current is δj H = σ H (−δE y x + δE x ŷ).Employing δE = −∇δV , we see that the corresponding divergences turn out to be: Furthermore, since σ H depends only on y, we nd div δj H = (∂ x δV )(∂ y σ H ). The charge continuity equation, aided by the assumption of the local capacitance model, then tells us that Dierent Fourier modes along x-axis will be decoupled, motivating the ansatz δV (x, y) = e at e 2πnxi h(y).Our differential equation then becomes: where q x = 2πn.

B. Case of α = 0
We rst consider the case of uniform Hall conductivity but nonzero domain-wall thickness.Equation ( 15) now turns into The regions immediately surrounding the domain wall and the regions far away from it should be considered separately, turning this problem into a WKB-type computation.The interval [0, 1] is broken into four regions, two bulk domain regions and two thin (extending a few l dw ) regions around the domain walls at 1 4 and 3 4 .The equation is separately solved in each region, then dierent solutions are patched at various boundaries.We will assume a power series expansion for eigenvalue a and eigenfunction h(y): a = a i l i dw , h(y) = h i (y) l i dw .
Note that E 2 0 is symmetric under reection through either the center of domain wall or the center of the domain.That means that we can pick the solution to be symmetric or antisymmetric with respect to reection through the centers of domains or through the domain walls.There are four cases to consider, and we will rst briey mention those two where the solution is antisymmetric under the reection through domains' midpoints.
Going order by order in the expansion, one can show that all contributions are identically zero.We now turn to the more interesting case of midpoint-symmetric solutions.
Solving for h i 's in the bulk and in the boundary layers, and matching these solutions yields the terms in the eigenvalues expansion.We note that since h i in the boundary layers are most conveniently written as function of z = (y−y wall )/l dw , the matching conditions would be Assuming the quartic form (Eq. 8) for the Lyapunov function g(E), we have worked out the rst few leading terms in the expansion of a sym , the eigenvalue corresponding to the mode that is invariant under reection through both domains wall and domain midpoints ) and a asym , the eigenvalue corresponding to the mode that is invariant under reection through domain midpoints but acquires a minus sign when reected though the domain walls (that is, ).As the computation is relatively straightforward, we will omit the details of it here and present the nal result: These expressions illustrate two points.As expected, for high-q x modes, both eigenvalues are negative, simply reecting the fact that modulation of potential on wavelengths smaller than domain wall thickness are strongly suppressed.For relatively small q x , though, a sym > 0.
In particular, considering q x = 2π, time evolution of the eigenmode (see lower right panel of Fig. 5) takes us from the potential that has two horizontal walls to the solution that we know to be stable, the one with two horizontal and two vertical walls, i.e. panel a) of Fig. 1.
Finally, we emphasize that the periodic conditions along y caused the zeroth order (in q 2 x l dw ) term to vanish.In a Corbino setup, the zeroth order term is negative, thus causing both eigenvalues to be negative for thin enough domain walls.This implies that the solution with two horizontal walls would be locally stable in the Corbino geometry.
C. Case of l dw = 0 In the limit of zero domain wall thickness, E 2 0 = 1 everywhere but on the domain walls, where ∂ y E 2 0 diverges.
Setting l dw = 0 in Eq. 15 must then be accompanied by the prescription of how to deal with h(y) on the domain walls.This is not a straightforward task.Instead, we shall take a stab by assuming a simple boundary condition, namely that h = 0 as one approaches the domain walls.This boundary condition will suce to determine the eigenvalues, as ( 15) is now simplied, away from the domain walls, to a second order dierential equation, (Since the positions of domain walls are xed, we have introduced an oset parameter β to study the general case where the domain wall is not necessarily located at an extremum of σ H .) We shall be interested in understanding the behavior of eigenvalue a(α) for small α; in particular, we will be interested in the most physically relevant eigenvalue, with the largest (most positive) real part, as α → 0, a(α) → 0. Note that the requirement that h vanishes at each domain wall essentially decouples dierent domains.The goal is to work out the expansion of eigenvalue a and eigenfunction h(y) in powers of α: a = ãj α j and h(y) = hj (y)α j .One readily sees that real part of eigenvalue (and eigenfunction) only contains even powers of α, while the imaginary part has only odd powers.
Starting with ã0 = 0 and h0 (y) = 1, we can get h1 (y) by integrating twice the dierential equation it obeys, The assumed boundary condition, h 1 ± 1 4 = 0, then xes ã1 .In a similar manner, h2 , ã2 , and higher order terms in the expansion can be obtained.The rst few terms in the eigenvalues expansion are as follows: The next term in the expansion would not be given explicitly here, but we nd it useful to give numerical values of the a i s here: in the case of β = 0, n = 1, one nds ã1 = 25.13,ã2 = −.8987,ã3 = .0149,ã4 = −.0041.We also nd that the rst two terms in the expansion of real part of a, ã2 (β)α 2 + ã4 (β)α 4 < 0 for all β and for α smaller than about 8.36.
We nally comment that making shift β → β + This last result, of course, is a consequence of our assumption that h = 0 at the domain walls, but it is reminiscent of the type III periodic solutions in our numerical results, where charge is moving in one directions in one of the domains, and in the opposite direction in the other domain.In this light, it is not at all surprising that the solutions of this kind show up at fairly large α, where the fact that the domain wall thickness is actually nonzero is not qualitatively important.This gives us reassurance that setting h (y) = 0 on domain walls was a reasonable guess.
D. Case of l dw , α = 0 We are now ready to address the case where the Hall conductivity is nonuniform and the domain walls have nonzero thickness.We restrict the analysis to the cases where α is at most O(1) and we consider the mode with the lowest q x , q x = 2π.The analysis of the previous subsection then implies that, for α/l dw 1, we might expect a(α, l dw ) < 0. The eigenmode then is a decaying one, and we nd that the solution with two horizontal walls is stable.Below some α c (which is a function of n H and l dw ), however, the mode will start growing, and we nd that putting (just about any) perturbation on top of the two-wall solution and time-evolving the resulting guess would take us away from it.Indeed, for α/l dw 1, we expect a(α, l dw ) ≈ a sym (l dw ) > 0 (see Eq. 16).
There is, additionally, another special value of α, α d , below which the eigenvalues are real and unequal to each other, whereas above it the eigenvalues are complex and are conjugate to each other.At that special point, the eigenvalues are, of course, equal to each other.We expect that α d /l dw is O(1).We, however, can't estimate a(α d , l dw ); thus we can't establish whether α c > α d or the other way around.It is plausible, though, that nonstationary solutions would exist for α ∈ [α d , α c ], as, in that case, we would have oscillatory runaway modes.
As for the functional form of the eigenmodes, we expect them to be mostly symmetric between two domains for α < α d .If our guess for the boundary conditions in the previous subsection were a reasonable one, we would expect the modes to have support in one of the domains only for α/l dw 1.

VII. SIMULATION DETAILS A. Numerical Setup & Methodology
The continuity equation was discretized on a triangular lattice of grid points and then evolved in time.For the linear analysis part, the initial guess consisted of a two-wall solution with perturbation that only had one Fourier component in the direction parallel to the walls (along x).Though the precise form of y-dependence of the perturbation did not matter much, we have tried, among other things, a function that is constant on both domains and one that vanishes in one domain but is constant in the other one.In Fig. 4, we have plotted the root-mean-square amplitude A(t) ≡ i ∂ρi ∂t

2
, with the sum being taken over all grid points, as a function of t, for a particular solution.As can be seen in the plot, for small values of t the amplitude increases exponentially over time, while also exhibiting oscillations.
From this portion of the plot, we can infer both the real and imaginary parts of the eigenvalue a l , which governs behavior in the linear regime.At larger times, the exponential growth has saturated, but oscillations persist.From this portion of the plot, we can identify a second frequency, a sat which characterizes oscillations in the saturated regime, as indicated in the gure.
When we were looking for various possible nonlinear solutions, the initial guesses corresponded to the potential at each site being assigned a random number between -0.5 and 0.5.Most of the analysis was carried out on a lattice with mesh size 1 55 , though resolution studies have also been carried out on grids with mesh about half as big and quarter as big.Working with the rougher mesh has restricted us to domain wall thicknesses l dw no smaller than about .03.

B. Linear Regime Results
We rst comment on the solutions and then on the eigenvalues.Recall that the prediction was for eigenmodes to be largely localized to one domain for large α and be symmetric between the two domains around α = 0.And this is indeed what is seen in Fig. 5.The presence of free constants in the perturbative expansions we obtained for the eigenmodes prevents us from making more detailed comparison between theory and simulation.
The eigenvalues inferred from numerical simulations in the case of l dw = .028,n H = 1 are shown in Fig. 6.For α = 0 case, we get an excellent quantitative agreement between the simulation results and the contributions of the rst few terms in perturbative expansion for the eigenvalue a, shown in Eq (16).As predicted, the eigenvalue remains real for a while before turning complex, and these complex ones do indeed come in pairs.
For this particular domain wall thickness, we nd, using the notation of Sec.VI D, α d ≈ 0.07, α c ≈ 0.7.We have found though, that if l dw is suciently big, the order is reversed and one would have α d > α c .Note that the cur-a =1.0 a =.7 a =0 a =.3 vature of a(α) for α > α d is quite close to ã2 of Eq. 17 (solid line of Fig. 5).We can't predict the vertical oset between the simulation results and the solid linethis is a direct consequence of our inability to predict α d and a(α d ) which we have referred to earlier.
As for the imaginary part of the eigenvalue, the agreement between the simulation results and the perturbative expansion is excellent.It is noteworthy that not only is the agreement excellent in the linear regime but even by the time we get to the saturated regime, the agreement is still very good.Indeed, it was the this striking linear (in α) behavior exhibited by the full nonlinear solutions that has motivated the whole expansion enterprise.

VIII. CONCLUSIONS AND EXPERIMENTAL PROSPECTS
We have seen that non-uniformities in the Hall  Real and imaginary parts of eigenvalue a.The left panel shows the real part of a extracted from linear regime simulations on rougher mesh (lled circles), ner mesh (crosses), as well as the rst term in perturbative expansion for the case of l dw = 0 (solid line).The dashed lines connecting crosses and lled circles are merely a guide for the eye.Right panel shows the imaginary part extracted from simulations in the linear regime (with lled circles corresponding to rougher mesh and crosses to ner mesh), the nonlinear regime (open circles), and the rst term in the perturbative expansion for the case of l dw = 0 (solid line).
ionized Si donors in the set-back doping layer, or from density variations that are deliberately imposed on the sample.Density variations of the latter type could be created by applying voltages to one or more gates, located at appropriate distances from the sample.
In our calculations, we have studied explicitly two simple models: a Corbino-type geometry, with a linear dependence of σ H on the distance y form one edge of the sample; and a model with periodic boundary conditions and a variation in σ H that has a sinusoidal dependence on one of the coordinates.In both cases, we have found time-dependent solutions for appropriate choices of the parameters.
The model with periodic boundary conditions may be useful as a rst step to examine the eects of random density uctuations in a macroscopic sample.In this case, we might identify the size L of our model with a region who size is given by the correlation length of the most relevant density uctuations in the physical sample.In our model we found that time-periodic solutions occurred for density uctuations falling within a certain range.We speculate that similar time-varying solutions may occur in the random system, with some regions of the sample trying to generate time-dependent solutions of varying frequencies, and other regions favoring decay to a time-independent solution.If a time-dependent solution occurs in a macroscopic disordered sample, presumably the power spectrum of the dynamic potential uctuations will contain many frequencies, and the domain-wall motion will be chaotic.Since our starting equations are only valid on a length scale greater than the cyclotron radius, of order 1 µm, we would not want to consider intrinsic density uctuations on the scale of the set-back distance, but would restrict ourselves to longer wavelength uctuations of an extrinsic origin.
The Corbino model studied in Section IV suggests a way in which time-dependence may by implemented in a controlled way, which may be most promising for experimental realization.Experimentally, one might aim to create a circularly symmetric Corbino annulus, with an electron density gradient along the radial direction.
The density gradient could be be produced by applying a voltage V G to a front gate, whose distance s from the 2DEG varied linearly with the radial coordinate r.If s is the dielectric constant of the spacer between the gate and the 2DEG, then the imposed variation in the equilibrium electron density would be given by ρ where ρ0 is the mean charge density at V G = 0. for the dielectric constant.If we assume that the electron density changes by a factor 1 + β across the width w of the annulus, then the drift velocity v in Eq. ( 14) has magnitude v = βσ xy /wC.If we assume that the domains have a typical length of order w, then the moving domain structure will lead to a time-dependence of the charge density, at any given point, with a characteristic period given by For the parameters considered here, if β corresponds to a density change of 1%, we nd τ p ≈ 15 ns.
As we have assumed that the domain size is set by the width w of the Corbino sample, to be consistent we should check that w = 100 µm does not exceed the cuto scale L max = ρ/CE c of Eq. 11.There are a few estimates of E c in the theoretical literature 8,18 ; for the typical experimental density n = 3×10 (A1) The high eld conductivity is then given by σ ∞ = σ c (1 + β 2 )/2.In the limit β → ∞, the function g β reduces to the quartic form g quartic (E) given by (8).We note that β = 0. We shall now examine the way in which several different assumptions for g aect the relative Lypunov cost of domain walls of dierent angles γ.We rst demonstrate that, for the quartic form used in our calculations, the Lyapunov energy of a domain wall is proportional to sin 3 γ.We will consider an innite system, with a single domain wall centered around y = 0.The problem is invariant under translations along x-axis, so that electric eld E, and, consequently, longitudinal current j l , are functions of y-coordinate only.Moreover, the x-component of eld is constant: E x = E c cos γ.The boundary conditions on the y-component are E y (±∞) = ±E c sin γ.We are looking for a time-independent solution of a continuity equation.Since both electric eld E and Hall conductivity σ H depend on y only, the Hall current has a vanishing divergence.Therefore, solving div j = 0 would only entail solving dj l y (y)/dy = 0, where the longitudinal current j l is given in Eq. ( 12).Since |E| = E c at y = ±∞, the current density vanishes far away from the domain wall.It follows that j l y must vanish for all y.Then, the equation to solve is This is one of the few nonlinear dierential equations that happens to be exactly solvable.The solution that satises both the boundary conditions and the requirement that the domain wall be centered around y = 0 is (A4) A closed form expression cannot be obtained for the Lyapunov density given by Eq.A1.In that case, we can solve j l y (y) = 0 numerically, then integrate to get δG(γ).Figure 7 shows the Lyapunov cost per unit length of a domain wall as a function of angle γ, normalized by its value at γ = π/2, for quartic Lyapunov density g quartic (E), as well as for g β=.1 (E) (used in Refs.19,20)  and g β=1 (E), which gives σ ∞ = σ c .
Two points need to be made.First, we note that the ratios δG(π/4)/δG(π/2), for both β = 0.1 and β = 1, are smaller than the value 2 −3/2 that is obtained from the quartic form.Thus as stated in Sec.III, Fig. 7 demonstrates that more realistic forms of Lyapunov density would break the degeneracy between having a single horizontal domain wall and a couple of diagonal ones, with diagonal walls (corresponding to π/4 domains) being preferred at least in these examples.That, of course, is helpful, because it is the time evolution of π/4 walls that produces dynamics.
The second point is that the energy cost of very small angle domain walls is actually universal, independent of the the precise form of g(E), if E c and σ c are specied.Indeed, if the angle γ 1, then even the smallest value that (E/E c ) 2 attains, cos 2 γ, is very close to 1. But, by construction, all g(E) look alike in the neighborhood of E = E c .Hence, the energy penalty for a small angle domain wall is given by (A4), irrespective of form of g(E).
, we discuss the phenomenological relation between the dissipative current and the local eld, provide a brief overview of the Lyapunov functional formalism, then use the local capacitance model to arrive at a continuity equation that can be solved for the position-dependent electrochemical potential.We consider the solutions to the continuity equation for two geometries, the rst being Corbino-like, with the potential held xed at two edges along one direction and periodic along the second direction.In a second geometry (torus), periodic conditions along both directions are used.Having established in Sec.II that only variations in the equilibrium charge density contribute to the divergence of the Hall current, in Sec.III we describe, for both geometries, the solutions to continuity equation for uniform equilibrium charge density.Sec.IV contains the central result of the paper; we show that, in a Corbino setup, linearly varying equilibrium charge density would give rise to periodic solutions.In Sec.V, we present the numerical evidence for the presence of nonstationary solutions on the torus.The following section, Sec.VI, is dedicated to gaining analytical understanding of how solutions of previous section arise by analyzing the linearized continuity equation.In Sec.VII, we provide some details on how simulations have been carried out and compare the predictions of the linearized analysis of the previous section with the numerical simulations.We comment on experimental prospects in our concluding section, Sec.VIII.A discussion of some alternative assumptions for the Lyapunov function, and their consequences for the relative costs of dierent domain-wall congurations, is given in an Appendix.

Fig. 1 InFigure 1 :
Fig. 1 In the simpler case of all-periodic boundary conditions, the ground state, shown in panel a) of Fig. 1, would have two horizontal and two vertical walls, with electric eld forming angle δ (tan δ = κ) with horizontal walls.Of course, there will be a degeneracy of these states due to translational invariance; the pattern can be displaced by a arbitrary amounts in the horizontal and vertical directions with no change in Lyapunov energy.The situation is more complicated in the case of the Corbino-like geometry, where we assume periodic boundary conditions in the x-direction, but the top and bottom boundaries are grounded.Results actually depend on the sample aspect ratio κ.For κ > 0.5, the case illustrated in panel b) of Fig. 1, the lowest energy state has a single horizontal domain wall at y = κ/2.(There is actually a doublet of ground states, related to each other by a global ip of the signs of the electric eld.)For κ = 0.5, there are two degenerate ground states, the rst again consisting of the single horizontal wall centered at the middle of the domain, with the second, shown in panel c) of Fig. 1, having four domain walls forming an angle of π 4 with the top/bottom boundaries, and the electric eld pointing along ±ŷ in the domains adjacent to the top and bottom boundaries, and along ±x in the domains that are fully in the bulk.This second solution can, of course, be displaced in the horizontal direction by an arbitrary amount, without penalty.For the case of κ < 0.5, there are an uncountable innity of ground states, which may be thought of as the combination of the single horizontal wall and the diagonal wall states, as illustrated schematically in panel d) of Fig. 1.The interactions between the walls will most likely break this continuum down to a nite set of states, but such investigation lies outside length that is much smaller than the sample size, and that, in such cases, one might get away with modeling a smaller piece of sample with periodic boundary conditions, rather than all of it with Corbino or Hall-bar boundary conditions.We consider here a sample of size L x = L y = 1, with periodic boundary conditions.As described in Sec.III, the equilibrium state for the case of uniform Hall conductivity has two horizontal and two vertical domain walls with eld-angle γ = π/4.In our analysis, we shall assume that the spatially-varying part of the Hall conductivity has the simple form σ H (r) = α sin 2πn H y, where n H is a positive integer.Our interests lie in trying to understand how the solutions depend on the non-uniformity parameter α, for dierent n H and domain wall thicknesses l dw .As we tune nonuniformity α away from 0, the Hall current develops a nonzero divergence in the equilibrium state.We nd that, so long as the α is small enough, the four-wall solutions can accommodate such nonuniformity by bending.We may imagine that the gradient of the σ H leads to a horizontal drag force", proportional to ẑ × ∇σ H , which can be counteracted, for small α, by the restoring force" produced by a distortion of the domain structure, away from the shape that minimizes the Lyapunov functional.By contrast, for large values of α, we nd only solutions which are translationally invariant along the x-direction, and have have two or more horizontal domain walls with γ = π/2.The divergence of Hall current trivially vanishes for such states.

Figure 2 :
Figure 2: Numerically-obtained electochemical potential V (r) for domain wall thickness l dw = .11,nH = 1.All V (r) shown are stationary.The transition between four-wall solution at α = 0 and the one with two walls for α = 0.3 is continuous; that is, it is not interrupted by the appearance of nonstationary solutions at the intermediate values of α.

Finally, the periodic
solutions of the third type, seen only for n H = 1, can be thought as the combination of two solutions of type II: in one half of the sample, the potential moves to the right, while in the other half it moves to the left, with more complicated processes happening at the boundary between the two regions.As examples of these behaviors (all but the rst one not shown in the gures), for α = .1,l dw = .028,we have found periodic solutions for n H = 1, 2, 3, with the periods being 2.49 (n H = 1, type I), 1.90 (n H = 2, type II), 1.64 (n H = 3, type I), and 1.21 (n H = 3, type II).

Figure 3 :
Figure 3: Numerically obtained electrochemical potential V (r, t) for domain wall thickness l dw = .028,nH = 1, α = 0.1 at dierent points during its period.The solution here tunnels from one α = 0 solution to the one related to it by translation by x/2.
from one domain to the other.This leaves the real part of eigenvalue unchanged while ipping the sign of the imaginary part.It follows that eigenvalues come in conjugate pairs, and that corresponding eigenfunctions must identically vanish in one of the domains.

Figure 4 :
Figure 4: Time trace of A ≡ r P i ˛∂ρ i ∂t ˛2, with the sum being taken over grid points, can be used to extract both the real and imaginary parts of frequencies.The transition from linear regime (where the real part of the frequency ω l is nonzero) to the saturated one (where the real part of ωsat vanishes) is apparent.

Figure 5 :
Figure 5: Plots of absolute value of eigenmodes of the linearized continuity equation for all-periodic boundary conditions, l dw = .028,and nH = 1.The eigenmodes are deviations from the solution with two horizontal walls; thin dashed lines indicate the location of these domain walls.Plus and minus signs indicate whether the deviation is positive or negative, while white color corresponds to zero deviation.For suciently big values of the nonuniformity α, all of the eigenmode's support is within one domain.On the other hand, for α = 0, the (unstable) mode is symmetric between two domains.
conductance σ H , produced by inhomogeneities in the equilibrium electron density, can lead, under certain circumstances, to time-dependent domain patterns in the microwave-induced Zero-Resistance states seen in 2DEGs.Inhomogeneities in the equilibrium density may result from disorder, such variations in the density of

1
was used in numerical calculations whose results are shown in Figs. 1 and 2 of Ref. 19 and Figs. 4 and 5 of Ref. 20

Figure 7 :
Figure 7: Normalized Lyapunov cost of domain wall per unit length, δG(γ)/δG(π/2), as a function of the angle γ between the eld in the bulk of the domain and the domain wall.Solid line is the exact result for quartic Lyapunov density.Triangles correspond to the Lyapunov density g β of Eq.A1 with the value β = 0.1, used in our previous work (Refs.19,20), while crosses show results for β = 1.
For example, if s were to vary between 1 µm and 2 µm, from the inner radius to the outer radius of the annulus, and if s ≈ 13, one would need a gate voltage V G = 0.09 V to produce a density dierence of 3.5 × 10 9 cm −2 across Let us consider a Corbino sample, whose width w is much smaller than its circumference L, say .1 mm and 10 mm respectively.In this regime, the curvature of the annulus is not important, and the periodic states should resemble those described in Sec.IV.To apply the calculations of Sec.IV, we estimate the capacitance per unit area C = 0 /d by equating the distance d to 1.5 µm, the average distance to the gate, and use s = 13 even if the spacer thickness is uniform by depositing several separately-contacted gates on top of the spacer and putting dierent voltages on them.We can use the above numbers to estimate the characteristic frequency of the charge oscillations that might be observed.
in a system with disorder, as it did in our model with periodic boundary conditions, where the nonuniform drift velocity v is in competition with the tendency of the Lyapunov function to favor an incompatible domain structure.The dimensionless parameter α, which controlled the transition from time-independent to time-dependent solutions in Sec.V, may be written as α = βσ xy /σ c , where β is here the fractional change in the electron density between its minimum and maximum values.The value of σ c depends on details of the microscopic model, but it is expected generally to be of the order of the longitudinal conductivity σ xx in the absence of microwave radiation.For the sample described in Ref.2, in a eld of ∼.11T, where the most prominent ZRS was seen once the sample was subjected to microwaves, we estimate a value for σ xx of order 5 × 10 −5 S, while σ xy ≈ 5 × 10 −3 S. It is the simplest possible function that has the following desired properties: it is an analytic function of E 2 , and hence an an isotropic analytic function of E; it has an absolute minimum at a non-zero value of E, and a local maximum at E = 0.The function is then determined by the position E c of the minimum and the curvature σ c at E = E c .However, if the quartic function is extended to values of E much larger than E c , the longitudinal dierential conductivity, d 2 g/dE 2 is seen to increase without limit, going eventually as E 2 , which does not seem reasonable.It might be more reasonable to assume that at large values of E, the longitudinal dierential conductivity should approach a constant value σ ∞ .Although we are not directly interested in the behavior of g(E) for E > E c , the assumption that g is analytic does introduce a link between this behavior and the behavior at E < E c .In particular, if we assume a nite value for σ ∞ , this implies that d 2 g/d(E 2 ) 2 → 0 for large E 2 , while it is positive at E = E c , which at least suggests that d 3 g/d(E 2 ) 3 is likely to be negative at E = E c .If we assume that g has the form of a cubic polynomial in E 2 for E < E c , this would imply that g(0) − g(E c ) is larger than in case of a pure quartic form, for the same values of E c and σ c .In previous work, (Refs.19,20),we assumed the fol- 11cm −2 , the highest estimate would be on the order 20 V/cm.Experimentally, while E c has not been measured, some of the results of Refs.3,4 can be interpreted to imply E c ∼0.1 V/cm.Cautiously, using the higher estimate for the critical eld and the same estimate for capacitance per unit area C as above, we arrive at L max ≈0.6 cm.So, in a Corbino sample that is 100 µm wide, the size of the domain should indeed be set by the the Corbino ring's width.By varying the sample width, the spacer thickness, or the density gradient the time scale of interest, τ p , can be tuned over a wide range.w is given by E c w.For a domain of width w = .1mmandEc = .1V/cm,thiswouldimply voltage varying by 1mV.The parameter σ c in the Lyapunov function g(E) does not enter directly in our estimate for the oscillation period in the Corbino geometry.However, it should play a role 2 .