Photonic crystal for graphene plasmons

Photonic crystals are commonly implemented in media with periodically varying optical properties. Photonic crystals enable exquisite control of light propagation in integrated optical circuits, and also emulate advanced physical concepts. However, common photonic crystals are unfit for in-operando on/off controls. We overcome this limitation and demonstrate a broadly tunable two-dimensional photonic crystal for surface plasmon polaritons. Our platform consists of a continuous graphene monolayer integrated in a back-gated platform with nano-structured gate insulators. Infrared nano-imaging reveals the formation of a photonic bandgap and strong modulation of the local plasmonic density of states that can be turned on/off or gradually tuned by the applied gate voltage. We also implement an artificial domain wall which supports highly confined one-dimensional plasmonic modes. Our electrostatically-tunable photonic crystals are derived from standard metal oxide semiconductor field effect transistor technology and pave a way for practical on-chip light manipulation.


Supplementary Note 1: Near-field images at other gate voltages and regions:
Supplementary Figure 1 shows the plasmonic response of the photonic crystal at all gate voltages. The general trends are the same as Figure

Supplementary Note 3: Gold-launched SPPs:
The high gate voltage g = −90 V image in Figure 2c shows gold-launched propagating SPPs in the photonic crystal 1 . The p wavelength fringes are generated due to the interference of laser beam #1 and #2 in Supplementary Figure 5. The incoming laser beam simultaneously reaches the tip (beam #1) and gold launcher (beam #2). Laser beam #1 is back-scattered directly by the tip and registered on the detector.
Laser beam #2 first reaches gold launcher which launches SPPs of wavelength p . The SPPs travel for a distance x until it reaches the tip and scattered out for detection. The two beams interfere on the detector with a phase difference of Δ = 2 p , due to propagation phase delay in graphene. The phase delay caused by light propagating in air can be neglected as p ≪ 0 , where 0 is the free-space laser wavelength. This creates the observed oscillatory fringes in the near-field image with p wavelength.
Supplementary Figure 5: Detection scheme for gold-launched SPPs. The detected optical signal is the interference between beam #1 and #2. Beam #1 is the direct scattering from the tip. Beam #2 first reaches gold antenna and propagate in the graphene device as SPPs and finally scatter out into far field light at the location of the tip.

Supplementary Note 4: Tip-launched SPPs:
The intermediate gate voltage g = −50~− 70 V images in Figure 2c show hexagonal pattern of dark spots. This hexagonal pattern originates from coupling of the tip dipole source and the localized plasmonic mode in the photonic crystal. The z ( ) pattern launched by a point source is described by the retarded real-space Green's function at the excitation frequency , Here ( ) is the plasmonic eigen-mode wavefunction, is its eigen-frequency, the sum runs over all possible plasmonic states, and denotes the Cauchy principle value. The scattered near-field signal ( , ) is related to z directly under the tip by ( , )~( , ; ) ∝ LDOS( , ) (2) Where LDOS( , ) is the local plasmonic density of states at position and frequency . Thus, the near-field image represents a local map for LDOS.
A more realistic treatment is to model the tip as a point dipole source with the effective separation a = 30 nm (tip radius) from the sample surface. The corresponding Green's function is given by 2 where ( ) = | | 2 −2| | (4) is a form-factor. Accordingly, the scattered near-field signal ( , ) is a twicerepeated convolution of ( , ′ ; ) with the function ( ) = ∫ 3 ( ), the Fourier transform of ( ). If the spatial variation of the plasmonic Green's function is slow compared to the tip radius, then the scattered near-field signal should still contain the local information about LDOS. (For more sophisticated models that go beyond the point-dipole approximation, see ref. 2

Supplementary Note 5: Optical constants for h 11 BN
The infrared dielectric function for isotopically enriched hBN has been reported to have exceptionally narrow phonon line width already at room temperature 3 . To extract the dielectric function of the isotopically enriched hBN used in this study, we performed broadband FTIR measurement and fit the reflectance spectra using a multilayer model consisting of hBN micro-flakes, SiO2 and Si substrate 4 . The dielectric function of amorphous SiO2 and Si substrate is separately measured and can be described with multiple Lorentzians, in accordance with previous studies 1 . The dielectric function of the hBN is parameterized with the standard "TOLO" form for polar materials 3 : where ∞ , , , and is the high-frequency permittivity, frequency, the TO and LO phonon frequency and the phonon damping.
The infrared reflection spectra covering both the upper (UR) and lower Reststrahlen (LR) bands of hBN micro-flakes on SiO2/Si substrate were collected using a FTIR microscope (Bruker LUMOS). We have performed the measurement both at room temperature (T = 300 K) and low temperature (T = 60 K), see Supplementary Figure 6. In Supplementary Figure 6 the black squares are experimental data and the red dashed lines are the fitting. The fitting parameters for hBN are listed in the figure.
In the reflectance spectra, the broad peak near 1080 cm −1 is related to the SiO2 phonons and the sharper peak near 1360 cm −1 is due to hBN in-plane phonon. From the room temperature data (Supplementary Figure 6 left panel) we can confirm that the hBN understudy is of isotope-enriched 11 B by comparing both the phonon frequency and phonon damping with previous reports 3 . Remarkably, the narrow line-width for h 11 BN gives rise to almost total reflection (R = 1) at , much higher than that of natural hBN (R ~ 0.55 at 1 ). With decreasing temperature, the phonon damping is further reduced and the phonon position blue-shifted slightly, also consistent with previous results.
Supplementary Figure 6: Reflectance spectra of the hBN/SiO2/Si heterostructure at T = 300 K (left) and T = 60 K (right). In both panels the black squares are experimental data and red dashed lines are fitting using Lorentzian oscillator model. Fitting parameters related to the isotopically enriched hBN is listed as insets.

Supplementary Note 6: Band structure simulation method:
In this section, we elaborate in detail on the numerical method that is used to calculate the plasmonic band structures. The plasmonic excitation is described by the coupled integral equations of the electric potential ( , ) and the carrier density ( , ): ,where �( , , ′ ) is the effective 2D screened Coulomb interaction, �( , , ′ ) is the density response function of graphene, and ( , ) is the external drive potential 5 . If we use the other set of physical quantities, capacitance and conductivity, the equivalent description is given as: Here, we used the basic constitutive relations and the continuity equation to derive Combining the above two equations and eliminating ( , ) by substitution yields: The plasmonic oscillations are self-sustained mode satisfying ( , ) = 0. From this point and below, we will use � instead of � for compact description of plasmonic oscillations: For 2D periodic system, these quantities are evaluated only over a discrete set of momenta, = + , for Bloch states with a Bloch momentum ( is the set of reciprocal lattice vectors). Then, the self-oscillation condition can be casted into a matrix nonlinear eigenvalue problem: , where = � ∥ / ⊥ is the confinement factor of the hyperbolic hBN layers, = ( , ), and the graphene is located at = 0. This ansatz is obtained by requiring ⋅ = ⋅ ( ) = − ⋅ ( ) = 0 , thereby assuming quasi-electrostatic limit. For example, in the patterned SiO2 superlattice layer, ( ) and satisfy , where is 2D gradient. It is convenient to solve for obtain the matrix expression for the dynamic capacitance �̃� using the Gauss' law: The final expression for the dynamic capacitance can be found in ref. 6 . The discussions on the density response function �( , , ′ ) , its semiphenomenological expression in the case of periodic doping, and the non-local electron responses are well given in the supplementary material of 6 . Here, we simply provide a comparison between the local Drude model and the non-local semi-phenomenological model. For the local Drude model, the density response function is given as: , where � ( ) is the Fourier component of the Fermi-level distribution ( ) = ∑ � ( ) ⋅ . For the non-local semi-phenomenological model, extra terms follow: The second term is a local correction from the virtual inter-band transition, and the third term is the leading-order non-local correction from the intra-band transition. Supplementary Figure 7 clearly depicts that the non-local effects introduce an overall blue shifts in the band structure. Also, on a closer inspection, the blue shift at K point is greater than that at M point due to 4 dependence of the non-local term, which reduces the size of the complete bandgap. Lastly, we provide a few comments on the electrostatic simulation of the Fermilevel distribution ( ) . The local Fermi-level on graphene is given as = ℏ � , where is the carrier density. We take the carrier density distribution ( ) from electrostatic simulations using COMSOL Multiphysics, see Supplementary Figure 8, and calculate ( ) from it. For the values of the fermi velocity , we actually take into account the so-called ' -renormalization', which is a mean-field many-body effect that gives -dependent . We took a phenomenological function of ( ) from 7 , which is based on experimental measurement.
-renormalization is also experimentally studied in 8 . We simply assumed that the local Fermi-level distribution is given as Supplementary Figure 8: Electrostatic simulation with COMSOL Multiphysics. The color map represents the z-component of the electric field , the arrow plots the electric field vectors, and the white lines are electric field lines.

Supplementary Note 7: Density of states and near-field signal simulations:
The plasmonic density of states (DOS) or the spectral density/function is defined as It is computationally challenging to extract all the eigen-frequencies at every given gate voltage with its dense variation. Thus, we introduce an approximately proper quantity: , where the response matrix [ ] is given as: It is apparent that [ ] is not properly defined at the plasmonic resonances, because [A] -1 has vanishing eigenvalues. This leads to a diverging trace of [ ]. Actually, in the limit of vanishing permittivity dispersion and we ignore the nonlocal effect, we have the following exact relation: Then, near the vicinity of each resonances ( ∼ , ), it is approximated as: As we included the permittivity dispersion and the nonlocal effects, the accuracy of this approximation will decrease; however, it is still a good approximation near the resonances. In Supplementary Figure 9, we provide the simulation of the total DOS as a function of frequency along with the band structure. The total DOS is defined as , where the Block momentum summation is done over the Brillouin zone. As expected, we observe highly suppressed DOS inside the bandgap, and divergent behaviors near the saddle points in the band structure 9 . This result illustrates that our system can be used as a plasmonic platform for tunable Van Hove singularities. Supplementary , where , ( ) is the electric field distribution of an eigen-resonance. In Supplementary Figure 10, we calculated LDOS at several frequencies for the structure simulated in Supplementary Figure 9. We provide the square root of LDOS, and the eigen-mode wavefunctions are evaluated at the surface of the top hBN layer. It is clear that LDOS simulations agree with our experimental results that the bright spots feature appears at the lower bands (high gate voltages) and the dark spots feature appears at the upper bands (low gate voltages). The LDOS doesn't vanish inside the bandgap due to the finite loss ( = 1 meV = 8 cm −1 ) we introduced for the simulation. Examining more carefully the LDOS pattern inside the bandgap, we even find that the residual pattern in the bandgap appears to be the dark spots pattern, which also agrees with the experiments.