Quantum oscillations from a pair-density wave

A pair-density wave state has been suggested to exist in underdoped cuprate superconductors, with some supporting experimental evidence emerging over the past few years from scanning tunneling spectroscopy. Several studies have also linked the observed quantum oscillations in these systems to a reconstruction of the Fermi surface by a pair-density wave. Here, we show, using semiclassical analysis and numerical calculations, that a Fermi pocket created by first-order scattering from a pair-density wave cannot induce such oscillations. In contrast, pockets resulting from second-order scattering can cause oscillations. We consider the effects of a finite pair-density wave correlation length on the signal, and demonstrate that it is only weakly sensitive to disorder in the form of $\pi$-phase slips. Finally, we discuss our results in the context of the cuprates and show that a bidirectional pair-density wave may produce observed oscillation frequencies.


I. INTRODUCTION
The underdoped cuprate high-temperature superconductors are known to harbor a variety of electronic orders. 1 Among them, charge-density waves (CDW) have attracted attention over recent years owing to a series of experimental observations. 2 Much more illusive is the pair-density wave (PDW) that is associated with a spatially oscillating superconducting order parameter of zero mean. 3 Such a state has its origin in the superconducting FFLO phase, 4,5 which can emerge in a magnetic field. Subsequently, the PDW was conjectured to occur without explicit time-reversal breaking in the cuprates. [6][7][8] Numerical studies suggest that the PDW may be energetically close to the uniform superconducting state, 6,[9][10][11][12] and recent scanning tunneling spectroscopy gives evidence that it is realized within halos surrounding Abrikosov vortex cores in Bi 2 Sr 2 CaCu 2 O 8+x . 13 Electronic orders that break translational symmetry naturally lead to redistribution of spectral weight in momentum space. If the latter contains gapless Fermi segments, as is the case in the underdoped cuprates, this can result in the formation of Fermi pockets, which give rise to quantum oscillations. It is commonly believed that the observed quantum oscillations in the high-temperature superconductors are due to such a scenario. 14 In particular, the established presence of a bidirectional CDW in these materials offers a natural route for formation of electron-like pockets with the required area to match the oscillations' frequency. 15,16 Nevertheless, the fact that the correlation length of the bidirectional CDW is shorter than the cyclotron radius of the closed orbit responsible for the oscillations casts doubt on this picture. 17 Other options, including unidirectional CDW, [18][19][20] and coexisting d-wave superconductivity with d-density wave, 21,22 or loop current order, 23,24 were also considered.
A PDW is appealing from the perspective of generating quantum oscillations in that it provides a superconducting state with gapless Fermi arcs. 25 This fact has led to suggestions that a PDW is the source of the oscillations in the pseudogap regime. [26][27][28] However, these studies differ in their reconstruction schemes. Refs. 26 and 27 considered pockets, which we dub first-order pockets, that are generated by first-order scattering and comprised of Fermi segments shifted by the PDW wave vector Q x (or Q y in the case of bidirectional order). Yet, their ability to produce oscillations was subsequently questioned by Ref. 28 on the basis that they include both electron-like and hole-like pieces. Instead, other pockets, designated by us as being second order, were suggested. These are constructed from 2Q x and 2Q y translated Fermi segments, all having the same character.
Motivated by the above disagreement, the lack of detailed and approximation-free calculations, and by the experimental situation we carry out a theoretical investigation of quantum oscillations from a PDW. We show, using both semiclassical analysis and numerical calculations, that the coupling of the Bogoliubov quasiparticles to the superfluid velocity field around vortices, v s (r), plays a pivotal role in establishing the structure of the local density of states (DOS) and hence of the oscillations spectrum. If one ignores this coupling then every Fermi pocket, regardless of its order, supports a periodic semiclassical orbit along constant energy contours E(k), where k is the gauge invariant crystal momentum. However, in the presence of the coupling semiclassical motion takes place along constant E(k ± mv s / ) electron-like and hole-like segments, respectively. Hence, first-order pockets that contain both types of segments generally do not sustain closed momentum-space orbits in the presence of a magnetic field. Moreover, in real space the motion of electrons and holes occurs under the influence of opposite effective magnetic fields, further hindering the formation of periodic orbits. In contrast, a semiclassical wave packet maintains its character around a secondorder pocket and is able to execute periodic motion that is largely free of these problems. The outcome is similar to the one arising from a 2Q x,y CDW, except for the electron-hole mixing that is present near PDW scattering points. This mixing leads to broadened oscillations at a frequency that nearly obeys the familiar Onsager relation to the area of the pocket.
The effect of PDW phase disorder can be analysed semiclassically along the lines of Ref. 17. It takes the form of a Berry phase, to which every Bragg scattering event off the PDW contributes the local PDW phase. The result is a generic suppression of the oscillatory signal by a Dingle factor that decays exponentially with the inverse PDW correlation length. However, we note an interesting possibility that arises from the fact that the oscillatory DOS originates from second-order pockets, where each scattering between Fermi segments contributes twice the local PDW phase. Hence, disordering of the PDW by π-phase slips leads to a Berry phase, which is a multiple of 2π and has no effect on the quantization. We show that while such complete immunity is not exact the DOS is indeed much less sensitive to disorder of this type.
The next three sections are devoted to substantiating the above results by defining a model Hamiltonian, analyzing it semiclassically and diagonalizing it numerically. We then discuss the relevance of our findings to the cuprates. We argue that it is possible that the quantum oscillations in YBa 2 Cu 3 O 6+x and HgBa 2 CuO 4+δ are due to a yet undiscovered PDW that is bidirectional with half the wave vector of the CDW that was detected in these systems at low magnetic fields. We also calculate the expected oscillation frequency from a PDW of the nature suggested by the observations in Bi 2 Sr 2 CaCu 2 O 8+x .

II. MODEL
We consider a model of a bidirectional d-wave PDW residing on the bonds of a square lattice with unit lattice constant and size L × L. The PDW wave vectors in the α = x, y directions are taken to be commensurate Q α = 2π/λ α , with integer λ α . We allow for independent positional disorder in the x and y PDW components (via phases θ x and θ y ), but assume that they share the same superconducting phase, φ, which is disordered by common vortices. We neglect PDW amplitude fluctuations due to vortices or otherwise. Incorporating the effects of a transverse magnetic field, B = Bẑ, by a Peierls substitution we are led to study the following Bogoliubov-de Gennes (BdG) Hamiltonian where Ψ † r = (c † r↑ , c r↓ ), τ are the Pauli matrices, and Here, A r,a = (e/ c)A r+a/2 ·a, and the d-wave form factor is given by ∆ α ±x = ∆ α , ∆ α ±ŷ = −∆ α . For the phase field we use the decomposition φ r+a/2 = φ r + ∇φ r+a/2 · a/2 ≡ φ r + ∇φ r,a /2, where the site field φ r is defined mod 2π but the bond gradient field ∇φ r+a/2 is single valued. As a result the phase factor in Eq. (3) is well defined and one also finds that φ r+a = φ r + ∇φ r,a mod 2π. A similar decomposition is used for θ α . Next, we remove the site variables φ r from H by the single valued unitary transformation, 29 and introduce the superfluid velocity field v s r+a/2 = ( /2m)[∇φ r+a/2 + (2e/ c)A r+a/2 ]. Consequently, in terms of v s r,a = (m/ )v s r+a/2 · a, the Hamiltonian becomes where h r,a (v s ) and ∆ r,a (0) are defined by the appropriate substitutions into Eqs. (2) and (3).

A. Solution for constant fields
Anticipating the semiclassical treatment of the next section, we are interested in the Bloch states of the above Hamiltonian when the "external fields" take constant values, i.e., A r+a/2 = A, v s r+a/2 = v s , and θ α r+a/2 = θ α . These states serve us in order to construct a wave packet whose motion we analyze in the case where the fields vary slowly in space. Given the translational symmetry of the problem at hand we expand where Ψ † k = (c † k↑ , c −k↓ ). The sum over the gauge invariant crystal momentum k = (2π/L)(n xx + n yŷ ), runs over the reduced Brillouin zone (RBZ), n α = −(N α − 1)/2, · · · , (N α − 1)/2, with N α = L/λ α the number of unit cells in the α direction. The Λ = λ x λ y Bragg vectors are q = m x Q x +m y Q y with m α = −(λ α −1)/2, · · · , (λ α − 1)/2. Expressed in momentum space the Hamiltonian ξ k = −2t(cos k x + cos k y ) − µ, and ∆ α k = ∆ α (cos k x − cos k y ).
The spectrum of H is discussed in Appendix A. There we show that in the presence of v s it consists of doubly degenerate bands whose number, Λ k , may change over the RBZ but sums up to the total number of degrees of freedom, owing to Λ −k = Λ − Λ k . The excitations are created from the ground state, |g , by the quasiparticle operators γ † kn± , where n is the band index, and take the form The periodic parts, ϕ kn± , of the Bloch states are given in Eq. (A8) and depend on v s and θ α .

A. Wave packet
Our goal is to obtain the DOS by semiclassical quantization of periodic cyclotron orbits. To this end, we follow Ref. 30 and consider a quasiparticle wave packet that is centered around k c in momentum space and constructed from Bloch states, Eq. (8), of the "local Hamiltonian" with η = ±. Here, W (k) is a real weighting function that is peaked around k = 0, with width much smaller than Q x,y and the typical wave vector of variations in the external fields. It is also assumed to be periodic over the RBZ and square-normalized k W 2 (k) = 1. The k-space Berry connection is given by with its r c dependence inherited from the dependence of ϕ on A,v s and θ α , and where the overlap is defined by a sum over a unit cell in real space We note that Eq. (A8) implies ϕ kn− = iτ 2 ϕ * kn+ , which can be used to show A kn+ = A kn− . Henceforth, we assume that the motion takes place within a single band and suppress the band indices.
Our choice for the phase of the wave packet in Eq. (9) is tacitly related to the identity of its position in real space. While there is no ambiguity with regards to this question in the case of normal electronic systems the issue is less clear for superconductors. 31,32 Most generally, one may consider a position operator of the formr = r rΨ † r RΨ r , where R is a hermitian matrix, and require that : Ψ|r|Ψ : ≡ Ψ|r|Ψ − g|r|g = r c . For R = τ 3 the operator is related to the center of mass (or charge), while R = I corresponds to the spin center. Because the charge of a superconducting quasiparticle may change along its travel in momentum space, whereas its spin is conserved (barring spin-orbit coupling), it was suggested that the spin center should be identified with the center of the wave packet. 32 Here, we would like to offer another argument in favor of using R = I. In order to achieve self-consistency of the semiclassical treatment, the wave packet (9) needs to be concentrated in real space such that the local Hamiltonian constitutes an approximate generator of its dynamics. Ideally, this would mean that for H = r H r [F (r)], where F (r) represents the collection of external fields, : Ψ|H|Ψ : = : Ψ| r H r [F (r c )]|Ψ :, at least when H r is expanded to first order in r − r c . We note from Eq. (5) that when θ α is constant ∂H r /∂r c ∝ a Ψ † r+a Ψ r . Thus, the above requirement is fulfilled in the absence of PDW phase disorder if : Ψ|r|Ψ : = r c , for R = I. We demonstrate in Appendix B that this expectation value holds for wave packet (9).

B. Dynamics
The semiclassical dynamics of the coordinates k c and r c is derived from the Lagrangian L =: and the second Berry connection is also independent of η.

C. First-order scattering
Quantum oscillations originate from the alternating presence and absence of states at the chemical potential, as the magnetic field is varied. For a weak PDW, which is our focus here, such states correspond to periodic solutions of Eqs. (15) and (16), for which the motion is largely along sections of the unperturbed Fermi surface. On these sections the eigenstates are either "electron-like" with character close to c † ηk,η↑ |g , and energy ξ(k c + mv s / ), or "hole-like" given approximately by c −ηk,−η↑ |g with energy −ξ(k c − mv s / ), see Eq. (7). The PDW may cause scattering of the wave packet between Fermi segments at points that are connected by an integer combination m x Q x + m y Q y , where m x + m y is the order at which the process appears in perturbation theory. Odd-order scattering connects segments of opposite character, while even-order scattering preserves the nature of the segment. Refs. 26 and 27 claimed to find periodic orbits that emerge from first-order scattering. We proceed to show that this claim is erroneous and is caused by neglecting the effects of the superfluid velocity on the motion of quasiparticles.
To this end, we consider the semiclassical dynamics away from the scattering points in an ordered system (θ α = 0), where the Berry curvatures vanish. On an electron-like segment the energy is a function, E(k e ), of the variable k e = k c + mv s / . Consequently, Eqs. (15) and (16) may be cast into the forṁ with the effective magnetic field B eη =B+η(mc/e)∇ rc × v s . For v s that is generated by a collection of vortices at positions R n one finds B eη = ηB + (1 − η)(hc/2eB) n δ(r c − R n )B. Note that the equations imply a reversed r c motion for down-spin excitations (η = −) as compared to up-spin excitations (η = +). This is consistent with the fact that for down electrons r c is the inverted center of mass position. Conversely, on a hole-like segment the energy is a function of the variable k h = k c − mv s / and the equations of motion becomė If one follows Ref. 26 and sets v s = 0 then k e = k h = k c and B eη = B hη = B. As a result, Eqs. (18) and (20) imply that motion in k c -space takes place on constant ξ(k c ) contours, or more generally, according to Eq. (16) along constant E(k c ) contours. The equations also indicate that the sense of rotation is changed across each As a result a (dark grey) pocket appears, constructed from shifted Fermi sections. Motion in rc-space proceeds along a rotated and scaled version of the kc orbit. In contrast, when accounting for the effects of vs the sense of rotation is the same on all segments and one may conclude that a (light grey) pocket emerges due to 1 → 2 orbit. However, the picture is misleading since the two segments are drawn as function of two different parameters: kc ± mvs/ . Moreover, such a pocket would produce an open skipping trajectory for rc.
Bragg scattering. Assuming that B is small enough to disregard magnetic breakdown this condition determines the shape of the consequent pockets, see Fig. 1. The latter coincide with the pockets that emerge from diagonalizing H for B = v s = 0. Concomitantly, Eqs. (19) and (21) tell us that r c also executes periodic motion that is derived from that of k c by a π/2 rotation and scaling by l 2 B = c/eB. Upon semiclassical quantization these periodic orbits would give rise to Landau levels and hence to oscillations. However, this is an artefact of the approximation v s = 0.
In the presence of v s , and as long as the semiclassical trajectory misses the vortex cores, B eη = −B hη . Thus, motion on electron-like and hole-like segments follows constant ξ(k e ) and ξ(k h ) contours, respectively, with the same sense of rotation for both cases. Superficially, this may lead to the conclusion that closed k-space pockets may still form, as suggested by the light gray pocket in Fig. 1. However, it should be noted that the constancy of the ξ contours is defined relative to two different variables: k c ± mv s (r c )/ . Therefore, generically, one does not expect that solving the coupled equations for k c and r c would yield a closed orbit. Even if a closed cycle is formed for k c , the fact that B eη = −B hη would lead by Eqs. (19) and (21) to an open skipping orbit for r c . Indeed, as we demonstrate by numerical calculations below, there is no evidence for quantum oscillation from first-order scattering off a PDW.

D. Second-order pockets
In light of the preceding discussion we are led to consider the possibility of orbits that maintain their character during semiclassical evolution. Under such circumstances we may expect the entire motion to unfold as function of a single momentum variable, k e or k h , and the real-space trajectory to evolve under a unique effective magnetic field. We will show that this expectation is partially borne out by the following analysis. Since each scattering off the PDW changes the identity of the orbit, it is necessary to construct it from Fermi segments connected by second (or in general even) order scattering. Here we will do so for a bidirectional PDW, see Fig. 2, although it can also be done for a unidirectional PDW.
For the remaining discussion, let us focus on an electron-like excitation. Away from the scattering points and for weak PDW the semiclassical dynamics is given to a good approximation by Eqs. (18) and (19). To understand the behavior near the scattering points we refer to Fig. 2 and examine the vicinity of the upper tip of the diamond pocket. There, the PDW mixes a state |1 = c † k,σ |F s on segment 1 with a state |2 = c † k+2Qx,σ |F s on segment 4 via an intermediary state |3 = c −k−Qx,−σ |F s . Here, σ = ± is the spin and |F s is the filled Fermi sea in the absence of the PDW. The reduced Hamiltonian within this subspace is where We are interested in the eigenstate of H eff that approaches |1 in region a away from the tip and |2 in region b, as defined in Fig. 2. The desired state is where and where the phase χ is to be chosen such that it approaches 2θ x for ∆ξ |P |, and 0 for ∆ξ −|P |. For concreteness, we pick More pertinent, from the perspective of elucidating the semiclassical dynamics, is the fact that the corresponding energy depends significantly on both k e and k h at the vicinity of the tip. Hence, while the k-space orbit essentially coincides with constant E(k e ) contours along the arcs of the diamond, this is no longer true near the scattering points. We note that the situation is different from the one encountered for reconstruction scenarios due to particle-hole orders, such as CDW. 17 There, at least in an ordered system, the k-space orbit is given precisely by the Fermi surface, thus leading to the Onsager relation. We will not analyse in detail the semiclassical motion near the tip, but instead assume that periodic solutions do exist, at least for a significant set of initial conditions. Our assumption is backed by numerical calculations, which we detail below, that demonstrate quantum oscillations, albeit with broadened peaks and a slightly shifted frequency compared to the v s = 0 case. We attribute this behavior to the spread in the semiclassical trajectories due to tip effects.
To proceed with semiclassical quantization of the periodic orbits we identify the conjugate momenta from the Lagrangian, Eq. (12), They enter the Bohr-Sommerfeld condition where n is integer and µ is the Maslov index of the trajectory. We first treat the ordered case θ α = 0, for which the Berry connections vanish owing to the fact that |ν is real, and Assuming that no vortex cores are encountered and using Eq. (19) to express dr c = l 2 B dk e ×ẑ, the quantization rule reads where A k is the area swept by the k e orbit and N v is the number of vortices encircled. Being a multiple of 2π the vortex contribution does not affect the quantization condition. The latter takes the familiar Onsager form, except that A k , as noted above, may deviate somewhat from the area of the pocket calculated for v s = 0.

E. PDW phase disorder
Before continuing with the analysis, let us note that the self-consistency argument of Section III A, which favors constructing the wave packet such that : Ψ|Ψ † r Ψ r |Ψ : = r c , fails when θ α is not constant since then ∂H r /∂r c ∝ a Ψ † r+a Ψ r . Nevertheless, we assume that in the limit of ∆ → 0, where the effects of the scattering events are concentrated at small regions of phase space, this choice is still optimal. Taking this point of view, we analyse the effects of PDW phase disorder within the same semiclassical framework used so far. We then contrast its predictions with numerical calculations of the model in the next section.
In the presence of PDW phase disorder the Berry connections do not vanish. With the help of Eq. (27) they evaluate near the upper tip to Their contribution to S/ , coming from the orbit section that crosses the tip from point (k i , r i ) on segment 4 to point (k f , r f ) on segment 1, is where [t i , t f ] is the time interval of the motion. The last equality is valid in the limit of vanishing PDW amplitude where (d/dt)(∆ξ/∆E) = 2δ(t 1 ), with t 1 the time at which the wave packet scatters at position r c = r 1 .
Collecting the contributions from the other scattering points yields −2η[θ x (r 1 )−θ y (r 2 )−θ x (r 3 )+θ y (r 4 )]. Since the amplitude of the fundamental harmonic of the DOS is proportional to exp(iS/ ), its average over the fluctuations of the PDW phases is suppressed by a Dingle factor. 17 In terms of the correlation lengths ξ x,y of θ x,y it takes the form where ∆ kx,y are the k-space distances between the leftright and bottom-up scattering points, respectively. We conclude this section by noting an interesting consequence of the fact that each scattering event adds ±2ηθ α to S/ . Counter to the general case, in the particular set of configurations where the PDW is disordered via π-phase slips this contribution is a multiple of 2π and therefore has no effect on the DOS. Because of the reservation expressed at the beginning of the discussion it is unclear to what extent this conclusion and Eq. (40) should be trusted. In the next section we provide numerical evidence, which show that although these results are not exact they hold to a good degree of approximation. While we do not have a microscopic reasoning for why disordering of the PDW should progress via creation of π-phase slips we note that it offers a concrete model, which supports quantum oscillations in a highly disordered system without affecting the frequency of the oscillations. This stands in contrast to first order scattering, e.g. from a CDW, where phase discommensurations do not suppress the signal but leads to its appearance at higher harmonics. 17 FIG. 3. DOS as function of inverse magnetic field for a unidirectional d-wave PDW. First-order scattering leads at B = 0 to the red pocket, as indicated by the lower inset. Turning off vs results in the DOS oscillations depicted in the upper inset, whose frequency is related by the Onsager rule to the area of the pocket. As shown by the main figure these oscillations are eliminated once one includes the coupling to the superfluid velocity due to a vortex liquid. New, faster oscillations appear at high fields, which are due to magnetic breakdown and correspond to an orbit around the original Fermi surface, depicted in blue in the lower inset.

A. Model and method
In order to establish the validity of the semiclassical results we use numerical methods to study the BdG Hamiltonian, Eq. (5). We do so for the cases of a bidirectional d-wave PDW with ∆ x = ∆ y = ∆, a unidirectional dwave PDW with ∆ x = ∆ and ∆ y = 0, and for a unidirectional s-wave PDW, for which the pairing term in Hamiltonian (1) reads ∆ r cos(Q x · r + θ r )Ψ † r Ψ r . We carry out the calculations on a rectangular square lattice with L x sites and open boundary conditions in the x-direction. We apply periodic boundary conditions along the y-direction, such that Q y L y = 2πn, where L y is the number of sites around the cylinder and n is integer. To conform with the boundary conditions we use the Landau gauge A = (0, Bx, 0). The superfluid velocity field is constructed by summing over N v = (2e/hc)BL x L y randomly placed vortex configurations, as described in Appendix C. Note that the allowed magnetic fields in our finite lattice are constrained by the fact that N v is integer.
Since quantum oscillations originate from the electronlike part of the DOS, 33 our goal is to calculate where Tr e is the partial trace over the electron-like part of Ψ. We concentrate on the B-dependence of the zeroenergy DOS, ρ(0), and present this quantity in the following, while using δ = 2.5×10 −4 . We have checked that the full DOS, defined by the trace over all the entries of Ψ, yields qualitatively similar results. In order to observe DOS oscillations one needs to handle rather large systems. Typically, we use systems of size L x ×L y = 960×72 and average over 2000-4000 vortex liquid realizations. For such sizes calculation of ρ(0) by straightforward diagonalization is tasking. Instead, we utilize the recursive Green's function method, 34 which allows to calculate the diagonal terms of the Green's function with computational cost that scales linearly with L x .

B. First-order scattering
To check the ability of first order pockets to induce DOS oscillations we have considered a number of models where such pockets arise at zero magnetic field. A representative example is presented in Fig. 3, corresponding to a period 6 unidirectional d-wave PDW with ∆/t = 0.4 and µ/t = −3.2. The lower inset depicts in red the pocket with area A k = 0.025A BZ , where A BZ is the area of the full Brillouin zone. The upper inset shows the DOS generated from Hamiltonian (5) upon substituting v s = 0. Sharp oscillations are visible, at a frequency that follows the Onsager relation F = (A k /A BZ )φ 0 /a 2 , where φ 0 = hc/e. They originate from closed orbits around the pocket as depicted in Fig. 1. We have found similar oscillations in all cases where first order pockets are present, as long as we set v s = 0. Equally generic, is the disappearance of the oscillations once v s is included in the Hamiltonian, as demonstrated by the main panel of Fig. 3. We have found no exception to this statement in any system with first order pockets. Note that Fig. 3 does show oscillations at high fields. However, they appear at a frequency F a 2 /φ 0 = 0.064, which is very nearly the ratio between the area of the original Fermi surface and A BZ . Hence, we associate the signal with a periodic orbit around the original Fermi surface due to magnetic breakdown. We have observed DOS oscillations from magnetic breakdown in many cases and found that they tend to commence at lower fields in the presence of vortices, as compared to the case v s = 0.

C. Second-order pockets
Similar to the situation with first-order pockets, second-order pockets also induce sharp DOS oscillations under the approximation v s = 0. An example is shown in the inset of Fig. 4, which corresponds to a system hosting a period 6 bidirectional d-wave PDW with ∆/t = 0.5 and µ/t = −0.966. The PDW generates a diamond-shaped pocket, as in Fig. 2, whose area is precisely related to the observed oscillation frequency via the Onsager rule.
However, in contrast to first-order pockets, the oscillations stemming from second order pockets survive the inclusion of v s , as illustrated by the main panel of Fig.  4. The DOS peaks are broader and slightly shifted compared to the signal for v s = 0. We attribute these effects to deviations of the semiclassical orbit from motion along constant energy contours near the tips of the pocket, as discussed in section III D. Nevertheless, the Onsager relation between the frequency of the DOS oscillations and the area of the pocket continues to hold to a good approximation. This fact is demonstrated by Fig. 5 where   FIG. 6. DOS oscillations for a system with a bidirectional d-wave PDW of the same parameters as in Fig. 4, but whose phase is smoothly disordered in space with a correlation length ξ. we plot F vs. A k for the model mentioned above, as well as for three other models. These include a period 8 unidirectional d-wave PDW with ∆/t = 0.3 and µ/t = −2.5 and two period 8 unidirectional s-wave PDWs, one (1) with ∆/t = 0.4 and µ/t = −3, and the second (2) with ∆/t = 0.25 and µ/t = −2.5.

D. PDW phase disorder
Next, we consider the effects of phase disorder. We start by calculating the DOS for the bidirectional d-wave PDW, whose parameters are given above, in the case where its phases vary smoothly in space along a single FIG. 8. DOS oscillations for a system with a bidirectional d-wave PDW of the same parameters as in Fig. 4, but whose phase is disordered via a random array of π-phase slips, resulting in a correlation length ξ. direction, θ α r+a/2 = θ(x). We generate the phases using a one-dimensional random walk with a step size that is normally distributed. The standard deviation of the step is chosen to be 2/ξ, such that ξ is the resulting correlation length for exp(iθ). In addition, we smooth every phase configuration by convoluting it with a Gaussian of width 2. Fig. 6 depicts the DOS for various values of ξ, where each trace is an average over 2000-4000 realizations that differ both by their vortex distributions and by their PDW phase fields. The decay of the oscillations with increasing disorder is evident. In order to quantify it we split the 1/B range into sections, each containing 4 oscillations. We then plot, in Fig. 7, the amplitude of the Fourier transform peak resulting from each section after normalizing it by the corresponding amplitude for the clean system. We only present peaks that are significantly discernible from the background. Based on our semiclassical result, Eq. (40), we expect an exponential decay of the signal with (ξB) −1 . Note, however, that owing to the one-dimensional nature of the disorder that we use, only scattering events in the x-direction involve a phase difference and contribute to the Dingle factor. Consequently, the latter becomes exp(−4l 2 B ∆k y /ξ). For the diamond pocket created by the PDW considered here, this leads, in units where φ 0 = 1, to exp(−0.584/ξB). This is about 10% slower than the best fit to the observed decay, depicted by the solid line in Fig. 7.
Finally, we consider phase disorder in the form of a one-dimensional array of π-phase slips. The latter is generated by a random walk, where a π jump is introduced at each step with probability (2ξ) −1 . The resulting θ α r+a/2 = θ(x) configuration is again smoothed out by convoluting it with a Gaussian of width 2. The results for the DOS are presented in Fig. 8 and show that the quantum oscillations are considerably more robust against disorder of this type, as expected on the basis of our semiclassical analysis.

V. DISCUSSION
As alluded to in the Introduction, circumstantial evidence for a PDW state exists only in underdoped Bi 2 Sr 2 CaCu 2 O 8+x (Bi2212). It takes the form of bidirectional charge modulations around Abrikosov vortex cores with approximate periods 4a and 8a, where a is the lattice constant. The appearance of two simultaneous CDWs with wave vectors Q and 2Q is a natural consequence of Q-PDW modulations existing in parallel to a uniform superconducting order. 28,35 In contrast, to date there are no measurements showing quantum oscillations from a reconstructed Fermi surface in this compound. The opposite situation holds for underdoped YBa 2 Cu 3 O 6+x (YBCO), YBa 2 Cu 4 O 8 and HgBa 2 CuO 4+δ (Hg1201), where quantum oscillations from small Fermi pockets have been detected, 3 but there is no evidence to a PDW (very recently, 36 high magnetic field-resilient superconductivity has been observed in YBCO, which may involve a PDW). It is then natural to ask whether the observed oscillations in the latter group of materials can be generated by a PDW, and what are the expected characteristics of PDW-induced quantum oscillations in Bi2212.
X-ray scattering has detected a bidirectional CDW in both YBCO and Hg1201. 2 There is a significant correlation, which follows the Onsager relation, between the area of the pocket obtained by folding the Fermi surface via the CDW wave vectors and the frequency of the quantum oscillations. 37 However, the correlation length of the bidirectional CDW in YBCO reaches only about 30 lattice constants, 38,39 while quantum oscillations commence at a magnetic field where the cyclotron radius is approximately three times longer. The Dingle factor associated with such a ratio between the two length scales would totally suppress the quantum oscillation signal. 17 This problem is even more pronounced in Hg1201. 37 In YBCO, a CDW with the same period and considerably longer-ranged correlations appears under strong magnetic fields. 38,39 Nevertheless, this CDW is unidirectional and it is not clear how it can give rise to pockets of the desired area. We have argued that the longer-range CDW likely nucleates around vortices, where superconductivity is strongly suppressed, and is oriented due to its Coulomb coupling to the chain layers. 40,41 Here, we would like to raise the possibility that the unidirectional CDW is actually a subsidiary order to a bidirectional PDW with twice the period, which competes with uniform superconductivity. Such a PDW is capable of producing nodal "diamond" pockets via second order scattering, as shown in Fig. 2, that match the observed oscillation frequencies. As mentioned, a PDW with wave vector Q can also combine with uniform superconductivity to produce a CDW at the same Q. To the best of our knowledge, such a signal has not been observed thus far in YBCO, 42 but it may be too weak to detect at high fields where the uniform order is effectively quenched. On the same note, one may contemplate the possibility that a bidirectional PDW is also the parent order of the bidirectional CDW. However, the fact that experiments show no signatures of an additional CDW with double the period under conditions where uniform superconductivity is still strong, 43 most likely rules it out.
Applying the above scenario to Hg1201 would lead to the expectation that strong and ordered bidirectional PDW should appear around vortex cores at higher fields than the highest field used so far (∼ 16 T) but below the onset field of quantum oscillations in this material (∼ 55 T). Such a PDW would likely be accompanied by a CDW of half the period (the same as of the observed low-field CDW), which is expected to be long-ranged and bidirectional, owing to the absence of orienting chains in the tetragonal Hg1201.
Finally, if a bidirectional PDW of period 8 exists in Bi2212, as the experiment suggests, 13 then it can lead to the formation of second-order pockets. In Fig. 9 we plot such pockets calculated for the case of a strong PDW ∆ = 0.5t, and using the tight-binding dispersion of Ref. 44. They are expected to induce DOS oscillations at a frequency of approximately 830T. We note, however, that if the PDW shares the very short correlation length of the observed CDWs within the halos then the oscillations would be massively damped. Furthermore, reducing the strength of the PDW in our calculations leads to small gaps that may be breached by magnetic breakdown and lead to different frequencies. ACKNOWLEDGMENTS We thank S. Kivelson for useful comments. This research was supported by the Israel Science Foundation (Grant No. 701/17).
We note that for v s = 0 the number, Λ k , of positive eigenvalues need not equal that of negative eigenvalues. The eigenvectors form the columns of U † The spectra of H(k) and H(−k) are related by symmetry. The blocks H qq (k), Eq. (7), obey the relation H * −q−q (−k) = −τ 2 H qq (k)τ 2 , from which it follows that where Given an eigenvector |ν −k,n of H(−k) corresponding to an eigenvalue E (positive or negative) we can apply Eq. (A2) to |ν −k,n * and find that T 2 |ν −k,n * is an eigenvector of H(k) with energy −E. Hence, we conclude that for non-degenerate E > 0 and up to a phase T 2 |ν −k,n * = |ν k,Λ k +n . Similarly, T 2 |ν −k,Λ −k +n * = |ν k,n , for nondegenerate E < 0. By proper orthogonalization these statements can be made true also for the degenerate case.

Appendix B: Expectation values of |Ψ
In order to calculate : Ψ|r|Ψ : for the position operator, r = r rΨ † r Ψ r , we use Eqs. (A4) and (A8) to find g|γ k η Ψ † r Ψ r γ † kη |g = g|Ψ † r Ψ r |g δ kk + η Λ L 2 e i(k−k )·r ϕ † k + (r)ϕ k+ (r), (B1) where η = ±. Combining this result with Eq. (9) gives : Ψ η |r|Ψ η : = iη k,k δ kk ∇ k W (k − k c )W (k − k c ) × e iη(k −k)·(rc−A k ) ϕ k + |ϕ k+ = r c , (B2) where we have used re ik·r = −i∇ k e ik·r and integrated by parts over k to arrive at the first equality. The final result then follows from the periodicity of W and the definition of A k . For the purpose of evaluating the semiclassical Lagrangian we are in need of : Ψ η |i d dt |Ψ η :. Here, we draw attention to the contribution originating from the implicit time dependence of γ † kη : Ψ ± | Λ 1/2 L r e i[k±( m vs− e c A)]·r × ± d dt e c A − m v s · r +ṙ c · ∇ rc Ψ † r ϕ k+ (r)|g : ϕ † k− (r)Ψ r |g : . (B3) Using similar steps to those leading to Eq. (B2) one finds that the first term evaluates to r c ·(d/dt)[(e/ c)A− (m/ )v s ], while the second gives rise to the Berry connection A r . Adding the remaining contributions and a total time derivative, which does not affect the equations of motion, results in Eq. (12).

Appendix C: Construction of the vortex liquid configuration
Constructing a vortex configuration on the cylinder requires a smooth phase-gradient field whose circulation vanishes around every plaquette, except the one containing the core. To this end, we start by considering a 2L x × 10 4 L y lattice with open boundary conditions, whose geometric center is defined as the origin. The phase on each site, φ r , takes the value of the site's azimuthal angle. To keep the phase gradient smooth it is defined by ∇φ r+a/2 = φ r+a −φ r −2πΘ(φ r+a −φ r −π)sgn(φ r+a −φ r ), (C1) where Θ is the step function. Next, we sum over 1000 replicas of this configuration, each with a vortex center that is successively translated by 10L y along the ydirection. The configuration obtained by closing the central L y strip on itself fulfills the above requirements, up to a very small flux on plaquettes crossing the seam. The vortex liquid is obtained by superposing N v such vortex configurations with randomly placed cores. It is used to calculate the v s field via v s r+a/2 = ( /2m)[∇φ r+a/2 + (2e/ c)A r+a/2 ]. Finally, to avoid any residual total current through the system we subtract the spatial average of the velocity from every y-bond.