Optical imprinting of superlattices in two-dimensional materials

Hwanmun Kim ,1,2,* Hossein Dehghani ,1,3,* Hideo Aoki ,4,5 Ivar Martin,6 and Mohammad Hafezi 1,2,3 1Joint Quantum Institute, NIST and University of Maryland, College Park, Maryland 20742, USA 2Department of Physics, University of Maryland, College Park, Maryland 20742, USA 3Departments of Electrical and Computer Engineering and Institute for Research in Electronics and Applied Physics, University of Maryland, College Park, Maryland 20742, USA 4Department of Physics, The University of Tokyo, Hongo, Tokyo 113-0033, Japan 5National Institute of Advanced Industrial Science and Technology (AIST), Tsukuba, Ibaraki 305-8568, Japan 6Materials Science Division, Argonne National Laboratory, Argonne, Illinois 60439, USA

In this paper, we propose a method to create superlattice structures in a 2D material by shining spatially periodic laser beams, as schematically shown in Fig. 1. We illustrate the idea with an example of monolayer graphene irradiated by a circularly polarized beam with a superlattice structure, where the beam amplitude is spatially periodic. To demonstrate the tunability of this superlattice structure and unique physics originating from the superlattice, we first study the case of a square superlattice and explore the topological phase transition induced by varying the superlattice size. Then, we investigate the topological phase transitions, when the square superlattice is sheared to a stretched hexagonal one. In particular, we examine the relationship between this topological phase transition and the role of lattice geometry in creating complex tunneling phases. Further, we demonstrate the possibility of creating more exotic lattices by superposing multiple lattices, with an example of tuning between a hexagonal and a kagome lattice where the flat bands can be obtained. These flat bands particularly can harbor strongly correlated phenomena in Floquet systems.

II. GRAPHENE WITH SPATIALLY PATTERNED LIGHT
Let us consider a monolayer graphene with the interatomic distance a and the tight-binding energy t between the nearest neighbors. The low-energy description for this monolayer FIG. 1. A 2D material irradiated by a spatially periodic CP light with frequency ω. Here, we use the example of a monolayer graphene. The superposition of multiple CP Gaussian beams generates a periodic amplitude pattern A 0 (r) with translation vectors L 1 and L 2 . |L 1 | = l. Upper inset: We denote the interatomic distance of the graphene as a and the tight-binding energy between the nearest neighbors as t. Lower inset: Each Gaussian beam has a peak amplitude A 0 and a half waist w (black lines). The overall beam amplitude (red line) results from the superposition of the Gaussian beams. graphene under the electromagnetic field A(r, t ) is given by where σ x , σ y , σ z are Pauli matrices acting on sublattice degrees of freedom, v = (3/2)ta is the Fermi velocity at Dirac points, and τ z = ±1 is the valley index [41]. In particular, if we shine the CP beam with spatial amplitude pattern A(r, t ) = A 0 (r)e iωt (x + iŷ) + c.c. (Fig. 1), the effective Floquet Hamiltonian to the first order in ω −1 becomes [31,[42][43][44][45][46][47] We denote the peak amplitude of A 0 (r) as A 0 . Then, Eq.
(2) becomes a valid description when frequency ω is high enough (ω evA 0 ) and the amplitude varies in length scale larger than a (A 0 /max{|∇A 0 (r)|} a). For brevity, we seth = 1 from here on.
We specifically study the superlattice structure created by a spatially periodic amplitude |A 0 (r)| = |A 0 (r + L 1 )| = |A 0 (r + L 2 )|. While the 2D material with spatially modulated beams has been studied in the different contexts [48][49][50], here we investigate the generation of a superlattice with spatially periodic beams. In particular, to make the beam experimentally relevant, we consider the superposition of CP Gaussian beams positioned on the superlattice, where w is the radius of each Gaussian beam. This beam configuration is achievable with recent progress in beamshaping technologies [20][21][22][23][24][25]. For the cases |L 1 |, |L 2 | = l a, the Brillouin-zone folding occurs on a momentum scale 1/l. Furthermore, the hybridization of Floquet sidebands is suppressed for v/l ω so that the low-energy description is captured by Eq.
(2) (see Appendix A). We obtain Bloch eigenstates |ψ m,k and eigenenergies E m,k , where m is the band index and k is the crystal momentum within the Brillouin zone set by reciprocal lattice vectors of L 1 and L 2 . Note that Eq. (2) preserves particle-hole symmetry (σ x H * eff σ x = −H eff ) and therefore the energy spectrum is symmetric with respect to the zero energy. Also, σ y H eff σ y = H eff | τ z →−τ z , so two valleys have the same spectrum and eigenstates up to a unitary operation, σ y . This also ensures that both valleys have the same Chern number. For brevity, let us only consider the τ z = 1 valley from now on.

III. ILLUMINATION OF SQUARE SUPERLATTICE
We first consider the simplest case of a square superlattice, L 1 = lx and L 2 = lŷ. Before directly diagonalizing Eq. (2), we can make some speculations. First of all, the contribution from the spatial average of |A 0 (r)| opens up the gap around the zero energy ( b ) as in the case of the graphene under the CP uniform light, where the Chern number, C 1 , of the first band above E = 0 is nonzero [31][32][33]48,51]. C 1 remains nonzero for small l, as far as the maximum kinetic energy within the Brillouin zone, which is of the order of v/l, is much larger than the spatial Fourier components of the σ z term in Eq. (2), which is of the order of e 2 v 2 A 2 0 /ω. On the other hand, as l → ∞, the contribution of the kinetic term becomes negligible and therefore the bands become flat. Also, the Bloch wave functions look similar regardless of k and therefore the bands become topologically trivial. Therefore, there must be a topological phase transition where C 1 changes from a nonzero value to zero as we increase l. This topological transition would occur at a superlattice size that makes the two energy scales e 2 v 2 A 2 0 /ω and v/l comparable to each other. For a succinct description of this phase transition, we use the rescaled superlattice size so that the critical superlattice size χ c is O(1). Here, χ represents the ratio of the effective superlattice potential over the kinetic energy.
To study the detail of this topological phase transition, we numerically diagonalize Eq. (2) as shown in Fig. 2(a). Along with the energy spectrum, we present the Chern number C of each band calculated based on Ref. [52]. In Fig. 2, we set A 0 = 0.006(ea) −1 , ω = 0.06t, and w/l = 0.3. With these parameters, we can check that the topological phase transition occurs at χ c = 0.965, which is close to 1. This topological transition accompanies the direct gap closing at k = M and the band inversion between the first-and second-lowest positive-energy bands. To see this, we compare the particle and current densities of the lowest positive-energy band's wave function at the direct-gap closing point. Here, for the Bloch wave function of the mth band, ψ (r) = r|ψ m,k , the particle and current densities are given by The comparison of n(r) and j(r) before (χ = 0.8) and after (χ = 4) the transition point shows a drastic change in the wave function, which signifies that the band inversion has occurred in the phase transition. In the current density plot, one can also find that the circulation direction of the electron flips as the band inversion occurs. This phenomenon can also be captured in the calculation of the mth band contribution to the orbital magnetization [53][54][55], where |u m,k = e −ik·r |ψ m,k and H k = e −ik·r H eff e ik·r . In Fig. 2(b), one can see that M orb of the lowest positive band shows the sign flip at the phase transition point, agreeing with the observation in the current density plots. We also remark that even if this topological phase transition theoretically exists regardless of the Gaussian beam size, it is desirable to keep w comparable to l for experimental realizations since a fainter superlattice will imply a smaller direct band gap.
This topological phase transition could be experimentally detected in several ways. The change in C 1 causes the difference in the Hall current carried by the chiral edge state, and such difference can be revealed by transport measurements, similar to Ref. [35]. For the bulk property, one can measure the orbital magnetization, where the sudden jump would be observed at the phase transition shown in Fig. 2(b).
As the superlattice size χ increases, the electrons become localized at the local minima of |A 0 (r)|. This provides an explanation for the exponential suppression of the bandwidth of the lowest positive-energy band (δE ) in χ [ Fig. 2(c)]. For well-localized electrons, the dynamics can effectively be described by a tight-binding model, and the tunneling energy of that model is approximately given by the WKB integrals. This integral decays exponentially with the distance between the superlattice sites, so the bandwidth decreases exponentially as well. The band gaps where the details of this band-gap scaling are explained in the Appendix B.

IV. SUPERLATTICE SHEARING
To further investigate the role of the superlattice geometry, let us shear the square superlattice by angle θ so that L 1 = lx and L 2 = l (tan θx +ŷ). From the perspective of the Floquet Chern insulator created by uniform CP light, in a large superlattice size limit where the tight-binding description is valid, we might interpret the electron tunneling between superlattice sites as the chiral currents around the strongly irradiated region. That is, the paths that these chiral currents flow would give the major contribution to the path integral from one superlattice site to another. In this viewpoint, two superlattice sites can have a complex tunneling phase between them if the system has no reflection symmetry along the line connecting the two sites [ Fig. 3(a)], which is analogous to Ref. [56]. Then we can see that the tunneling terms of the tight-binding model for the square lattice (θ = 0 and θ = π/4) are real. At angles close to θ = tan −1 (1/2), the localized electrons form a hexagonal superlattice under a uniform strain and can have complex tunneling phases between the next-nearest neighbors. Then we can construct a tight-binding model for the lowest positive band similar to the Haldane model [57], as explained in the Appendix C. Similar to the Haldane model, a complex tunneling phase in the next-nearest-neighbor tunneling makes C 1 nonzero at this angle. With these considerations, we can predict successive topological phase transitions as we increase θ from 0 to π/4. We obtain the phase diagram numerically in Fig. 3(b) by calculating the Chern number of the lowest positive-energy band for each value of χ and θ . As we predicted, we can observe the successive topological phase transitions at χ larger than a certain value, which corresponds to the phase transition point described in Fig. 2. Another salient feature is that the C 1 = 1 regime very sharply blows up toward the angle θ = tan −1 (1/2), at which the χ region for C 1 = 1 diverges. This can be explained by combining the fact that the size of tunneling strengths decreases exponentially with the distance between the superlattice points and another fact that the Dirac cones can disappear and the topologically trivial gap opens in the extreme strain (see Appendix C). We can also see that the topological phase transition also accompanies the gap closing and the band inversion, as shown in the particle density plots [ Fig. 3(c)].

V. HEXAGONAL LATTICE TO KAGOME LATTICE
To engineer favorable features such as flatter bands, we can create an even more complicated superlattice by superposing different kinds of lattices. For instance, we consider the superposition of the triangular lattice beam A tri (r) and the hexagonal lattice beam rA hex (r), where r is the amplitude ratio of the two lattices (Fig. 4). When the contribution from the hexagonal lattice beam is negligible, the localized electrons form a hexagonal superlattice and the lowest part of the positive-energy spectrum can be explained by a two-band model. As r increases, electrons are confined to a kagome superlattice [58] and the lowest part of the positive-energy spectrum can be explained by a three-band model including a flat band. Note that slight gaps are observed in both the twoband model for the hexagonal superlattice and the three-band model for the kagome superlattice. The gap in the two-band model can be explained with the Haldane model with complex phases in the next-nearest-neighbor tunneling, as shown in Fig. 3(a). The gap in the kagome lattice comes from the complex phase in the nearest-neighbor tunneling [47,59]. At   FIG. 3. (a) We shear a square lattice by angle θ. Tunneling between two sites can be understood as the flow of chiral edge currents around each Gaussian CP beam. If the system has reflection symmetry around the line connecting the two sites, this tunneling should be real. Otherwise, the tunneling can have a complex phase. As examples, the next-nearest-neighbor tunnelings for the θ = 0, π/4 case and the θ = tan −1 (1/2) case are presented. (b) The Chern number of the lowest positive-energy band C 1 is shown as a phase diagram between the shearing angle θ and the superlattice size χ . (c) Energy spectra for χ = 2.4 at selected angles are shown, where the colors of low-lying bands represent the Chern numbers. The particle density in units of l −2 is plotted for angles before and after the phase transition.
FIG. 4. Superposition of the triangular lattice beam A tri (r) and the hexagonal lattice beam rA hex (r). As we increase the ratio r, we effectively change the electron superlattice from the hexagonal lattice to the kagome lattice. Energy spectra for χ = 5.4 at selected values of r are shown where the colors of low-lying bands represent the Chern numbers. By zooming in the spectrum, we can check the gaps in the two-band model and the three-band models in the lowest part of the spectrum. r = 0, we can see that the third band is nearly flat, while it is gapped well from the other bands. This flat band can be potentially used to stabilize strongly correlated phases.

VI. EXPERIMENTAL FEASIBILITY
For numerical calculation, we have set A 0 = 0.006(ea) −1 , ω = 0.06t, and w/l = 0.3 for Figs. 2 and 3. With the typical values of t = 3 eV and a = 0.142 nm for the monolayer graphene, these parameters of the laser field correspond to the field amplitude 7.6 × 10 6 V/m, the beam frequency 43.5 THz, and beam spot size 0.1 μm (FWHM). This is similar to the beam frequency in a recent experiment [35] while the peak intensity is about 4% of the beam used in the same experiment. With these parameters, the typical size of the gap ( b in Fig. 2) is 4 meV. Figure 4 uses A 0 = 0.0015(ea) −1 and ω = 0.06t, while w/l = 0.3 and w/l = 0.15 for A tri (r) and A hex (r), respectively. Finally, we remark that due to the injection of photons into the system, heating effects could eventually destroy the nontrivial topological behavior that is initially formed. Therefore, we only consider the prethermal regime where electron-electron and electron-phonon scatterings can be ignored [45]. In the past few years, the existence of this transient regime has been convincingly demonstrated in several pump-probe experiments [34,35,60].

VII. OUTLOOK
By considering Coulomb interaction in our nearly flat and topologically nontrivial bands, one could potentially induce strongly correlated phases such as fractional Chern insulators [3,[61][62][63], superconductors [7,13,[17][18][19]64,65], or magnetic phases [9,[14][15][16]. Moreover, by irradiating with frequencies comparable to the bare tunneling strength instead of the highfrequency regime considered here, higher-order terms become relevant [47] and, therefore, one can induce a wider class of  structures. While we focus on the Dirac semimetal system in this paper, our scheme can also be applied to other 2D materials such as semiconductors [66]. Our approach can be combined with other methods, such as surface acoustic waves in a solid-state platform [67], for trapping, cooling, and controlling charged particles, and for simulation of quantum many-body systems. Finally, these ideas could be used to engineer a new class of dielectric materials for potential applications in optical devices [68].

APPENDIX A: FLOQUET EFFECTIVE HAMILTONIAN IN HIGH-FREQUENCY REGIME
Let us consider the Hamiltonian given by Eq.
(1) with A(r, t ) = A 0 (r)e iωt (x + iŷ) + c.c. Then, we can write the time-dependent Hamiltonian as where σ ± = (σ x ± iσ y )/2. For this Hamiltonian, the nonzero temporal Fourier components H q = (ω/2π ) 2π/ω 0 H (t )e −iqωτ dτ are H 0 = v(τ z p x σ x + p y σ y ) and H ±τ z = 2evτ z A 0 (r)σ ± . Then, the effective Hamiltonian in the high-frequency regime is [31,[42][43][44][45][46][47] The description in terms of Eq. (A2) is valid as long as H q ω for every q. The condition ω evA 0 ensures that H q=±1 ω. For H 0 ω, we require v/l ω and the parameters we use in our paper satisfy this condition. Yet, one may wonder if the band structure is affected by the hybridization of different Floquet sidebands [36,47] since the driving frequency that we consider in this paper (ω = 0.06t) is much smaller than the original bandwidth of the graphene which is of the order of t. To see how much our . Therefore, if we consider the case that oscillating terms are slowly turned on, the relevant quasimodes should have high overlaps with the zeroth Floquet sideband, which is quantified by For the comparison of the two descriptions given by Eqs. (A2) and (A3), we calculate the band structure plotted in the left of Fig. 2(a) with these two descriptions, respectively (see Fig. 5). For the band structure calculated with Eq. (A3), we represented the overlap with the zeroth Floquet sideband, p (0) s,k , for each state. In the energy much lower than ω/2, the spectrum calculated with the high-frequency expansion and the quasienergies of the Floquet eigenstates with high overlaps with the zeroth Floquet sideband agree with each other. As the energy approaches ω/2, the Floquet sideband hybridization due to the resonant process affects the band structure. Therefore, one can use the high-frequency expansion description in Eq. (A2) for energies far smaller than the driving frequency.

APPENDIX C: TIGHT-BINDING MODEL FOR HEXAGONAL LATTICE UNDER A UNIFORM STRAIN
Let us consider the effective lattice model for the sheared lattice in the vicinity of angle θ = tan −1 0.5. By considering the terms up to the next-nearest neighbors, we can build a tight-binding model similar to the Haldane model, Here, c (A/B) † m,n creates an electron in the sublattice A or B at the unit cell (m, n) and t i=1,2,3 (s i=1,2,3 ) is the nearest-(nextnearest)-neighbor tunneling amplitude, as shown in Fig. 6. This model can be thought of as a hexagonal lattice under a uniform strain. By considering the inversion symmetry of the corresponding pairs of lattice sites, we can find that Im(t 1 ) = Im(t 2 ) = Im(t 3 ) = 0. Now we can write the Bloch Hamiltonian of this tight-binding model as where the σ z = ±1 corresponds to the sublattice A or B. In the absence of the next-nearest-neighbor tunnelings, V = h z = 0 and the location of the Dirac points is determined by h x (k) = h y (k) = 0. If this equation has two solutions, we denote those solutions as k = ±K D . In Eq. (C2), we can see that the second-neighbor tunnelings solely determine the σ z component of the Bloch Hamiltonian and do not affect the σ x and σ y components. By turning on the second-neighbor tunnelings, we effectively turn on the mass term around each of the Dirac points ±K D . Since h z (−K D ) = −h z (K D ), the sign of the effective mass term is opposite at the two different Dirac points. Then, each Dirac point equally contributes 1/2 to the Chern number, just as in the Haldane model. Therefore, if the two Dirac points exist in the absence of the next-nearest-neighbor tunneling, the lowest positive band has nonzero Chern number when the next-nearest-neighbor tunneling is turned on. Regarding this condition, the equation h x (k) = h y (k) = 0 has two solutions as long as |t i − t j | < |t k | for every i = j = k = i. Since the tunneling strength decreases exponentially in the intersite distance, t 1 , t 2 , and t 3 become very different as the superlattice size gets larger. Then there is no Dirac point after some value of χ , as shown in the square-lattice case. Yet, at angle θ = tan −1 0.5, t 1 = t 2 , so that |t i − t j | < |t k | is satisfied as long as t 3 = 0, and therefore C 1 can remain nonzero at this angle.

APPENDIX D: GAUGE-INDEPENDENT CALCULATION OF ORBITAL MAGNETIZATION
We want to numerically calculate the orbital magnetization of the mth band expressed in Eq. (6), where A m,k = 2Im(∂ k x u m,k j |)∂ k y |u m,k j is the Berry curvature. To numerically calculate it, we first need to discretize the Brillouin zone and calculate the Bloch state |ψ m,k j . Although the orbital magnetization is gauge independent, we need local gauge fixing to make |ψ m,k differentiable. While this local gauge fixing works well for smooth A m,k , it can work badly for the system in the vicinity of the topological phase transition. To avoid this subtlety, let us find a way to calculate this quantity in a gauge-independent way. For Berry curvature A m,k , a method for gauge-independent calculation is known [52]. Similar to this method, we can calculate the first integral of Eq. (D1). For this, let us consider a square patch whose four corners are q l, j ≡ k l + (δk/2)(s jx + w jŷ ), where (s 1 , w 1 ) = (−1, −1), (s 2 , w 2 ) = (1, −1), (s 3 , w 3 ) = (1, 1), and (s 4 , w 4 ) = (−1, 1). Now, u m,q l, j = u m,k l + δk 2 s j ∂ k x u m,k l + w j ∂ k y u m,k l + δk 2 8 ∂ 2 k x u m,k l + 2s j w j ∂ k x ∂ k y u m,k l +∂ 2 k y u m,k l + O(δk 3  = 1 + δkRe 4 j=1 s j u m,k l ∂ k x u m,k l + w j u m,k l ∂ k y u m,k l +δk 2 Re u m,k l ∇ 2 k u m,k l + δk 2 Re j s j w j 2 u m,k l ∂ k x ∂ k y u m,k l + δk 2 4E m,k l 4 j=1 s j s ( jmod4)+1 ∂ k x u m,k l H k l ∂ k x u m,k l + s j w ( jmod4)+1 ∂ k x u m,k l H k l ∂ k y u m,k l + w j s ( jmod4)+1 ∂ k y u m,k l H k l ∂ k x u m,k l + w j w ( jmod4)+1 ∂ k y u m,k l H k l ∂ k y u m,k l + O(δk 3 ) = 1 + 2δk 2 Re u m,k l ∇ 2 k u m,k l + i 2δk 2 E m,k l Im ∂ k x u m,k l H k l ∂ k y u m,k l + O(δk 3  = e 4π 2 Im ∂ k x u m,k l H k l ∂ k y u m,k l δk 2 + O(δk 3 ), (D3) and this corresponds to the first integral of Eq. (D1) over the square patch that we considered. One can easily check that this expression is invariant under any gauge transformation, |u m,k → exp[iλ(k)]|u m,k , ∀λ(k), and does not require any local gauge fixing.