Hydrodynamic and Ballistic AC transport in 2D Fermi Liquids
Mani Chandra Gitansh Kataria Deshdeep Sahdev Ravishankar Sundaraman



mani@quazartech.com
Research Division, Quazar Technologies, Sarvapriya Vihar, New Delhi, India, 110016


Research Division, Quazar Technologies, Sarvapriya Vihar, New Delhi, India, 110016


Research Division, Quazar Technologies, Sarvapriya Vihar, New Delhi, India, 110016


Department of Materials Science and Engineering, Rensselaer Polytechnic Institute, Troy, NY 12180



Electron transport in clean 2D systems with weak electron-phonon (e-ph) coupling can transition from an Ohmic to a ballistic or a hydrodynamic regime. The ballistic regime occurs when electron-electron (e-e) scattering is weak whereas the hydrodynamic regime arises when this scattering is strong. Despite this difference, we find that vortices and a negative nonlocal resistance believed to be quintessentially hydrodynamic are equally characteristic of the ballistic regime. These non-Ohmic regimes cannot be distinguished in DC transport without changing experimental conditions. Further, as our kinetic calculations show, the hydrodynamic regime in DC transport is highly fragile and is wiped out by even sparse disorder and e-ph scattering. We show that microwave-frequency AC sources by contrast readily excite hydrodynamic modes with current vortices that are robust to disorder and e-ph scattering. Indeed, current reversals in the non-Ohmic regimes occur via repeated vortex generation and mergers through reconnections, as in classical 2D fluids. Crucially, AC sources give rise to strong correlations across the entire device that unambiguously distinguish all regimes. These correlations in the form of nonlocal current-voltage and voltage-voltage phases directly check for the presence of a nonlocal current-voltage relation signifying the onset of non-Ohmic behavior as well as also for the dominance of bulk interactions, needed to confirm the presence of a hydrodynamic regime. We use these probes to demarcate all regimes in an experimentally realizable graphene device and find that the ballistic regime has a much larger extent in parameter space than the hydrodynamic regime.


Hydrodynamic and Ballistic AC transport in 2D Fermi Liquids
    Ravishankar Sundaraman
    Accepted 23/03/2018
===========================================================





§ INTRODUCTION


Charge transport in conductors is typically dominated by electrons scattering against phonons and defects, resulting in momentum relaxation (MR) over a time scale τmr∼ 10^-14-10^-15 s and a corresponding length scale lmr∼ 10^-2-10^-3 μm. In comparison, the time scale τmc∼ 10^-12 s of momentum-conserving (MC) scattering (primarily due to electron-electron (e-e) scattering[Note that electron-electron Umklapp scattering relaxes momentum, while small-angle electron-phonon scattering at low temperatures does not; we work with momentum-conserving and momentum-relaxing time-scales to avoid this ambiguity.]) is negligible <cit.>, with the corresponding length scale lmc∼ 1 μm. Electrons scattering against phonons or defects give rise to the characteristically diffusive Ohmic transport. However, in clean systems with weak electron-phonon scattering, the length scale of momentum relaxation can approach the device scale[We are considering devices with dimensions slightly smaller than lmr but much larger than the phase coherence length scale so as to preclude interference effects.] (e.g., lmr∼ a few μm in graphene <cit.>). Electrons in this ballistic regime scatter predominantly against device boundaries. If e-e scattering is then increased upon adjustment of experimental conditions such as carrier concentration and/or temperature such that lmc is made sufficiently smaller than the device scale, electron transport can be hydrodynamic. In this regime of slow MR and fast MC scattering[Our discussion pertains to Fermi liquids. Hydrodynamics is then a regime of weakly interacting (so that quasiparticle excitations are well-defined), but fast MC interactions.], electrons are expected to move collectively as in a fluid whose dynamics are governed by macroscopic conservation laws, i.e., the Navier-Stokes equations. Candidate materials for hydrodynamic charge transport require large lmr. These include (Al,Ga)As heterostructures <cit.>, GaAs <cit.>, PdCoO_2 crystals <cit.>, WP_2 <cit.> and graphene <cit.>.

Signatures and novel consequences <cit.> of a hydrodynamic regime in a Fermi liquid (also referred to as a viscous regime) have recently been the focus of intense activity, since the regime arises from enhanced MC interactions. Calculations based on fluid models indicate that DC charge transport by viscous electrons obeys a nonlocal current-voltage relation and produces a negative nonlocal resistance <cit.>, in sharp contrast to the local current-voltage relation of the Ohmic regime; clearly, interactions do not simply renormalize the conductivity. Perhaps even more striking is the possibility of generating current vortices in a hydrodynamic flow <cit.>, which are absent in an Ohmic regime. 

Distinct from the hydrodynamic regime is the ballistic regime, which arises when both MR and MC interactions are rare. The absence of MR interactions implies that this regime also conserves momentum in the bulk. This results in a degeneracy in DC transport between the hydrodynamic and ballistic regimes. We show that both have negative resistances (fig. 2c) and, remarkably, even current flows ordered into vortices (fig. 3). Thus, a negative resistance only indicates a non-Ohmic regime (hydrodynamic/ballistic) resulting from the absence of MR interactions, not the dominance of bulk MC interactions. Experiments <cit.> currently resolve the hydrodynamic-ballistic degeneracy by changing temperature, carrier concentration or device/contact geometry, and fitting against an expected hydrodynamic or ballistic response to the altered conditions<cit.>.

In this paper, we introduce and exploit AC transport as a powerful technique for studying non-Ohmic regimes. We show that each of hydrodynamic, ballistic and Ohmic regimes can be directly identified using spatiotemporal correlations in AC transport without any change in experimental conditions. Further, we show that vortex formation can be accessed in AC transport with much less fine-tuning compared to DC transport. A snapshot of the key results is presented in fig. 1.

We illustrate each of these points by considering graphene (for concreteness) in a generic device geometry. We show that DC sources fail to generate vortices in the presence of even a modest amount of MR scattering (τmr∼ 10 ps), and show that these are easily generated by switching to AC sources of experimentally accessible frequencies (∼ GHz). In fact, vortex dynamics are crucial for AC transport to proceed both in the hydrodynamic and ballistic regimes. The associated current reversals occur via repeated vortex formation and mergers through reconnection and annihilation; mechanisms seen widely in classical two-dimensional fluids.

We obtain our results by solving for the dynamics of quasiparticles using a deterministic high-resolution numerical scheme which converges much faster than statistical particle methods. The kinetic approach naturally gives rise to all three transport regimes by simply varying τmr and τmc. In particular, the ballistic regime cannot be accessed by effective fluid models and requires solving for the full time-dependent non-equilibrium distribution function over the entire device. This is made possible using an efficient implementation of our numerical scheme in the package bolt <cit.> that exploits the massive computational power of GPU computing clusters.

The paper is structured as follows. In <ref>, we describe the kinetic model which we solve for the graphene device outlined in <ref>. We then study DC transport for this device in <ref> and highlight the challenges of using DC signatures to discern the different regimes. Next, we proceed to AC transport and present vortex dynamics and the distinct spatiotemporal correlations of the regimes in <ref>. In <ref>, we exploit these correlations to map all the regimes in the {τmc, τmr} phase space for the device and contact geometry considered. Finally, we provide an argument in <ref> to understand why vortex ordering occurs in both the hydrodynamic and ballistic regimes and then conclude in <ref>. Additionally, in appendix <ref>, we show how our kinetic calculations incorporate the effect of long-range Coulomb interactions and also show how they appear in fluid models in appendix <ref>. The appendices also discuss two different approximations used to compute current-voltage characteristics and verify that they are consistent with each other.





§ KINETIC MODEL
 

We consider graphene (a) well above the charge neutrality point, where quasiparticle excitations are well-defined, and (b) over length scales (∼μ m) where quantum interference effects are washed out. Transport is then described by the Boltzmann equation that governs the evolution of a charge carrier distribution f(𝐱, 𝐩, t) in the 4-dimensional phase space of spatial 𝐱≡ (x, y) and momentum 𝐩≡ (px, py) coordinates,

    ∂ f/∂ t +  𝐯·∂ f/∂𝐱    = - f - fmr_0/τmr -  f - fmc_0 /τmc

where the velocity 𝐯 = ∂ℰ/∂𝐩 and ℰ(𝐩) is the band energy. For doped graphene in the upper band, ℰ(p) = vF  p 𝐯 = vF 𝐩̂, where vF≈ 10^6 m/s = 1 μm/ps is the Fermi velocity. The terms on the right are the MR and MC collision operators, parametrized in a relaxation time approximation by τmr and τmc respectively. In writing (<ref>), we are working in the μ/(k_BT) ≫ 1 regime where the lower band is not needed. 

Electrical transport is set up through current injection at device boundaries which gives rise to a chemical potential gradient. We show in appendix <ref> that backreaction -e𝐄·∂ f/∂𝐩 (e≡ charge magnitude) of the self-consistent electric fields 𝐄 is incorporated at linear order, even though the term is not explicitly present in (<ref>).

The collision operators relax f(𝐱, 𝐩, t) to stationary and drifting local Fermi-Dirac distributions, fmr_0(μmr, Tmr) and fmc_0(μmc + 𝐩·𝐯d, Tmc) respectively. The spatiotemporal Lagrange multipliers {μmr(𝐱, t), μmc(𝐱, t)}, {Tmr(𝐱, t), Tmc(𝐱, t)} and 𝐯d(𝐱, t) are needed for charge, energy and momentum conservation respectively. These are solved for by imposing the matching conditions, 

    ⟨ fmr_0 ⟩ = ⟨ f ⟩ = ⟨ fmc_0 ⟩
    ⟨ℰ(p) fmr_0 ⟩ = ⟨ℰ(p) f ⟩ = ⟨ℰ(p) fmc_0 ⟩
    ⟨𝐩f⟩ = ⟨𝐩fmc_0⟩

where ⟨⟩ = N/(2 πħ)^2 ∫ d^2p and N = 4.

The model thus evolves the four-dimensional electron distribution function f(𝐱, 𝐩, t) by additionally solving the six Lagrange multiplier constraints (<ref>, <ref>, <ref>) at every time step. The free parameters in the model are τmc, τmr. We treat these as numerical inputs in units of picoseconds and do not invoke any functional dependence on thermodynamics quantities such as temperature and number density. However, the model incorporates temperature smearing of the Fermi surface, although we find that all results are independent of this effect. We numerically integrate this computationally expensive system on a GPU cluster using bolt <cit.>, a fast, massively parallel high-resolution solver for kinetic theories based on a finite volume method to achieve 𝒪(Δ x^2, Δ p^2, Δ t^2) accuracy, where Δ x, Δ p and Δ t are the sizes of discrete elements in real space, momentum space and time respectively.



 §.§ Current-Voltage Relationship
 

A key quantity of interest is the  current-voltage relationship. After f(𝐱, 𝐩, t) has been solved for, the current 𝐣(𝐱, t) is easily computed using 𝐣(𝐱, t) = -e⟨𝐩 f ⟩. However, to obtain the voltage V(𝐱, z, t), one needs to solve the 3D Poisson equation,

    ∇·(ϵr∇ V) = 4π e n δ(z - z_0)

where n(𝐱, t) = ⟨ f ⟩ is the 2D charge carrier density, z_0 is the location of the 2D sample in the perpendicular direction and ϵr(𝐱) is the dielectric function. 

To proceed, we assume a graphene field-effect transistor geometry with a dielectric substrate ϵr of thickness d ≪ L, where L ∼μm is a lateral device scale. The mean carrier density n_0 can be set to a desired level by applying a backgate voltage V(𝐱, z=0) = Vg on the substrate (located at z=0) beneath the dielectric. This is given by,

    n_0 = C Vg/e

where C = ϵr/(4π d) is the capacitance per unit area.

Further, we require the relationship between charge inhomogeneities and in-plane voltage fluctuations. For the field-effect geometry, an analytic solution of the 3D Poisson equation in the 2D plane of the device has been computed by Tomadin and Polini <cit.>. To 𝒪(d/L), this solution simplifies to

    Δ V ≈ -e/CΔ n

where Δ V and Δ n are the in-plane voltage and charge density differences respectively.  Note that while (<ref>) and (<ref>) appear similar, they relate different quantities. The approximation (<ref>), referred to as “gradual channel” in <cit.> and “local capacitance” in <cit.>, thus relates f obtained by solving (<ref>) to the voltages measured. As explained in appendix <ref>, this computation is a post-processing step and is not required for the evolution of (<ref>) to linear order.









§ SETUP
 

We consider a 5 μm × 10 μm graphene device, with drive contacts 1 μ m wide at the center of the left and right edges (see fig. 1a). We assume an electron density of n = 10^12 cm^-2 , and ideal Ohmic contacts i.e., the Fermi level of the contact metal is the same as the electron chemical potential in graphene (at the chosen carrier density). 
At the contacts, we impose Dirichlet boundary conditions that implement a current source/sink, with the distribution at both contacts set to a shifted Fermi-Dirac (fmc_0) with the drift velocities 𝐯Ld = 𝐯Rd = (v(t), 0), where L and R denote the left and right contacts respectively and v(t) is a time-dependent magnitude. For DC calculations, we set v(t) = v_0 = 10^-4  vF, corresponding to a current injection of ∼ 0.1 μ A. For AC calculations, we set v(t) = v_0 sin(2 π f t), where f is the source frequency. We present results for f = 10 GHz, but they are valid over a wide range of frequencies as discussed later.

On the device boundary outside of the contacts, we impose perfect reflection on the electrons (specular scattering). This corresponds to “free-slip” boundaries in the parlance of fluid models, as opposed to “no-slip" boundaries. The question of what boundary conditions are correct is an open one, but there is increasing evidence in support of free-slip boundaries because of the suppression of the Gurzhi effect <cit.> in graphene, which becomes dominant only in the presence of no-slip boundaries <cit.>.



§ DC TRANSPORT
 

We first examine signatures of a hydro regime in DC transport. For fast MC and very slow MR interactions, {τmc, τmr} = {0.2, 50} ps τmr/τmc = 250, we get current vortices flowing against the drive (fig. 2a); the distinctive features of a fluid. Further, the nonlocal resistance computed using voltage measured far from the driving leads divided by the injected current is negative (fig. 2c), as shown using fluid models <cit.>. Now consider {τmc, τmr} = {0.2, 10} ps τmr/τmc = 50, still expected to be deep in the hydro regime. However, the vortices no longer appear (fig. 2b). We have verified their absence everywhere in the domain down to 25 nm (∼ 0.1 × e-e
mean free path); well below the length-scale at which a
hydrodynamic description is expected to apply. In addition, the nonlocal resistance is negative
only locally, within 0.5 μm of the driving leads (fig. 2(b,c)). A further reduction in τmr to 5 ps (τmr/τmc=25) leads to a disappearance of this negative resistance as well (fig. 2c). 



 §.§ Hydro-Ballistic Degeneracy
 

Consider now the ballistic regime set by the parameters {τmc = ∞, τmr≃ L/v_F = 5} ps, i.e. zero MC interactions and weak MR interactions. The ballistic flow is qualitatively indistinguishable from a hydrodynamic flow (both shown in fig. 3), and produces a negative resistance all along the edge (fig. 2c). 

Thus, we see that not only are DC sources inefficient at exciting hydrodynamic behaviour in parameter regimes that are clearly MC dominated but also fail to distinguish between ballistic and hydrodynamic regimes. Note that the specific thresholds between regimes we report here are
for the device and contact geometry shown in fig. 1a, and these would vary for other geometries. In particular, the geometry proposed by Torre et al. <cit.> to investigate hydrodynamic effects has contacts on the same side of the device and is much more conducive to vortex generation in DC transport <cit.>. Our choice of geometry, adopted from Levitov & Falkovich <cit.>, has been made so as to demonstrate vigorous vortex production even in the least favorable setup, as we shall see in the next section.





§ AC TRANSPORT
 

The situation improves dramatically if we use an AC source with frequency f ≪ vF/L = 200 GHz. As long as lmr≡τmrvF≳ L, MC scattering (either specular boundary or bulk e-e) excites collective modes involving vortices that are continuously generated and destroyed at the rate f.
This condition in AC is far more enabling for vortex generation in the hydrodynamic regime (lmc≪ L) than in DC which requires D ∼ vF√(τmcτmr)/2 ≳ L/(√(2)π), where D is the vorticity diffusion length <cit.>. Therefore, vortex generation occurs in AC even for parameter regimes where current streamlines appear distinctly Ohmic in DC.
For example, with {τmc, τmr} = {0.2, 10} ps,
lmr≈ 2· L enabling vortices in AC as shown in fig. 4(b-d),
but D ≈ 0.6· L/(√(2)π) resulting in Ohmic-like DC transport in fig. 2(b).
Indeed, vortices form in AC even for the marginal case of
{τmc, τmr} = {0.2, 5} ps,  where lmr≈ L
(fig. 1a).



 §.§ Vortex dynamics
 

Fig. 4(a-d) show the flow structure of the hydro mode excited by a 10 GHz source for {τmc, τmr} = {0.2, 10} ps. As the source is about to change sign near the half-cycle, vortices form symmetrically at the left and right contacts resulting in a quadrupolar mode (fig. 4b). These vortices, with the same sign of vorticity, grow and merge through reconnection (fig. 4c); as in 2D classical fluids. There is now a dipolar mode in the device with vortices in the top and bottom halves having opposite signs of vorticity (fig. 4d). These then annihilate in the middle, and allow the flow to reverse. In contrast, Ohmic AC transport proceeds through wave-fronts (fig. 4g) that originate from the contacts and travel into the device (fig. 4(h)). We note that ballistic AC transport also proceeds through vortex dynamics, albeit with an altered choreography; vortices form at the top/bottom boundaries and move inwards.



 §.§ Spatiotemporal correlations
 

The persistent time-dependence of an AC source produces several useful spatiotemporal correlations. We consider the phase ϕ(Id, V) between a nonlocal voltage V(t) measured by contacts on the edge (top of fig. 1a), and the current source Id(t) (center of fig. 1a). An Ohmic regime is defined by a local current-voltage relationship; an injected current causes local changes in voltage (∝ Id) which then propagate into the device (fig. 4(g,h)). Therefore, the measured V lags the source Id (fig. 5a) and the phase ϕ is negative.

A transition from a negative to a positive phase (fig. 5a) signifies the breakdown of a local current-voltage relation and the onset of a non-Ohmic regime (fig. 5b) with a non-local current-voltage relation (see <cit.> for hydro).
This arises whenever MR interactions are weak, which is a necessary
but not sufficient condition for a hydrodynamic regime.
Hydrodynamic transport also requires strong MC interactions that impose a local equilibrium
in the bulk, unlike specular scattering at the boundary in the ballistic regime.

This key difference in the ballistic regime
is directly captured by the two-point correlation
ϕ(V_1, V_2) = ϕ(Id, V_2) - ϕ(Id, V_1),
with a V_1(t) measured closer to the source (schematic shown in fig. 1a).
When voltage gradients flow through bulk interactions
(MC or MR), V_2 must lag V_1 and ϕ(V_1, V_2) < 0.
This is evident from the negative slope of ϕ(Id, V)
in both Ohmic and hydrodynamic regimes (fig. 5b).
The ballistic regime which only has boundary scattering
exhibits a positive slope in ϕ(Id, V) (fig. 5b),
so that ϕ(V_1, V_2) > 0 and V_2 closer to the
top boundary now leads V_1.



 §.§ AC vs DC signatures


The ϕ(Id, V)>0 diagnostic of non-Ohmicity is similar to negative resistance in DC, with both indicating a nonlocal current-voltage relation. However, ϕ(V_1, V_2) has no analogue in DC; the slope of resistance does not show any robust pattern (fig. 2c). In addition, the spatial extent of the phase correlations is much greater in a hydrodynamic regime. For {τmc, τmr} = {0.2, 10} ps, the negative resistance is confined to within 0.5 μm of the driving leads (fig. 2c). In contrast, the phase at 10 GHz for even {0.2, 5} ps is positive over the entire edge (fig. 5b), and there is no negative resistance anywhere for these parameters in DC (fig. 2c).





§ PHASE DIAGRAM
 

We now use the phases ϕ_1 ≡ϕ(Id, V_1), ϕ_2 ≡ϕ(Id, V_2) and ϕ_12≡ϕ(V_1, V_2), measured using contacts shown in fig. 1a, to map all regimes in the {τmc, τmr} parameter space for our device. We consider a wide range of these parameters, all of which satisfy τmc/τmr < 1; parameters where the hydrodynamic regime could potentially arise.

Fig. 6a shows ϕ_2>0 for all τmr≳ 3 ps (same for ϕ_1, not shown), in accordance with the τmr≳ L/vF = 5 ps condition for non-Ohmic AC transport. Further, we find ϕ_12 < 0 for all parameters in fig. 6a, indicating an Ohmic to hydrodynamic transition. To now transition from hydrodynamic to ballistic, we vary τmc for a fixed τmr=5 ps in fig. 6b. We indeed find non-Ohmicity (ϕ_1, ϕ_2 > 0) over the entire range of τmc but bulk scattering (ϕ_12<0) only for τmc≲ 1.7 ps ≈ 0.3 L/vF; consistent with the lmc≡τmcvF≪ L requirement of a hydrodynamic regime. Thus we need both τmr≳ 3 ps and τmc≲ 1.7 ps, quantifying the weakness and strength of MR and MC interactions needed for a hydrodynamic regime in our device geometry.

Both conditions on MR and MC interactions are within reach of current devices<cit.>, with the requirement of weak MR further mitigated by adjusting the lateral scale such that L ≲ lmr≡τmrvF. As discussed in <ref>, this condition is much weaker than that required for vortices in DC where the the minimum τmr scales as L^2. The critical τmr in both AC and DC is shown in fig. 7 for a fixed τmc = 0.2 ps. This τmc is well within the limit ≲ 0.3 L/vF required for a hydrodynamic regime for the smallest device considered (L = 1.8 μm), once the criterion on τmr is satisfied. The regime boundaries for the device geometry are summarized in fig. 1d. Notably, the parameter space of the ballistic regime is much larger than that of the hydrodynamic regime.







 §.§ Consistency check
 

We additionally consider the transition from Ohmic to ballistic by setting τmc = ∞ and varying τmr. A correct identification of this transition requires a deep consistency within the correlations; each of ϕ_1, ϕ_2 and ϕ_12 must change sign from negative to positive at the same τmr. Fig. 8 of the phases versus τmr indeed shows this at τmr = 2.3 ± 0.1 ps, with the small spread due to finite width (1 μm) contacts placed a finite distance (1 μm) apart in our setup.



§ VORTICES IN HYDRO AND BALLISTIC REGIMES
 
The hydro regime is interaction dominated whereas the ballistic regime occurs in the near absence of interactions. However, the current flow patterns in both these regimes are strikingly similar. Both regimes have flows organized into vortices, and therefore have a nonzero vorticity (ω≡∇×𝐯d). Note that the vorticity is identically zero in the Ohmic regime since 𝐯d∝∇ V.

Let us ask the question, can we deduce flow ordering in both hydrodynamic and ballistic regimes without solving the governing equation (<ref>)? We provide here an argument using Landau theory that the ordering principle for the flow is the same in both regimes.

To proceed, we require an energy functional constrained by the symmetries of the system. Consider no MR interactions (τmr = ∞) and fast MC interactions resulting in vanishing electron viscosity ν∼τmc→ 0. We can then invoke a well-known result that applies to inviscid 2D fluids, the conservation of “enstrophy”,

    F[𝐯d] = ∫ d^2r |∇×𝐯d|^2

where 𝐯d is the local carrier drift velocity. The conservation of enstrophy follows from momentum conservation of a fluid element in two-dimensions <cit.>.

We interpret F[𝐯d] as a Landau energy functional with a local order parameter ω. A vortex ordered state then has F > 0 with F = 0 denoting no ordering. A variation of F with respect to the velocity field 𝐯d gives,

    δ F = ∮ d𝐬· (δ𝐯d×ω) + ∫ d^2r δ𝐯d·(∇×ω)

The first term on the right vanishes with δ𝐯d = 0 on the boundaries, corresponding to the imposition of boundary conditions. Demanding that the variation of F vanish,

    δ F = 0    ∇×ω = 0

Now applying the curl operator on (<ref>),

    ∇^2 ω = 0

This is in fact the time-independent fluid momentum conservation equation (<ref>) in the limit τmr→∞, after an application of the curl operator. Clearly, (<ref>) admits solutions with F > 0 depending on the boundary conditions, i.e, device and contact geometry. Therefore, we see that the conservation of enstrophy (<ref>) allows for flow ordering in the form of vortices.

Now we note that the enstrophy (<ref>) is also conserved in the ballistic regime and thus we expect the same ordering in this regime as well. Indeed, the ballistic regime preserves an infinite number of invariants because the momentum of every carrier is conserved, thus conserving the momentum of a macroscopic fluid element. In contrast, the hydro regime only conserves the momentum of the fluid element; the momenta of individual carriers is thermalized by MC interactions.

The conservation of enstrophy (<ref>) implies that any local vorticity injected from the boundaries persists in the device. The device and contact geometry adopted here is an example where vorticity is zero globally but nonzero locally, near the drive contacts. Since the invariant measure (<ref>) is positive definite, this geometry allows for vortex formation. A counterexample is a wire geometry where the injected vorticity is exactly zero. Note that the arguments presented in this section only show that flow ordering is possible provided a conducive geometry and do not inform on the interaction timescales required.



§ DISCUSSION AND CONCLUSION
 

The hydrodynamic transport regime in Fermi liquids has long been theorized, but its robustness has always been in question. In particular, the specific requirements on defect densities, and e-ph, e-e coupling needed to access this regime in various device geometries have been open issues. Indeed, we find that the regime is very delicate in DC transport, requiring unreasonably low momentum-relaxing scattering or specially designed geometries in order for its effects to manifest themselves. An even bigger concern is its high degeneracy with the ballistic regime in which we find nearly identical current vortex structures and negative resistances.

We have shown that both issues can be resolved by switching to AC transport. First, the requirement on momentum-relaxing scattering for the hydrodynamic regime is greatly mitigated, with the criterion simply being lmr≳ L as compared to lmr≳ (0.2/lmc) L^2 in DC. Hydrodynamic modes rich in vortices are therefore easily excited in AC transport, with current reversals occurring through intricate vortex dynamics. Note that the condition on lmr in AC is independent of lmc whereas it becomes increasingly stringent in DC as lmc→ 0!

Second, there exist strong spatiotemporal correlations in AC transport that are able to uniquely and directly identify all transport regimes (as summarized in fig. 1c). One of these correlations, the phase between a nonlocal voltage and the current drive tests for the locality of the current-voltage relation and can thus differentiate between Ohmic and non-Ohmic regimes. Further, the phase between two nonlocal voltages tests for the presence of bulk interactions and can therefore distinguish between the hydrodynamic and ballistic regimes. These correlations extend over the entire device, allowing for great flexibility in selecting the positions of probe contacts in experiments. In contrast, we have shown that negative resistances in hydrodynamic DC transport are tightly localized to the immediate vicinity of the drive contacts (in the narrow parameter range where they appear at all).

Using the correlation signatures, we have mapped out the regime boundaries in AC transport (shown in fig. 1d) for the graphene device depicted in fig. 1a. The hydrodynamic regime emerges for lmr≳ L and lmc≲ 0.3L, with the region lmc > 0.3 L occupied by the ballistic regime. Transport for all lmr≳ L, i.e., both the hydrodynamic and ballistic regimes, occurs through vortex dynamics, a reflection of the high degree of degeneracy between the two regimes. However, the dynamics have different choreographies: For the device shown, the vortices in the hydrodynamic regime form at the drive contacts, whereas they form at the top and bottom boundaries in the ballistic regime. We emphasize here that the voltage-voltage correlation being used to discern the two regimes probes the nature of signal propagation within the device and is largely independent of the dynamics of vortices.

A surprising outcome of our calculations is the presence of current vortices in the ballistic regime; structures which are typically associated with a fluid. We trace this back to the fact that both non-Ohmic regimes conserve momentum and present an argument using Landau theory that vortices arise directly as a consequence of enstrophy conservation, which in turn results from momentum conservation in two-dimensions.

Our results have been obtained using bolt <cit.>, a package we have developed to directly solve for quasiparticle transport at the kinetic level with band structures and collision operators as inputs. The kinetic framework is crucial to access the ballistic regime, which is not amenable to effective fluid models. The package uses a high-resolution numerical scheme that deterministically discretizes the four-dimensional phase space of 2D systems. While the computational complexity involved is massive, it is nevertheless within reach of present day GPUs. The package exploits these with sufficient efficiency so as to allow time-dependent device simulations to be routinely carried out in a matter of hours, with realistic device dimensions and interaction time scales.

Finally, while our calculations were performed for a Fermi liquid, our phase diagnostics in AC transport may also hold for strongly interacting quantum systems that admit a hydrodynamic description <cit.>, such as the Dirac fluid which arises at charge neutrality in graphene <cit.>. Specifically, the current-voltage phase correlation is sensitive to any nonlocality in the current-voltage relation, regardless of its precise form. So it will immediately detect the presence of a hydrodynamic regime, whether this originates from a collection of rapidly interacting quasiparticles or in a strongly interacting system devoid of quasiparticles. Similarly, the voltage-voltage phase correlation can be used to ensure the dominance of interactions and rule out ballistic effects.



§ LONG-RANGE COULOMB INTERACTIONS
 
In writing the governing equation (<ref>), we have omitted the term -e𝐄·∂ f/∂𝐩 on the left hand side. This term incorporates backreaction of the self-consistent long-range Coulomb interactions, and couples the Boltzmann equation (<ref>) to the 3D Poisson equation (<ref>). Here we show that the dynamical effects of this term are accounted for at linear order through a renormalized chemical potential.

We begin by assuming that the system is driven at constant temperature using a small background chemical potential gradient ∇μ≡ (∂μ/∂ x, ∂μ/∂ y) ∼ϵ (μF/L), where μF is the Fermi level, L is the device scale and ϵ≪ 1 is a dimensionless bookkeeping parameter. This gradient results in a distribution f = f_0(ℰ, μ, T) + δ f, where f_0 is the background local equilibrium distribution and δ f ∼𝒪(ϵ) is a small perturbation. We then have the force due to a self-consistent electric field in the plane of the graphene strip, -e𝐄(𝐱) = e∇ V(𝐱, z=z_0) ∼𝒪(ϵ), where z=z_0 is the vertical location of the strip. Substituting f in (<ref>), now with -e𝐄·∂ f/∂𝐩 included, the LHS is

    ∂δ f/∂ t + 𝐯·∂δ f/∂𝐱 - (e𝐄 + ∇μ)·∂ f_0/∂𝐩 - e𝐄·∂δ f/∂𝐩

where we have used 𝐯·∂ f_0/∂𝐱 = - ∇μ· (∂ f_0/∂𝐩); this transformation makes it explicit that ∇μ is a long-range force.

As seen from (<ref>), in the presence of Coulomb interactions, the effective long-range force to 𝒪(ϵ) is -∇μ̅≡ -∇(μ - eV). Further, (<ref>) and (<ref>) are the same at linear order. The last term in (<ref>) which is not included in (<ref>) is only 𝒪(ϵ^2) and is therefore ordered out.

The equivalence at linear order between (<ref>) and (<ref>) allows for a convenient computational scheme in which we can solve the Boltzmann equation in a manner decoupled from the Poisson equation, while still taking into account the effect of long-range Coulomb interactions. We simply begin with an equilibrium f_0 set to a desired doping and setup transport by imposing a non-equilibrium f at the location of the driving contacts. The price to pay however, is that we can only measure μ̅ from the simulation and not V (or μ).

A current-voltage experiment in fact measures V. To obtain this, we still need to solve the 3D Poisson equation with the appropriate boundary conditions and sourced by the 2D carrier density n(𝐱) = ⟨ f ⟩, where f is computed from (<ref>) (or equivalently (<ref>)). Notably, this computation is a post-processing step and is not required for the evolution. As described in <ref>, we use an approximate solution to the 3D Poisson equation (<ref>) for the FET geometry that is widely used in experiments to obtain V.



 §.§ Drift vs Diffusion currents
 

We now proceed to check the relative contribution of each term in the total long-range force -∇μ̅. The term -∇μ gives rise to a diffusion current whereas e∇ V generates a drift current. Using (<ref>) to express ∇ V in terms of ∇ n and writing ∇μ = (∂μ/∂ n) ∇ n, we have,

    ∇μ̅    = ∇μ - e ∇ V 
        = (e^2/CQ + e^2/C) ∇ n

where CQ = e^2 ∂ n/∂μ is the quantum capacitance <cit.>. We evaluate CQ for graphene using μ = ħ vF√(π n) at an equilibrium carrier density n_0 corresponding to the Fermi energy μF. Using C = ϵr/(4π d) for a single gated device with a dielectric substrate ϵr and thickness d, the ratio of the drift to the diffusion term is,

    8·αee· (kFd) ≫ 1

where αee = e^2/(ϵrħ vF) ∼𝒪(1) is a dimensionless number called the graphene fine-structure constant and kF = μF/(ħ vF) is the Fermi wavenumber. With a typical d∼ 100 nm, we have kFd ≫ 1. From (<ref>), we see that the drift current dominates.



 §.§ Two different approximations to compute V

With the diffusion current being negligible, a simple way to obtain voltages is to approximate -e∇ V ≈∇μ̅, thus avoiding the Poisson equation altogether. Recall that μ̅ (written as μ in the main text) is obtained at every time step as a solution of the constraints (<ref>, <ref>, <ref>) that match f to a local equilibrium f_0. The voltage thus obtained is indeed consistent with that obtained using the approximate solution to the Poisson equation (<ref>), since infact (<ref>) has been originally used in <ref> to arrive at the conclusion that the diffusion current is negligible.



§ INCOMPRESSIBLE FLUID MODELS
 
Several calculations of current-voltage relations based on fluid models (e.g., <cit.>) assume that the electron fluid is incompressible, and therefore δ n = 0. However, in a kinetic calculation, density perturbations are inevitable since they appear at linear order, δ n = ⟨δ f ⟩∼𝒪(ϵ). Crucially, in this paper, these compressible fluctuations δ n are used to compute V using an approximate solution (<ref>) to the 3D Poisson equation. Two questions arise: (1) How are voltages calculated when incompressibility is assumed?, and (2) Are they consistent with the approach to compute voltages followed in this paper which explicitly tracks compressible fluctuations?

Consider linearized fluid models described by the charge (<ref>) and momentum (<ref>) conservation equations,

    ∂ n/∂ t    = - n_0∇·𝐯d
    ∂𝐯d/∂ t   = - 1/mn_0∇ P - e/m𝐄 + ν∇^2𝐯d - 𝐯d/τmr

where n_0 is the background charge density, P is the pressure, m is an effective mass, and ν is the kinematic viscosity of the electron fluid. The last term in (<ref>) is a momentum loss term parametrized by a phenomenological timescale τmr. The model requires a closure relation for the pressure P, which to a good approximation in the μ≫ kT limit is P ≃ nμ_F. Equations (<ref>) and (<ref>) are adopted from <cit.>, but with the time-derivatives and the pressure term included here.

The model thus has three unknowns, (n, 𝐯d, V), and therefore requires three equations; the fluid equations (<ref>, <ref>) along with the 3D Poisson equation (<ref>). However, the form of the force in (<ref>) allows for the unknowns to be reduced by absorbing P/(mn_0) and -(e/m) V into a single scalar Φ = P/(mn_0) - (e/m)V. The fluid system is then forced by the gradient of this effective scalar, analogous to μ̅ in the linearized kinetic equation (<ref>). 

Rewriting (<ref>) and (<ref>) in terms of Φ and taking the steady-state limit (∂/∂ t → 0),

    ∇·𝐯d    = 0 
    
    - ∇Φ + ν∇^2𝐯d    = 𝐯d/τmr

Note that the steady-state limit removes the dependence of n in (<ref>). There are now only two variables (𝐯d, Φ), for which just the equations (<ref>, <ref>) are sufficient. This is similar to the kinetic scheme discussed in appendix <ref> where the 3D Poisson equation is not necessary to compute f to linear order.

However, just as described in <ref>, the reduced system only yields Φ and not V. At this point, one can invoke incompressibility (n = n_0 ≡ constant) to set ∇ P = 0. This immediately reduces the unknowns to (𝐯d, V), thus fully determining the current-voltage characteristics. It is then not necessary to solve the 3D Poisson equation (<ref>), which only determines the background carrier density n_0.

The incompressible approximation is equivalent to stating that the diffusion current that arises due to compressible fluctuations is zero. This is consistent with <ref> where it is seen that although both diffusion and drift currents appear at linear order in a kinetic calculation, the ratio of the prefactors (<ref>) are such that diffusion is negligible. Further, this ratio has been computed using (<ref>), which we use to obtain voltages in the main text. The voltages thus obtained are therefore consistent with those obtained by fluid models in the incompressible limit (whenever the fluid models are applicable: lmc,lmr≪ L).



We thank Shyam Sankaran for design and performance optimizations in bolt that allowed for an efficient parameter exploration.
We thank Krishnendu Chatterjee and Joshua Matthew for insightful discussions.
MC, GK and DS thank Quazar Technologies for funding this research effort.
RS acknowledges start-up funding from the Department of Materials Science and Engineering and
computational resources at the Center for Computational Innovations at Rensselaer Polytechnic Institute.


5

[Ashcroft & Mermin(1976)]AshMer N. Ashcroft, & D. Mermin, Solid State Physics, Brooks/Cole, 1976

[Mayorov et al.(2011)]GrapheneMicronBallistic A. S. Mayorov, R. V. Gorbachev, S. V. Morozov, et al., Nano Letters 2011 11 (6), 2396-2399

[de Jong & Molenkamp(1995)]deJong M. J. M. de Jong, & L. W. Molenkamp, Phys. Rev. B 51, 13389 (1995)

[Molenkamp & de Jong(1994)]Molenkamp L. W. Molenkamp, & M. J. M. de Jong, Phys. Rev. B 49, 5038 (1994)

[Braem et al.(2018)]ScanningGateGaAsB. A. Braem, F. M. D. Pellegrino, A. Principi, et al., Phys. Rev. B 98, 241304(R) (2018)

[Gusev et al.(2018)]ViscousFlow2DGas G. M. Gusev, A. D. Levin, E. V. Levinson, and A. K. Bakarov, AIP Advances, 8, 025318 (2018)

[Moll et al.(2016)]PdCoO2 P. J. W. Moll, P. Kushwaha, N. Nandi, B. Schmidt, and A. P. Mackenzie, Science 351, 1061-1064 (2016)

[Gooth et al.(2018)]WP2 J. Gooth, F. Menges, N. Kumar, et al., Nature Comm. 9, 4093 (2018)

[Bandurin et al.(2016)]NegativeLocal D. A. Bandurin, I. Torre, R. Krishna Kumar, et al., Science 351, 1055-1058 (2016)

[Krishna Kumar et al.(2017)]SuperBallisticExpt R. Krishna Kumar, D. A. Bandurin, F. M. D. Pellegrino, et al., Nature Physics 13, 1182 (2017)

[Bandurin et al.(2018)]FluidityOnsetExpt D. A. Bandurin, Shytov, A., Levitov, L., et al., Nature Comm. 9, 4533 (2018)

[Berdyugin et al.(2019)]GrapheneMHD A. I. Berdyugin et al., Science 10.1126/science.aau0685 (2019)

[Gurzhi(1968)]Gurzhi R. N. Gurzhi, Soviet Physics Uspekhi, 11, 255 (1968)

[Guo et al.(2017)]SuperBallisticTheory H. Guo, E. Ilseven, G. Falkovich, & L. Levitov, Proceedings of the National Academy of Science, 114, 3068 (2017)

[Guo et al.(2017)]StokesParadox1 H. Guo, E. Ilseven, G. Falkovich, & L. Levitov, arxiv:1612.09239

[Lucas(2017)]StokesParadox2 A. Lucas, Phys. Rev. B 95, 115425 (2017)

[Shytov et al.(2018)]FluidityOnsetTheory Shytov, A., Kong, J. F., Falkovich, G., and Levitov, L., Phys. Rev. Lett. 121, 176805 (2018)

[Alekseev(2016)]NegativeMagnetoresistance P. S. Alekseev, Phys. Rev. Lett. 117, 166601 (2016)

[Principi & Vignale(2015)]ViolationWF A. Principi, & G. Vignale, Phys. Rev. Lett. 115, 056603 (2015)

[Scaffidi et al.(2017)]HydroHallViscosity1 T. Scaffidi, N. Nandi, B. Schmidt, A. P. Mackenzie, & J. E. Moore, Phys. Rev. Lett. 118, 226601 (2017)

[Delacrètaz & Gromov(2017)]HydroHallViscosity2 L. V. Delacrètaz, & A. Gromov, Phys. Rev. Lett. 119, 226602 (2017)

[Pellegrino et al.(2017)]HydroHallViscosity3 F. M. D. Pellegrino, I. Torre, & M. Polini, Phys. Rev. B 96, 195401 (2017)


[Chow et al.(1996)]PhononEmission E. Chow, H. P. Wei, S. M. Girvin, & M. Shayegan, Phys. Rev. Lett. 77, 1143 (1996)

[Govorov & Heremans(2004)]FermiJets A. O. Govorov, & J. J. Heremans, Phys. Rev. Lett. 92, 026803 (2004)


[Tomadin et al.(2013)]Tomadin2013 A. Tomadin, & M. Polini, Phys. Rev. B 88, 205426 (2013)

[Tomadin et al.(2014)]CorbinoDisc A. Tomadin, G. Vignale, & M. Polini, Phys. Rev. Lett. 113, 235901 (2014)

[Lucas(2016)]SoundModes A. Lucas, Phys. Rev. B 93, 245153 (2016)

[Semenyakin et al.(2018)]FalkovichAC M. Semenyakin & G. Falkovich, Phys. Rev. B 97, 085127 (2018)

[Levitov & Falkovich(2016)]LF2016 L. Levitov, & G. Falkovich, Nature Physics 12, 672 (2016)

[Torre et al.(2015)]Torre2015 I. Torre, A. Tomadin, A. K. Geim, & M. Polini, Phys. Rev. B 92, 165433 (2015)

[Pellegrino et al.(2016)]WhirlpoolsOrNot F. M. D. Pellegrino, I. Torre, A. K. Geim, & M. Polini, Phys. Rev. B 94, 155414 (2016)

[Falkovich & Levitov(2017)]LF2017 G. Falkovich, & L. Levitov, Phys. Rev. Lett. 119, 066601 (2017)

[Dyakonov & Shur(1993)]DyaShu M. Dyakonov and M. Shur, Phys. Rev. Lett. 71, 2465 (1993)

[Svintsov et al.(2012)]HydroElectronHolePlasma D. Svintsov, V. Vyurkov, S. Yurchenko, T. Otsuji and V. Ryzhii, Journal of Applied Physics 111, 083715 (2012)

[Svintsov(2018)]HydroToBallisticCrossoverDirac D. Svintsov, Phys. Rev. B 97, 121405(R) (2018)

[Gabbana et al.(2018)]PreTurbulenceProspects A. Gabbana, M. Polini, S. Succi, R. Tripiccione, and F. M. D. Pellegrino, Phys. Rev. Lett. 121, 236602 (2018)

[Chandra et al.(in prep)]bolt M. Chandra, S. Sankaran, P. Yalamanchili (in prep), <www.quazartech.com/bolt>

[Kraichnan(1975)]Enstrophy R. H. Kraichnan,  1975, Journal of Fluid Mechanics, 67, 155 

[Lozovik et al.(2015)]GrapheneCq Yu. E. Lozovik, A. A. Sokolik, and A. D. Zabolotskiy, Phys. Rev. B 91, 075416 (2015)

[Andreev et al.(2011)]HydroStronglyCorrelated A. V. Andreev, S. A. Kivelson, and B. Spivak, Phys. Rev. Lett. 106, 256804 (2011)

[Crossno et al.(2016)]DiracFluid J. Crossno, J. K. Shi, K. Wang, et al., Science 351, 1058-1061 (2016)

[Gallagher et al.(2019)]GrapheneQuantumCriticalConductivity P. Gallagher et al., Science 10.1126/science.aat8687 (2019)

[Lucas et al.(2016)]DiracFluidTheory A. Lucas, J. Crossno, K. C. Fong, P. Kim, & S. Sachdev, Phys. Rev. B 93, 075426 (2016)

[Müller et al.(2009)]GraphenePerfectFluid M. Müller, J. Schmalian, & L. Fritz, Phys. Rev. Lett. 103, 025301 (2009)



