Nuclear magnetic resonance line shapes of Wigner crystals in $^{13}$C-enriched graphene

TL;DR

This paper theoretically analyzes the NMR line shapes of Wigner crystal states in $^{13}$C-enriched graphene under strong magnetic fields, highlighting how different crystal states affect the NMR signals and the importance of linewidth in detecting these states.

Contribution

It introduces a detailed computation of NMR line shapes for various crystal ground states in graphene, considering hyperfine and orbital couplings, and assesses their detectability via NMR.

Findings

Wigner and bubble crystals have distinct NMR line shapes from uniform states if linewidth is small.

Orbital motion coupling contributes significantly to NMR frequency shifts in nonuniform states.

Small linewidths can obscure the signatures of crystal states in NMR measurements.

Abstract

Assuming that the nuclear magnetic resonance (NMR) signal from a $^{13}$C isotope enriched layer of graphene can be made sufficiently intense to be measured, we compute the NMR\ lineshape of the different crystals ground states that are expected to occur in graphene in a strong magnetic field. We first show that in nonuniform states, there is, in addition to the frequency shift due to the spin hyperfine interaction, a second contribution of equal importance from the coupling between the orbital motion of the electrons and the nuclei. We then show that, if the linewidth of the bare signal can be made sufficiently small, the Wigner and bubble crystals have line shapes that differ qualitatively from that of the uniform state at the same density while crystal states that have spin or valley pseudospin textures do not. Finally, we find that a relatively small value of the bare linewidth is sufficient to wash out the distinctive signature of the crystal states in the NMR line shape.

Equations

$$H=\xi \hslash {\nu }_F({\tau }_x{k}_x-\xi {\tau }_y{k}_y),$$

$$g{\mu }_B{B}_0=0.116B_0\text{ meV,}$$

$$⟨S_z\left(\mathbf{r}\right)⟩=\frac{1}{2}\hslash \left[⟨n_{+}\left(\mathbf{r}\right)⟩-⟨n_{-}\left(\mathbf{r}\right)⟩\right],$$

$$g\left(f\right)=A_0{e}^{-{\left(f-f_0\right)}^2/{\Gamma }^2},$$

$$I\left(f\right)=\frac{2}{N_c{A}_c}\int d\mathbf{r}g\left(f-f_0-\Delta f\left(\mathbf{r}\right)\right),$$

$$I\left(f\right)=\frac{1}{A_{WC}}{\int }_{WC}d\mathbf{r}g\left(f-f_0-\Delta f\left(\mathbf{r}\right)\right).$$

$$h_D=h_{orb}+h_{spin},$$

$$h_{orb}\left(\mathbf{r}\right)=-\frac{{\mu }_0}{4\pi }\hslash {\gamma }_n{I}_z{\left(\frac{\mathbf{r}\times {\mathbf{j}}_{\xi ,\sigma }}{r^3}\right)}_z,$$

$${\mathbf{j}}_{\xi ,\sigma }=-\xi e{\nu }_F\left[{\tau }_x\hat{\mathbf{x}}-\xi {\tau }_y\hat{\mathbf{y}}\right]$$

$$h_{spin,\xi }\left(\mathbf{r}\right)$$

$$h_{orb}\left(\mathbf{r}\right)=\xi C\frac{1}{r^3}({\tau }_{x}y+\xi {\tau }_{y}x),$$

$$C=-\frac{{\mu }_0}{4\pi }e{\nu }_F\hslash {\gamma }_n{I}_z.$$

$$H_{orb,n}\left({\mathbf{R}}_m\right)$$

$${\Psi }_{n,\xi ,\sigma }\left(\mathbf{r}\right)=\frac{1}{\sqrt{L_y}}\sum _X{e}^{-iXy/{\ell }^2}{\eta }_{\xi ,n}\left(x-X\right)c_{\xi ,n,X,\sigma },$$

$${\eta }_{-,n}\left(x\right)$$

$${\eta }_{+,n}\left(x\right)$$

$${\eta }_{-,0}\left(x\right)$$

$${\eta }_{+,0}\left(x\right)$$

$$H_{orb,n}\left({\mathbf{R}}_m\right)$$

$$D_{n,\xi ,X^{\prime },X}\left({\mathbf{r}-\mathbf{R}}_m\right)$$

$${\rho }_{\alpha ,\beta }\left(\mathbf{p}\right)=\frac{1}{N_{\phi }}\sum _{X,X^{\prime }}{e}^{-\frac{i}{2}{p}_x\left(X+X^{\prime }\right)}{\delta }_{X,X^{\prime }+p_y{\ell }^2}{c}_{\alpha }^{†}{c}_{\beta }$$

$$c_{\alpha }^{†}{c}_{\beta }=\sum _{\mathbf{p}}{\rho }_{\alpha ,\beta }\left(\mathbf{p}\right)e^{\frac{i}{2}{p}_x\left(X+X^{\prime }\right)}{\delta }_{X^{\prime },X+p_y{\ell }^2}.$$

$$\int d\mathbf{r}{\Psi }_{n,\xi ,\sigma }^{†}\left(\mathbf{r}\right)e^{-i\mathbf{p}\cdot \mathbf{r}}{\Psi }_{n^{\prime },{\xi }^{\prime },{\sigma }^{\prime }}\left(\mathbf{r}\right)$$

$$F_{n,n^{\prime }}\left(\mathbf{p}\right)$$

$$H_{orb,n}\left({\mathbf{R}}_m\right)$$

$$D_{n,\xi ,X^{\prime },X}\left({\mathbf{r}-\mathbf{R}}_m\right)$$

$$\Lambda \left(\mathbf{r}-{\mathbf{R}}_m\right)=\left(y+ix\right)-\left(y_m+i{x}_m\right)$$

$${\Phi }_{n{X}^{\prime };mX}\left(x\right)={\phi }_n\left(x-X^{\prime }\right){\phi }_m\left(x-X\right).$$

$${\phi }_n\left(x-X\right)=\frac{1}{L_x}\sum _t{e}^{it\left(x-X\right)}{\tilde{\phi }}_n\left(t\right),$$

$$\sum _{X^{\prime },X}{e}^{-i\left(X-X^{\prime }\right)y/{\ell }^2}{e}^{\frac{i}{2}{p}_x\left(X+X^{\prime }\right)}$$

$$\frac{1}{N_{\phi }}\sum _X{e}^{iX\left(p-q\right)}={\delta }_{p,q}.$$

Full text

NMR, Resistively-detected NMR, Wigner crystal, Skyrme crystal

I INTRODUCTION

In a perpendicular magnetic field $\mathbf{B}=B_{0}\widehat{\mathbf{z}}$, the eigenstates of the two-dimensional electron gas (2DEG) are quantized into Landau levels with energy $E_{n}=(n+1/2)\hslash\omega_{c},$ where $n=0,1,2,...$ is the Landau level index and $\omega_{c}=eB_{0}/m^{\ast}$ is the cyclotron frequency with $m^{\ast}$ the effective electronic mass. At integer fillings of these levels, the 2DEG exhibits the quantum Hall effectReviewQHE . However, for non-integer fillings and in certain ranges of filling factor $\nu=N_{e}/N_{\varphi}$ in each Landau level ($N_{e}$ is the number of electrons and $N_{\varphi}=S/2\pi\ell^{2}$ is the Landau level degeneracy with $S$ the area of the 2DEG and $\ell=\sqrt{\hslash/eB_{0}}$ the magnetic length), the electrons can form a Wigner crystalWigner * *i.e., a crystal of electrons. In Landau levels $n>1$, the Coulomb interaction can also favor the formation of bubble crystals (i.e., Wigner crystals with more than one electron per siteRevuebulles ; Wignerrevue ) or, near half-integer filling, a the formation of a stripe phaseStripe ). In GaAs quantum wells, these crystal states have been detected by a series of transport and microwave absorption experimentsMicrowave . More recently, resistively-detected nuclear magnetic resonance (RD-NMR) has been used to study the Wigner and bubble crystals in Landau levels $n=0,1$ in these quantum wellsTiemann . This technique uses the strong hyperfine contact interaction between the electron and the nuclear spins in GaAs quantum wells to obtain a signal. In crystal states, the spatial modulation of the electronic density leads to a modulation of the hyperfine coupling and thus to a modulation of the nuclear frequency shift (the Knight shift). In consequence, the effective electronic Zeeman energy, which is what is ultimately detected in RD-NMR, also varies spatially. The resulting NMR signal is a reflection of the electronic density distribution in the crystalTiemann . From the spectral line shape of the NMR signal, it is thus possible to distinguish between different types of crystal statesTiemann ; CoteNMR .

Magnetically induced Wigner crystal were also predicted to occur in grapheneCoteSkyrmeGraphene ; GoerbigCW ; Joglekar and in multilayer grapheneYoshioka . In graphene, the Landau level energies are given by $E_{n}=\mathrm{sgn}\left(n\right)\sqrt{2\hslash v_{F}^{2}e\left|n\right|B_{0}},$ where the Landau level index $n=0,\pm 1,\pm 2,...$ and the Fermi velocity $v_{F}=3c\gamma_{0}/2\hslash,$ with $c=1.42$ Å the separation and $\gamma_{0}=3.12$ eV the hopping energy between adjacent carbon atoms. Since the Zeeman coupling is very small, Landau levels are usually considered as fourfold degenerate when counting spin and valley degrees of freedom. The additional valley degeneracy in graphene increases the diversity of crystal phases. For example, crystals with a valley pseudospin texture become possible ground states. In this paper, we consider, from a theoretical point of view, the possibility to detect these crystals using their signal in bulk NMR. We are motivated by the chiral nature of the conduction electrons in graphene which, as was shown recentlyDora , leads to a NMR behavior that is different from that of a normal metal. For example, we show that, when nonuniform states are considered, the interaction between the nuclear spin of a carbon 13C atom and the orbital motion of the conduction electrons is nonzero. It is of the same order as the hyperfine coupling between a nuclear spin and the spin of the conduction electrons. In graphene, both couplings cause a NMR frequency shift. They must be summed to compute the total signal. We study the line shape of the crystal phases predicted theoretically and compare them with the NMR signal obtained from a uniform state with the same average density of electrons. Among the predicted ground states, we find that Wigner and bubble crystals (i.e., states with no valley and/or spin texture) have NMR line shapes that differ qualitatively from those of the corresponding uniform states in part because the orbital shift is zero in the uniform states. Thus, in principle, bulk NMR could be used to identify the presence of the crystal states in the phase diagram of the chiral 2DEG in graphene.

Nuclear magnetic resonance cannot be used for ordinary graphene since the 12C atoms have no nuclear spin and the percentage of 13C atoms (which have a nuclear spin) is very small (of the order of $1\%$). It is possible, however, to increase the fraction of 13C atoms in graphene sheets and even to fabricate pure 13C graphene. Even then, the NMR signal may be too small to allow detection by bulk NMR since the number of nuclei is small in a two-dimensional system unless a prohibitively large sheet is considered. This problem was addressed previously in Ref. Dora, where it was concluded that bulk NMR detection in 13C-enriched graphene sheets should soon be within experimental reach. We will assume that such is the case and concern ourselves solely with the width and shape of the spectral line. We remark that RD-NMR would be another way to solve the intensity problem since it is more adapted to system with a small density of nuclear spinsGervais .

Although there is a number of theoretical calculations of NMR parameters in carbon nanotubes and grapheneFabio , there are at present no available NMR experimental results for a single graphene layer so that we have to rely on estimates for the NMR parameters. For the hyperfine contact and dipolar interactions, we use the values given in Ref. Dora, . For the bare NMR linewidth which is also unknown, our calculation allows us, by comparing the signal from the crystal and uniform state, to estimate the maximal value above which it is no longer possible to distinguish a given crystal state from the corresponding uniform state.

Our paper is organized as follow. We explain in Sec. II how the NMR line shape is calculated when crystal states are considered. In Sec. III, we derive the Hamiltonian for the coupling between the nuclear and electron spin and summarize our calculation for the coupling between the nuclear spin and the orbital motion of the electrons. In Sec. IV, we review how the calculation of the Fourier components of the electronic and spin densities of the crystals is done. The line shapes of the predicted crystal states in graphene are calculated in Sec. V. We conclude in Sec. VI.

II NMR LINE SHAPE FOR CRYSTAL STATES

Graphene is a single sheet of carbon atoms arranged in a honeycomb lattice which can be described as a triangular Bravais lattice with a basis of two carbon atoms denoted $A$ and $B$. The triangular lattice constant $a_{0}=\sqrt{3}c,$ where $c=1.42$ Å is the separation between to adjacent carbon atoms. The Brillouin zone of the reciprocal lattice has $6$ valleys, but only two of them are inequivalent. We use $\xi=\pm 1$ for the valley index and choose the two inequivalent valleys to be at wave vectors $\mathbf{K}_{\xi}=\frac{2\pi}{a_{0}}\xi\left(\frac{2}{3},0\right)$. In the sublattice basis $\left(A,B\right),$ and for wave vectors $\mathbf{k}$ near $\mathbf{K}_{\xi}$ (i.e., in the continuum approximation), the noninteracting electronic Hamiltonian is given by

[TABLE]

where $v_{F}=3c\gamma_{0}/2\hslash=9.\,\allowbreak 06\times 10^{5}$ m$\cdot$s*-1* is the Fermi velocity with $\gamma_{0}=2.8$ eV the nearest-neighbor hopping energy. In Eq. (1), $\tau_{x},\tau_{y}$ are Pauli matrices in sublattice space (we later use $\mathbf{\sigma}$ for the Pauli matrices in spin space). It follows from Eq. (1) that in the continuum approximation, electrons are massless and chiralGrapheneReview .

In a perpendicular magnetic field $\mathbf{B}=B_{0}\widehat{\mathbf{z}},$ the kinetic energy is quantized into Landau levels with energy $E_{n}=\mathrm{sgn}\left(n\right)\sqrt{2\hslash v_{F}^{2}e\left|n\right|B_{0}}$ where the Landau level index $n=0,\pm 1,\pm 2,...$ and $\mathrm{sgn}$ is the signum function. Each Landau level $n$ contains four sublevels with quantum numbers $\left(\xi,\sigma\right),$ where $\sigma=\pm 1$ is the spin index. In the absence of Coulomb interaction, these four levels are almost degenerate owing to the smallness of the Zeeman coupling which is given by

[TABLE]

with $B_{0}$ in Tesla and a g-factor given by $g=2.$ For Landau level $n=1$, $E_{1}=32.88\sqrt{B_{0}}$ meV while the Coulomb energy $e^{2}/4\pi\varepsilon_{r}\varepsilon_{0}\ell=22.91\sqrt{B}$ meV for $\varepsilon_{r}=2.45$ appropriate for graphene on a SiO2 substrateJang .

We denote the filling factor of level $n$ by $\nu_{n}=2\pi n_{e}\ell^{2}\in\left[0,4\right],$ where $n_{e}$ is the averaged electronic density *in this level. *The average spin density (in the direction of the applied magnetic field) is given by

[TABLE]

where $\left\langle n_{\sigma}\left(\mathbf{r}\right)\right\rangle$ is the areal density of each spin species and $\left\langle\cdots\right\rangle$ denotes a ground state average (we set $T=0~{}$K in all our calculations).

Following previous worksCoteNMR ; Tiemann , we write $g\left(f\right)$ for the intrinsic line shape of the 13C nucleus and assume that it has the Gaussian form

[TABLE]

where $\Gamma$ is the intrinsic linewidth, $f_{0}$ the bare resonance frequency and $A_{0}$ the amplitude. Since we are not interested in the absolute amplitude of the signal, we set $A_{0}=1$ in our calculations. The 13C isotope has spin $I=1/2$ and its nuclear magnetic moment is given by $\mathbf{\mu}_{n}=\hslash\gamma_{n}\mathbf{I}_{n},$ where the gyromagnetic factor $\gamma_{n}/(2\pi)=10.705$ MHz/T. In an external magnetic field $B_{0},$ the bare nuclear precession frequency is $f_{0}=\left(\gamma_{n}/2\pi\right)B_{0}.$ This frequency is shifted by the interaction $h_{D}$ between the 13C nuclei and the electron gas. The difference between the measured and bare frequencies is $\Delta f\left(\mathbf{r}\right)$. We compute the interaction $h_{D}$ in the next section and show that the resulting frequency shift is proportional to the electronic density for the coupling with the orbital motion of the electrons and to the spin density for the Knight shift. The line shape $I\left(f\right)$ of the crystal state is obtained by summing the contribution from all nuclei in the graphene sheet. Since we work in the continuum approximation, we assume that these nuclei form a continuous background and so the intensity per carbon atom is given by

[TABLE]

where $A_{c}$ is the area of a primitive unit cell of the graphene lattice (which contains two carbon atoms), $N_{c}$ is the total number of carbon atoms and the integral is over the graphene sheet. The spatial variations of the electronic density in a non uniform state leads to an inhomogeneous broadening of the intrinsic line shape. In this paper, we concentrate on Wigner crystal phases so that $\Delta f\left(\mathbf{r}\right)$ is a periodic function. It is then only necessary to integrate over a unit cell area ($A_{WC}$) of the Wigner crystal i.e.,

[TABLE]

The frequency shift $\Delta f\left(\mathbf{r}\right)=\Delta f_{orb}\left(\mathbf{r}\right)+\Delta f_{spin}\left(\mathbf{r}\right)$ has two contributions, the first one, $\Delta f_{orb}\left(\mathbf{r}\right),$ comes from the orbital motion of the electrons while the second one, $\Delta f_{spin}\left(\mathbf{r}\right),$ comes from the interaction between the electronic spin density and the spin of the nuclei.

III ELECTRON-NUCLEUS INTERACTIONS

The interaction Hamiltonian $h_{D}$ is the sum of the two contributions

[TABLE]

where $h_{orb}$ is the coupling between a nucleus and the orbital motion of an electron and $h_{spin}$ is the coupling between a nucleus and the electronic spin density around it. The hyperfine coupling $h_{spin}$ is itself the sum of a dipolar and a Fermi hyperfine contact interactions. Although the $\pi-$bands conduction electrons in graphene occupy $p_{z}$ orbitals, the contact interaction is nevertheless not zero because of contributions from the valence $\sigma-$band states which are due to exchange-polarization effectsYazyev .

Due to the massless character of the electron in graphene, the standard textbook formSlichter of the hyperfine interaction $h_{orb}$ needs to be reformulated. As shown in Ref. Dora, , the coupling between the magnetic moment $\hslash\gamma_{n}I_{z}$ of a nuclei at $\mathbf{R}=0$ and the magnetic field $\mathbf{B}_{orb}\left(0\right)$ due to the orbital motion of one electron is given in graphene by

[TABLE]

where

[TABLE]

is the electric current of an electron with valley and spin indices $\xi,\sigma.$

The Hamiltonian $h_{spin}$ for an electron at $\mathbf{r}$ and an ion at $\mathbf{R}=0$ is unchanged with respect to its standard form and is given bySlichter

[TABLE]

where (and in this expression only) $\mathbf{r}$ is a three-dimensional vector and $\delta\left(\mathbf{r}\right)$ the three-dimensional Dirac delta function. Moreover $\mathbf{B}_{n}$ is the magnetic field of the nucleus of spin $\mathbf{I}$, $\mathbf{S}_{\xi}$ is the spin of an electron in valley $\xi$ and $\mu_{0}$ is the vacuum permeability. The Bohr magneton $\mu_{B}=e\hslash/2m_{e},$ where $m_{e}$ is the bareYafet electronic massSlichter . The electron magnetic moment is $\mathbf{\mu}_{e}=-g\mu_{B}\mathbf{S}$ where $S=1/2$ is the electronic spin ($\mathbf{I}$ and $\mathbf{S}$ are unitless in our notation). The first term in the square brackets in Eq. (III) is the isotropic hyperfine Fermi-contact interaction while the second term is the dipole interaction.

III.1 Interaction between electronic orbital motion and nuclear spin

In Ref.Dora, , it was shown that the orbital coupling does not contribute to the frequency shift for a uniform state in zero magnetic field. In this section, we show that there is a contribution from this coupling for a non uniform state and in the presence of a transverse magnetic field. Our starting point is the HamiltonianDora of Eq. (8) for the coupling of the magnetic moment $\hslash\gamma_{n}\mathbf{I}$ of a nuclei at position $\mathbf{R}=0$ to the magnetic field $\mathbf{B}_{orb}\left(0\right)$ due to the orbital motion of one electron. We then rewrite this Hamiltonian in second quantized form and generalize it to nonuniform states.

Using Eq. (9) for the electric current in graphene, the orbital coupling can be written as

[TABLE]

where the constant

[TABLE]

Assuming that Landau level mixing can be neglected, the second quantized form of $h_{orb,\xi}\left(\mathbf{r}\right)$ for the interaction between the magnetic moment of one nucleus at $\mathbf{R}_{m}$ and all conduction electrons in Landau level $n$ is given by

[TABLE]

where the field operator $\Psi_{\xi,n,\sigma}(\mathbf{r})$ for an electron in Landau level $n$ with valley index $\xi$ and with spin $\sigma$ is given by [we use the Landau gauge $\mathbf{A}=\left(0,B_{0}x,0\right)$]

[TABLE]

where $n$ and $X$ are the Landau level and guiding-center indices respectively and $L_{y}$ is the length of the 2DEG in the $y$ direction. The spinors (in the sublattice basis) are given by

[TABLE]

for $n>0.$ In these expressions, $\varphi_{n}\left(x\right)$ are the eigenstates of the one-dimensional harmonic oscillator. For the special case $n=0,$ the spinors are

[TABLE]

For the orbital coupling, we then have

[TABLE]

where

[TABLE]

We can use the operators $\rho_{\alpha,\beta}\left(\mathbf{p}\right)$ defined by

[TABLE]

to write the product

[TABLE]

The relation between the operators $\rho_{\alpha,\beta}\left(\mathbf{p}\right)$ and the field operators is

[TABLE]

where the function

[TABLE]

with $L_{n}^{m}\left(x\right)$ the generalized Laguerre polynomials. The average values $\left\langle\rho_{\alpha,\beta}\left(\mathbf{p}\right)\right\rangle$ can be interpreted as the order parameters of the crystal state. Inserting Eq. (30) in Eq. (III.1), we get

[TABLE]

Because of the particular form of the spinors in Landau level $n=0,$ we have immediately $H_{orb,n=0}=0.$ Thus, we need only evaluate the last term in Eq. (III.1) for $n\neq 0.$ We get

[TABLE]

where the function

[TABLE]

with $\mathbf{R}_{m}=\left(x_{m},y_{m}\right)$ and

[TABLE]

We define the Fourier transform

[TABLE]

so that we can write

[TABLE]

using the identity

[TABLE]

With the change of variables $\mathbf{r}\rightarrow\mathbf{r}+\mathbf{R}_{m},$ we get

[TABLE]

where $\Lambda\left(\mathbf{r}\right)$ has now been redefined as

[TABLE]

Another change of variables: $t\rightarrow t+p_{x}$ in the second line of $H_{orb,\xi,n}$ gives

[TABLE]

where

[TABLE]

and $\theta_{\mathbf{p}}$ is the angle between the vector $\mathbf{p}$ and the $x$ axis.

The function $\Upsilon_{\left|n\right|>0}\left(\mathbf{p}\right)$ can be written as

[TABLE]

where

[TABLE]

Using the definitions of the normalized eigenfunctions for the one-dimensional harmonic oscillator:

[TABLE]

we get for the functions $\Upsilon_{n}$

[TABLE]

In a uniform state, the average $\left\langle\rho_{\xi,n,\sigma;\xi,n,\sigma}\left(\mathbf{p}\right)\right\rangle\neq 0$ for $\mathbf{p}=0$ only and $\Upsilon_{\left|n\right|>0}\left(\mathbf{p}=0\right)=0$ so that the orbital shift vanishes in this case.

The orbital coupling can be written in terms of an effective magnetic field $B_{z,eff}$ defined by

[TABLE]

This additional magnetic field leads to a frequency shift

[TABLE]

where the constant

[TABLE]

We remark that the filled levels below Landau level $n$ do not contribute to the orbital shift since a filled state implies a uniform electronic density and $\Upsilon_{n}\left(\mathbf{p}=0\right)=0$. Also $\Delta f_{orb,n}$ depends on the sign of $n$ as expected since the current can be viewed as assumed by the holes for $n<0$ and by the electrons for $n>0.$

III.2 Interaction between the electronic spin and the nuclear spin

The coupling between the electronic spin and the spin of an ion at $\mathbf{R}_{m}=0$ is given by

[TABLE]

where

[TABLE]

is the isotropic Fermi-contact interaction and the dipole contribution is given by

[TABLE]

The term $\mathbf{a}^{\mathrm{dip}}$ in Eq. (57) is a traceless matrix given, using Eq. (III) (with $\mathbf{r}$ a 3d vector in these two equations), by

[TABLE]

where $i=x,y,z$.

The many-body spin coupling between all electrons in Landau level $n$ and one ion at $\mathbf{R}_{m}$ is

[TABLE]

The coupling $H_{spin,n}$ is sensitive to the microscopic electronic spin density near each nucleus, something that is not captured by the continuum approximation. For example, the spin coupling depends on the periodic part $u_{n,\mathbf{k}}\left(\mathbf{r}\right)$ of the Bloch function evaluated at the position of the nucleusSlichter . Thus, the spin coupling must be obtained from first principle calculations. One such calculation is done in Ref. Yazyev, for carbon nanotubes and graphene fragments. In this reference, a local axial symmetry is assumed so that the averaged $a_{xx}^{\mathrm{dip}}=a_{yy}^{\mathrm{dip}}=-\frac{1}{2}a_{zz}^{\mathrm{dip}}$ and the other averaged components of the $\mathbf{a}^{\mathrm{dip}}$ matrix are zero. If the nuclear spins are assumed completely polarized along the external magnetic field, then

[TABLE]

where $\left\langle S_{z,n}\left(\mathbf{R}_{m}\right)\right\rangle$ is the average effective electronic spin around each carbon atom from electrons in Landau level $n$ (including both valleys) and the nuclear spin $I_{z}=1/2$.

In our numerical calculations, we use the values of $A^{\mathrm{iso}}$ and $A_{zz}^{\mathrm{dip}}$ quoted in Ref. Dora, which are based on the calculations of Ref. Yazyev, . These values are $A_{zz}^{\mathrm{dip}}=146$ MHz and $A^{\mathrm{iso}}=-44$ MHz so that $A_{zz}^{\mathrm{spin}}=102$ MHz.

For a Wigner crystal, the spin polarization varies spatially and is much less that $S_{z}=1/2$ around a single carbon atom. The unit cell area of the graphene lattice is $A_{c}=3\sqrt{3}c^{2}/2=5.\,\allowbreak 24\times 10^{-20}$ m2. The area occupied by one carbon atom is $A_{c}/2$ and so $S_{z,n}$ in Eq. (III.2) is given by

[TABLE]

or, for an ion at $\mathbf{R}_{m},$

[TABLE]

where

[TABLE]

and $\mathbf{G}$ is a reciprocal lattice vector of the crystal.

The spin coupling is finally given by

[TABLE]

where the functions

[TABLE]

Using Eq. (52) to define an effective magnetic field for the spin coupling, we arrive at the following result for the Knight shift at the position of a nuclei:

[TABLE]

with the constant

[TABLE]

For graphene, the ground state of the uniform 2DEG has $\left\langle\rho_{z,n}\left(0\right)\right\rangle=-1$ for $\nu_{n}\in\left[0,2\right],$ where $\nu_{n}$ is the filling factor of Landau level $n,$ so that the shift $\Delta f_{spin}$ is positive when $\alpha_{spin}$ is given by Eq. (69).

IV ORDER PARAMETERS OF THE CRYSTAL STATES

We summarize in this section the calculation of the order parameters $\left\langle\rho_{\beta,\alpha}\left(\mathbf{G}\right)\right\rangle$ defined in Eq. (29) for the crystal phases. The formalism for this calculation is described in detail in Ref. CoteSkyrmeGraphene, for the case where the spin degree of freedom is neglected. In this paper, we include both valley and spin states. To simplify the notation, we use $\alpha,\beta$ for the combined valley and spin indices. All quantities are evaluated in Landau level $n$ and so we leave this index implicit hereafter.

The order parameters are related to the single-particle Matsubara Green’s function by

[TABLE]

where the Fourier transform of the Green’s function is given by

[TABLE]

with the Green’s function defined by

[TABLE]

with $T_{\tau}$ the time-ordering operator and $\tau$ the imaginary time.

The Green’s function is evaluated in the Hartree-Fock approximation with an Hamiltonian given by (with $\kappa=4\pi\varepsilon_{r}\varepsilon_{0}$)

[TABLE]

The term $E_{\alpha}=\sigma g\mu_{B}B/2$ is simply the Zeeman energy. The equation of motion for $G_{\alpha,\beta}\left(\mathbf{G,}i\omega_{n}\right),$ where $\omega_{n}$ are the Matsubara frequencies, is

[TABLE]

with the potentials

[TABLE]

and the Hartree and Fock interactions in Landau level $n$ given by

[TABLE]

where the form factor

[TABLE]

Eq. (74) is solved iteratively until convergence is reached for the order parameters.

V NMR SPECTRA OF SOME CRYSTAL PHASES

Because of the finite Zeeman coupling, the noninteracting uniform 2DEG in graphene is fully spin polarized for filling factor $\nu_{n}\in\left[0,2\right]$ and partially spin polarized for $\nu_{n}\in\left[2,4\right].$ Levels of same spin but different valley indices are degenerate. In graphene, the electron spin and magnetic moments are antiparallel so that the two lower noninteracting energy levels have spin $\sigma=-1.$

In this paper, we use a pseudospin language where the two valley states $K_{+},K_{-}$ are mapped onto the up and down valley pseudospin states respectively. It follows that the pseudospin vector field in real space is obtained from

[TABLE]

where the function $F_{n,n}\left(\mathbf{G}\right)$ is defined in Eq. (III.1) and the Cartesian components of $\mathbf{P}\left(\mathbf{r}\right)$ are given by

[TABLE]

For the spin field in real space, we have in the same way

[TABLE]

Finally, for the electronic density

[TABLE]

In the Hartree-Fock approximation, the electrons condense into distinct quantum Hall ferromagnetic states at filling factors $\nu_{n}=1,2,3,4$ and the (approximate) fourfold degeneracy is fully lifted. The ground state at integer filling is foundDoretto by (1) maximizing the spin and then (2) maximizing the valley pseudospin $\left|\mathbf{P}\right|$ to the extent allowed by rule (1). Thus, for the uniform 2DEG, the lower(upper) energy levels have spin per electron $s_{z}=-1/2$($s_{z}=+1/2$). Since there is no effective Zeeman coupling for the valley pseudospin $\mathbf{P}$, this vector has a full SU(2) symmetry. For $\nu=0$ (i.e., $\nu_{0}=2$), Hartree-Fock does not explain the experimental observations and a more sophisticated analysis that includes the electron-phonon interaction is neededKharitonov .

In the 2DEG, the Coulomb interaction in the presence of a quantizing magnetic field can also favor the formation of crystal states such as the Wigner crystalWigner . These states have been studied in a number of papers for grapheneCoteSkyrmeGraphene ; GoerbigCW ; Joglekar as well as for multilayer grapheneYoshioka . Several crystalline ground states were identified in Landau levels n=0,1,2,3:\triangular Wigner and bubble crystalsRevuebulles ; Cotebubble of electrons or holes and square crystals with a meronlike valley pseudospin texture. The spin degree of freedom was not considered in these previous works so that all these crystalline phases are fully spin polarized. In the present paper, we explicitly consider valley and spin degrees of freedom. For small values of the Zeeman coupling, this allows for the addition of a new phase in the phase diagram which is a crystal of spin skyrmions.

To compute the NMR line shape $I\left(f\right)$ of the crystal phases, we need a value for the linewidth $\Gamma$ in Eq. (4) i.e., the linewidth of the NMR line shape for a uniform state of the 2DEG in a 13C graphene sheet. Although there is a number of theoretical calculations of NMR parameters in carbon nanotubes and grapheneFabio , there are at present no available NMR experimental results for a single graphene layer to validate them. If we assume a linewidth limited only by the nuclear dipole-dipole interaction, then, for nuclei that are separated by a distance $r$ and have magnetic moment $\mu=\hslash\gamma_{n}I,$ a rough estimateSlichter of the local magnetic field they produce is $B_{loc}=\mu/r^{3}$ in cgs units. Since this local field can increase or decrease the external field, the spread in frequency that results is $\Delta f=2\left(\gamma_{n}/2\pi\right)B_{loc}$ for the FWHM. For graphene, this leads to a value $\Gamma\approx 1$ kHz independent of the magnetic field. Our strategy is to start with this value of $\Gamma$ and increase it until we find the maximal value above which there is no difference between the line shape of the crystal and that of the uniform state with the same average filling factor.

The intensity of the signal is of course of major importance experimentally. This issue was discussed in Ref. Dora, . In this paper, the authors estimated the limit of detection (LOD) parameter for a fully 13C enriched graphene layer using the results available for enriched fullerenes and carbon nanotubes. They concluded that, with the improvement in the synthesis of a fully 13C enriched graphene, NMR detection should be within reach. Assuming this to be the case, we will not concern ourselves with this problem in this paper.

In all the NMR line shapes shown in the figures below, the frequency is measured with respect to the resonance frequency $f_{0}.$ For comparison, we give the line shape of the crystal and that of the uniform state with the same filling factor. For a uniform state (US), the line shape is

[TABLE]

since $\Delta f_{orb}=0$ when there is no density modulation. Also, in such a state, $\Delta f_{spin}=-\alpha_{spin}B_{0}\left\langle\rho_{z,n}\left(0\right)\right\rangle$ is directly proportional to the spin polarization. In our calculations, we first set the Zeeman coupling $\Delta_{Z}.$ This, in turn, determines the value of the magnetic field $B_{0}$ and the electronic density when the filling factor is specified.

Figure 1 shows the NMR line shape of different Wigner crystals obtained with the Zeeman coupling (expressed in units of $\left(e^{2}/\kappa\ell\right)$) $\Delta_{Z}\equiv g\mu_{B}B/\left(e^{2}/\kappa\ell\right)=0.03.$ Figure 1(a) shows how the line shape changes with Landau level index $n=0,1,2$ for $\nu_{n}=0.2$ where the ground state (in the Hartree-Fock approximation) is an electron Wigner crystal. Note that the line shape of the uniform state is the same for all three crystals. Also, there is no orbital shift for $n=0$ while the signal for the other two Landau levels is a combination of spin and orbital shifts. For $n=0,$ the line shape is similar to that of the $n=0$ Wigner crystal in GaAs quantum wells obtained from RD-NMR. This is expected because the form factor function $F_{0,0}\left(\mathbf{G}\right)$ is the same in graphene and GaAs quantum wells when $n=0$ and the orbital shift is negligible in GaAsCoteNMR ; Tiemann . Figure 1 shows that the line shape changes qualitatively with Landau level index. In all cases shown in Fig. 1(a), the signal of the crystal is clearly distinguishable from that of the uniform state. This is also true for the other types of Wigner crystals found in graphene. For example, Fig. 1(b) shows the line shapes, in level $n=2,$ of the Wigner crystal at $\nu_{2}=0.2$ and of the bubble crystal (a Wigner crystal with two electrons per site) at $\nu_{2}=0.4$ while Fig. 1(c) shows the line shape, in level $n=1,$ of the Wigner crystal at $\nu_{1}=0.2$ and of the hole Wigner crystal at $\nu_{1}=3.8$ obtained by adding a filling fraction $0.2$ of holes to a fully filled (i.e., $\nu_{1}=4$) Landau level. All line shapes are very distinct from that of the corresponding uniform states. The crystals line shape tend to have long frequency tails and much larger width than the uniform states (almost $\approx 14$ kHz in some cases). However, this value is small in comparison with the width of the crystal line shape in GaAs quantum wellsTiemann ; CoteNMR which is about $30$ kHz. The hyperfine coupling is much larger in GaAs quantum wells than in graphene.

Contrary to the situation in GaAs quantum wells, there is no simple relationship between the histogram of density (or, more precisely, polarization) and the NMR line shape for the crystal states in graphene. This is due to the fact that orbital and spin shifts are of the same order and also to the fact that the orbital shift is not directly related to the electronic density or polarization. The relative importance of the spin and orbital shifts is shown in Fig. 2 where we plot the line shape obtained by keeping only the spin shift, only the orbital shift, and both shifts in the function $\Delta f\left(\mathbf{r}\right)$ in Eq. (6) for the Wigner crystal at $\nu_{1}=0.2.$ It follows that we cannot see sharp differentiating features between, say, the electron and hole crystals shown in Fig. 1(c). Moreover, we cannot give a simple explanation for the difference between the Wigner and bubble crystals in Fig. 1(b).

For $\nu_{n}<1,$ only one valley (with spin down) is occupied by the electrons. In GaAs quantum wells, the ground state in $n=0$ around $\nu_{0}=1$ is a crystal of spin skyrmions. In graphene, because of the valley degeneracy and the finite value of the Zeeman coupling, valley pseudospin skyrmions instead of spin skyrmions are expected to populate the ground state for $\nu_{0}$ or $\nu_{1}$ near $1$. Actually, numerical calculations show that merons (half a skyrmionSondhi ; Revueskyrmion ) are preferred over skyrmions when a crystal is formed so that for $\nu_{0}>0.5$, the electrons occupy both valleys leading to a state with complete valley pseudospin depolarization. Such a meron crystal state is depicted in Fig. 1 of Ref. CoteSkyrmeGraphene, for the case $\nu_{0}=0.8$ and reproduced in Fig. 4(b) below. The merons are arranged in a checkerboard configuration with four merons per unit cell. Nearest neighbors have opposite vorticity and opposite orientation of the valley pseudospin component $P_{z}$ at their center. The electronic density is a two-dimensional charge density wave. The amplitude of the density modulation is small. For example, for $\nu_{0}=0.6$, the density varies between $0.5/2\pi\ell^{2}$ and $0.7/2\pi\ell^{2}.$ Since electrons occupy the spin down levels only, the modulation of the spin component $S_{z}$ is also small. Consequently, the NMR line shape of the meron crystal is not that different from that of the uniform state as can be seen in Fig. 3(a). To see the detail of the modulations, a smaller value of $\Gamma$ is necessary. Fig. 3(b) shows the line shape from the same crystal obtained with $\Gamma=0.2,\Gamma=1$ and $\Gamma=2$ kHz. The distinct character of the meron crystal is visible for $\Gamma=0.2$ kHz while crystal and uniform states have the same line shape for $\Gamma=2.$ In the Wigner crystal (electron or hole) considered in Fig. 1, the spread in frequency is larger and it is only for $\Gamma\approx 8$ kHz that the difference between the two signals is lost. We remark that the same checkerboard type of meron crystal is found at filling factor $\nu_{1}=0.8$ and the line shape is again close to that of the uniform state. In higher Landau levels, this meron crystal disappears. For example, at $\nu_{2}=0.8,$ the ground state is a hole Wigner crystal. Our numerical approach, however, does not allow us to compute the crystal state for $\left|\nu_{n}-m\right|\lesssim 0.5$ for $m=1,2,3$ so that it is possible that the meron crystal phase exists in a very narrow range around integer filling. For isolated spin skyrmions, it was indeed shown in Ref. Luo, that the filling factor range over which they exist in $\left|n\right|=1,2,3$ decreases rapidly with Landau level index.

A new crystalline ground state that we found when considering the spin degree of freedom is a crystal with a skyrmionic spin texture. The electronic density $n\left(\mathbf{r}\right)$ and the spin texture in the $x-y$ plane of this crystal are shown in Fig. 4 where $\Delta_{Z}=0.01,n=1$ and $\nu_{1}=2.2.$ In this state, level $\left(K_{+},-\right)$ is fully filled while levels $\left(K_{-},-\right)$ and $\left(K_{-},+\right)$ are only partially filled, the former being more occupied than the latter (an equivalent ground state is obtained by interchanging $K_{+}$ and $K_{-}$). The spin texture occurs in valley $K_{-}$ . In this state, there are two vortices (the red square density patterns in Fig. 4) with the same vorticity but differing by a global rotation of $\pi$ of the in-plane spin vector. This configuration resembles the ”square lattice antiferromagnetic (SLA)” Skyrme crystal found in GaAs quantum wells around $\nu_{0}=1$BreySkyrme ; Barrett with the difference that the density profile on each site now reflects the character of the two orbitals that enters the $n=1$ spinor in graphene which is given in Eq. (17). A crystal with a spin texture is possible only if the Zeeman coupling is small. Figure 5 shows the NMR line shape for the spin Skyrmion crystal in $n=1$ and at $\nu_{1}=2.2$ at different values of the Zeeman coupling. For $\Delta_{Z}=0.03,$ [Fig. 5(a)], the spin texture is lost and the square lattice of Fig. 4 is replaced by a triangular lattice of spin up electrons with filling factor $0.2$ on top of a uniform state (with filling $2.0$) of spin down electrons. The corresponding NMR signal is similar to that of the Wigner crystals shown in Fig. 1. At $\Delta_{Z}=0.02$ and $\Delta_{Z}=0.01,$ the ground state is the spin skyrmion crystal of Fig. 4. The corresponding line shapes are shown in Fig. 5. The line shape gets closer to that of the corresponding uniform state (the dashed lines in Fig. 5) when the Zeeman coupling is reduced. At $\Delta_{Z}=0.02,$ Fig. 5 shows the line shape of both the triangular lattice (which is not the ground state) and the spin skyrmion crystal ground state. They are almost identical.

Figure 2 above shows that the spin and orbital shifts are equally important in graphene according to our calculation which is based on a graphene Hamiltonian written in the continuum approximation. This approximation, however, does not capture the rapid oscillations of the electron wave function near the carbon atoms. A first principle calculation would probably be necessary in order to estimate the validity of our approach for the orbital shift. We show in Fig. 6 the NMR line shape of the crystals considered in Fig. 1(a) when only the spin shift is kept. The signal from the crystal in $n=0$ is unchanged since there is no orbital shift in this case. For $n=1,2$ the broadening of the line shape is reduced and the signal is closer to that of the uniform state. Clearly, in graphene, the orbital shift is important in order for the NMR technique to distinguish the crystal from the uniform state.

The Hartree-Fock interaction in Eq. (77) is independent of the sign of the Landau level index so that the phase diagram in levels $n<0$ is identical to that in levels $n>0.$ Although the orbital shift in Eq. (III.1) depends on the sign of the Landau level index, the line shape of the crystals in levels $n>0$ can be simply related to the line shape of the hole crystal in levels $n<0.$ For example, the electron Wigner crystal at $\nu_{1}=0.2$ has the same line shape as the hole Wigner crystal at $\nu_{-1}=4-\nu_{1}=3.8.$ The spin polarization is the same in both crystals so that the spin shift is the same. The orbital shift has also the same sign since the current is carried by the holes for the crystal in $n=-1.$

We conclude this section by noticing that, for Landau level $n=1$ in graphene, the relative nuclear frequency shift for an ion at $\mathbf{R}_{m}=0$ can be written as

[TABLE]

using Eq. (III.1). It is interesting to compare this result with that for the 2DEG in GaAs. In this 2DEG, the orbital coupling hamiltonian is

[TABLE]

where $\mathbf{p}=m\mathbf{v}$ is the electron momentum and $m^{\ast}=0.067m_{e}$ is the effective electronic mass. Following the same steps as in Sec. 3(a) to compute the orbital shift, we find after a lengthy calculation that

[TABLE]

We checked numerically that the expression in the square brackets in Eqs. (V),(V) is of order $1$ for the crystal states so that the orbital shift is at least an order of magnitude larger in graphene than in GaAs at $B=10$ T. For $n=0,$ however, the relative frequency shift for level $n=0$ is not zero in GaAs contrary to graphene. It is given by

[TABLE]

which is however very small in comparison with the Knight shift in the GaAs 2DEG.

VI CONCLUSION

We have shown that in graphene both spin and orbital shifts contribute equally to the NMR line shape of the crystal phases. From all the crystal phases we studied, only the Wigner and bubble crystals have a line shape that is distinctively different from that of the uniform state with the same average filling factor. If our calculation does not overestimate the contribution of the orbital shift, a measurement of the NMR line shape would be one way to detect the formation of the electron crystal in graphene. This conclusion, however, is valid as long as the linewidth of the NMR signal in the uniform state is small. For the value of the hyperfine spin coupling that we used in our numerical calculations, we estimate that the maximal value of $\Gamma$ is around $8$ kHz. Since the spin and orbital coupling parameters are not precisely known, this $\Gamma_{\max}$ must be considered as only an estimate.

We have not considered the effect of quantum or thermal fluctuations in our analysis. The importance of these fluctuations was discussed in detail in connection with resistively-detected NMR measurement (RD-NMR) of Wigner crystals in GaAs quantum wellsTiemann . It is clear that they would tend to wash out even more the difference between the line shape of the crystal and that of the uniform state.

We assume in this paper that the NMR signal from the 13C isotope-enriched graphene layer is sufficiently intense and/or the experimental setup is sufficiently sensitive to allow the measurement of the change in NMR line shape induced by the formation of the crystal states. As discussed in Sec. V, this is a challenge for bulk NMR measurements where the signal is proportional to the number of polarized nuclei which is very small in a two-dimensional system such as graphene. A technique such as RD-NMR would probably be more appropriate for this type of measurement since it is well suited for systems with a small number of nucleiGervais . Our calculations of the NMR line shapes should apply to this technique as well.

Acknowledgements.

R. C. was supported by a grant from the Natural Sciences and Engineering Research Council of Canada (NSERC) and J.M. P. by scholarships from NSERC and the Fonds de recherche nature et technologies du Québec (FQRNT). Computer time was provided by Calcul Québec and Compute Canada.