Theoretical model and simulations
The model was described in detail in ref. 5. Here we provide further details of its derivation and application to our case. Atomic units are used unless otherwise stated.
Definition of the lattice
Hexagonal boron nitride is formed from two triangular sublattices A and B, which host boron and nitrogen atoms, respectively. Extended Data Fig. 1a shows the lattice up to next-nearest neighbour atoms. The two lattice vectors can be defined as
$$\begin{array}{cc} & {{\bf{a}}}_{1}=\frac{{r}_{0}}{2}\left(3,\sqrt{3}\right)\\ {\rm{and}} & {{\bf{a}}}_{2}=\frac{{r}_{0}}{2}\left(3,-\sqrt{3}\right),\end{array}$$
(1)
where \({r}_{0}\) is the distance between nearest neighbours. The distance from atom \(j\) to atom \(i\) can be written in terms of \({r}_{0}\) and the angle \(\alpha \) between those two atoms (Extended Data Fig. 1a):
$$\left[{r}_{{ij},x},{r}_{{ij},y}\right]={\left(\sqrt{3}\right)}^{m}{r}_{0}\left[\sin {\alpha }_{{ij}},\cos {\alpha }_{{ij}},\right],$$
(2)
where \(m=0,1\) for nearest neighbours and next-nearest neighbours, respectively. The Brillouin zone of hBN is shown in Extended Data Fig. 1b, along with its high-symmetry points.
Definition of the field
We start with a bicircular field vector, which results from the combination of two counter-rotating circular fields of frequencies \(\omega \) and \(2\omega \). We define it as
$$\begin{array}{c}{{\bf{F}}}_{\circlearrowright }(t)=[{F}_{x}(t),{F}_{y}(t)]=[{F}_{\omega }\sin (\omega t)-{F}_{2\omega }\sin (2\omega t+\varphi ),\\ {F}_{\omega }\cos (\omega t)+{F}_{2\omega }\cos (2\omega t+\varphi )].\end{array}$$
(3)
The field strengths of the fundamental and second-harmonic fields are \({F}_{\omega }\) and \({F}_{2\omega }\), respectively, and \(\varphi \) is the phase delay between the two fields. From the electric field, we define the vector potential \({\bf{A}}\left(t\right)=-\int \,{\rm{d}}t\,{\bf{F}}\left(t\right)\), so that
$$\begin{array}{c}{{\bf{A}}}_{\circlearrowright }(t)=[{A}_{x}(t),{A}_{y}(t)]=\left[\frac{{F}_{\omega }}{\omega }\cos (\omega t)-\frac{{F}_{2\omega }}{2\omega }\cos (2\omega t+\varphi ),\right.\\ \left.-\frac{{F}_{\omega }}{\omega }\sin (\omega t)-\frac{{F}_{2\omega }}{2\omega }\sin (2\omega t+\varphi )\right].\end{array}$$
(4)
With this definition, both the electric field and the vector potential rotate clockwise. When we switch the helicities of the two circular fields such that \(\left[{F}_{x}\left(t\right),{F}_{y}\left(t\right)\right]\to \left[-{F}_{x}\left(t\right),{F}_{y}\left(t\right)\right]\), we obtain
$$\begin{array}{l}{{\bf{F}}}_{\circlearrowleft }(t)=[-{F}_{\omega }\sin (\omega t)+{F}_{2\omega }\sin (2\omega t+\varphi ),\\ \,\,{F}_{\omega }\cos (\omega t)+{F}_{2\omega }\cos (2\omega t+\varphi )]\end{array}$$
(5)
and
$$\begin{array}{l}{{\bf{A}}}_{\circlearrowleft }(t)=\left[-\frac{{F}_{\omega }}{\omega }\cos (\omega t)+\frac{{F}_{2\omega }}{2\omega }\cos (2\omega t+\varphi ),\right.\\ \,\,\left.-\frac{{F}_{1}}{\omega }\sin (\omega t)-\frac{{F}_{2}}{2\omega }\sin (2\omega t+\varphi )\right].\end{array}$$
(6)
In this case, both the electric field and vector potential rotate anticlockwise. Between \({{\bf{F}}}_{\circlearrowright }\) and \({{\bf{F}}}_{\circlearrowleft }\), the field has rotated in space by 180° (compare the green curves in Extended Data Fig. 2a,d). However, the orientation of the vector potential is the same for \({{\bf{A}}}_{\circlearrowright }\) and \({{\bf{A}}}_{\circlearrowleft }\) (compare the purple curves in Extended Data Fig. 2b,e).
For generality, let us also consider the bicircular field defined as
$$\begin{array}{c}{\mathop{{\bf{F}}}\limits^{ \sim }}_{\circlearrowright }(t)=[{\mathop{F}\limits^{ \sim }}_{x}(t),{\mathop{F}\limits^{ \sim }}_{y}(t)]=[{F}_{\omega }\cos (\omega t)+{F}_{2\omega }\cos (2\omega t+\varphi ),\\ \,\,\,\,-{F}_{\omega }\sin (\omega t)+{F}_{2\omega }\sin (2\omega t+\varphi )]\end{array}$$
(7)
and, consequently,
$$\begin{array}{l}{\mathop{{\bf{A}}}\limits^{ \sim }}_{\circlearrowright }(t)=[{\mathop{A}\limits^{ \sim }}_{x}(t),{\mathop{A}\limits^{ \sim }}_{y}(t)]=\left[-\frac{{F}_{\omega }}{\omega }\sin (\omega t)-\frac{{F}_{2\omega }}{2\omega }\sin (2\omega t+\varphi ),\right.\\ \,\,\,\,\,\,\left.-\frac{{F}_{\omega }}{\omega }\cos (\omega t)+\frac{{F}_{2\omega }}{2\omega }\cos (2\omega t+\varphi )\right].\end{array}$$
(8)
When we perform the following switch of the circular field helicities \(\left[{\mathop{F}\limits^{ \sim }}_{x}\left(t\right),{\mathop{F}\limits^{ \sim }}_{y}\left(t\right)\right]\to \left[{\mathop{F}\limits^{ \sim }}_{x}\left(t\right),-{\mathop{F}\limits^{ \sim }}_{y}\left(t\right)\right]\),
$$\begin{array}{l}{\mathop{{\bf{F}}}\limits^{ \sim }}_{\circlearrowleft }(t)=[{F}_{\omega }\cos (\omega t)+{F}_{2\omega }\cos (2\omega t+\varphi ),\\ \,\,{F}_{\omega }\sin (\omega t)-{F}_{2\omega }\sin (2\omega t+\varphi )]\end{array}$$
(9)
and
$$\begin{array}{l}{\mathop{{\bf{A}}}\limits^{ \sim }}_{\circlearrowleft }(t)=\left[-\frac{{F}_{\omega }}{\omega }\sin (\omega t)-\frac{{F}_{2\omega }}{2\omega }\sin (2\omega t+\varphi ),\right.\\ \,\,\left.\frac{{F}_{\omega }}{\omega }\cos (\omega t)-\frac{{F}_{2\omega }}{2\omega }\cos (2\omega t+\varphi )\right].\end{array}$$
(10)
In this case, \({\mathop{{\bf{F}}}\limits^{ \sim }}_{\circlearrowright }\) and \({\mathop{{\bf{F}}}\limits^{ \sim }}_{\circlearrowleft }\) maintain the same spatial orientation, but the spatial orientation of the vector potential \({\mathop{{\bf{A}}}\limits^{ \sim }}_{\circlearrowright }\) is rotated by 180° with respect to that of \({\mathop{{\bf{A}}}\limits^{ \sim }}_{\circlearrowleft }\).
Regardless of the definition of the bicircular field that we use, it is the orientation of the vector potential relative to the crystal lattice (or Brillouin zone) that determines the band structure modification.
Laser-induced CNNN hopping
To understand the physical mechanism governing the band modification and the consequent valley polarization, we will use a two-orbital tight-binding model of gapped graphene, which captures the essential physics of the problem. The gapped-graphene lattice is the same as that of hBN, with one orbital per site. Those in sublattice A have on-site energy \({E}_{{\rm{A}}}\) whereas those in sublattice B have on-site energy \({E}_{{\rm{B}}}\). We will assume that, in the field-free state, next-nearest neighbour hoppings \({\gamma }^{({\rm{NNN}})}\) are negligible. However, upon interaction with the bicircular field, \({\gamma }^{({\rm{NNN}})}\) can be induced through virtual nearest neighbour hoppings \({\gamma }^{({\rm{NN}})}\). We will calculate this correction to lowest order in the field–matter interaction. To do so, note that a laser field modifies the nearest neighbour hopping between an orbital in sublattice A and an orbital in sublattice B according to the Peierls substitution,
$${\gamma }_{{\rm{AB}}}^{\left({\rm{NN}}\right)}\left(t\right)={\gamma }_{1}{e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{AB}}}\cdot {\bf{A}}\left(t\right)},$$
(11)
where \({\gamma }_{1}\) is the field-free nearest neighbour hopping, \({{\bf{r}}}_{{\rm{AB}}}\) is the distance from the site in sublattice B to the site in sublattice A and \({\bf{A}}\left(t\right)\) is the vector potential of the field. The laser-induced hopping term can be separated into one cycle-averaged term and one term that contains the dynamical corrections to this cycle average:
$${\gamma }_{{\rm{AB}}}^{\left({\rm{NN}}\right)}\left(t\right)={\gamma }_{1}\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{AB}}}\cdot {\bf{A}}\left(t\right)}\rangle +{V}_{{\rm{AB}}}\left(t\right),$$
(12)
where \(\langle \;\rangle \) means cycle-averaged. Therefore, we have
$${V}_{{\rm{AB}}}\left(t\right)={\gamma }_{1}{e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{AB}}}\cdot {\bf{A}}\left(t\right)}-{\gamma }_{1}\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{AB}}}\cdot {\bf{A}}\left(t\right)}\rangle .$$
(13)
This perturbation can lead to transitions between next-nearest neighbours (\({\rm{A}}\leftarrow {{\rm{A}}}^{{\prime} }\)) of the same sublattice through virtual nearest neighbours. To lowest order, the transition amplitude for such a second-order process reads (atomic units are used throughout)
$${{\mathcal{A}}}_{{\rm{A}}\leftarrow {A}^{{\prime} }}^{(2)}(t)=-\sum _{{\rm{B}}}{\int }_{{t}_{0}}^{t}{\rm{d}}{t}^{{\prime} }{\int }_{{t}_{0}}^{{t}^{{\prime} }}{\rm{d}}{t}^{{\prime\prime} }{e}^{{\rm{i}}({E}_{{\rm{AB}}}{t}^{{\prime} }+{E}_{{{\rm{BA}}}^{{\prime} }}{t}^{{\prime\prime} })}{V}_{{\rm{AB}}}({t}^{{\prime} }){V}_{{{\rm{BA}}}^{{\prime} }}({t}^{{\prime\prime} }),$$
(14)
where the energy difference \({E}_{{\rm{AB}}}={E}_{{\rm{A}}}-{E}_{{\rm{B}}}=\varDelta \), where \(\varDelta \) is the minimum bandgap energy (we take \({E}_{{\rm{A}}} > {E}_{{\rm{B}}}\)). The summation in principle runs along all possible intermediate states in sublattice B. However, note that only one nearest neighbour site can participate. For example, looking at Extended Data Fig. 1, we see that the transition from site A-11 to site A00 can happen only to first order through B01. Therefore, we can drop the summation. Using equation (13):
$$\begin{array}{l}{{\mathcal{A}}}_{{\rm{A}}\leftarrow {A}^{{\prime} }}^{(2)}(t)=-{\gamma }_{1}^{2}{\int }_{{t}_{0}}^{t}{\rm{d}}{t}^{{\prime} }{{\rm{e}}}^{i\varDelta {t}^{{\prime} }}\left[{{\rm{e}}}^{-{\rm{i}}{{\bf{r}}}_{{\rm{AB}}}\cdot {\bf{A}}({t}^{{\prime} })}-\langle {{\rm{e}}}^{-{\rm{i}}{{\bf{r}}}_{{\rm{AB}}}\cdot {\bf{A}}(t)}\rangle \right]\\ \,\,\,\,\,{\int }_{{t}_{0}}^{{t}^{{\prime} }}{\rm{d}}{t}^{{\prime\prime} }{{\rm{e}}}^{-{\rm{i}}\varDelta {t}^{{\prime\prime} }}\left[{{\rm{e}}}^{-{\rm{i}}{{\bf{r}}}_{{{\rm{BA}}}^{{\prime} }}\cdot {\bf{A}}({t}^{{\prime\prime} })}-\langle {{\rm{e}}}^{-{\rm{i}}{{\bf{r}}}_{{{\rm{BA}}}^{{\prime} }}\cdot {\bf{A}}(t)}\rangle \right].\end{array}$$
(15)
We perform a change of variables in the second integral, \(u=-{\rm{i}}\left[\varDelta {t}^{{\prime\prime} }+{{\bf{r}}}_{{{\rm{BA}}}^{{\prime} }}\cdot {\bf{A}}\left({t}^{{\prime\prime} }\right)\right]\), so that \({\rm{d}}u=-{\rm{i}}\left[\varDelta -{{\bf{r}}}_{{{\rm{BA}}}^{{\prime} }}\cdot {\bf{F}}\left({t}^{{\prime\prime} }\right)\right]\,{\rm{d}}{t}^{{\prime\prime} }\). For moderately strong fields and large gap materials, as in this case, we have that \(|{\bf{F}}\cdot {{\bf{r}}}_{{\rm{B}}{\rm{A}}}|\ll |\varDelta |\), so that we can neglect the second term in \({\rm{d}}u\) and write \({\rm{d}}u=-{\rm{i}}\varDelta \,{\rm{d}}{t}^{{\prime\prime} }\). In this way, the second integral is easily computed as
$$\begin{array}{c}{\int }_{u({t}_{0})}^{{u(t}^{{\prime} })}{\rm{d}}u\frac{{e}^{u}}{-{\rm{i}}\varDelta }-\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{{\rm{B}}{\rm{A}}}^{{\prime} }}\cdot {\bf{A}}(t)}\rangle {\int }_{{t}_{0}}^{{t}^{{\prime} }}{\rm{d}}{t}^{{\prime\prime} }{e}^{-{\rm{i}}\varDelta {t}^{{\prime\prime} }}\,=\,{\rm{i}}\frac{{e}^{-{\rm{i}}[\varDelta {t}^{{\prime} }+{{\bf{r}}}_{{{\rm{B}}{\rm{A}}}^{{\prime} }}\cdot {\bf{A}}({t}^{{\prime} })]}}{\varDelta }\\ -{\rm{i}}\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{{\rm{B}}{\rm{A}}}^{{\prime} }}\cdot {\bf{A}}(t)}\rangle \frac{{e}^{-{\rm{i}}\varDelta {t}^{{\prime} }}}{\varDelta }-{\rm{i}}\frac{{e}^{-{\rm{i}}\varDelta {t}_{0}}}{\varDelta }+{\rm{i}}\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{{\rm{B}}{\rm{A}}}^{{\prime} }}\cdot {\bf{A}}(t)}\rangle \frac{{e}^{-{\rm{i}}\varDelta {t}_{0}}}{\varDelta }\\ \,=\,\frac{{\rm{i}}}{\varDelta }{e}^{-{\rm{i}}\varDelta {t}^{{\prime} }}[{e}^{-{\rm{i}}{{\bf{r}}}_{{{\rm{B}}{\rm{A}}}^{{\prime} }}\cdot {\bf{A}}({t}^{{\prime} })}-\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{{\rm{B}}{\rm{A}}}^{{\prime} }}\cdot {\bf{A}}(t)}\rangle ]-\frac{{\rm{i}}}{\varDelta }{e}^{-{\rm{i}}\varDelta {t}_{0}}[1-\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{{\rm{B}}{\rm{A}}}^{{\prime} }}\cdot {\bf{A}}(t)}\rangle ].\end{array}$$
(16)
Substituting equation (16) into equation (15):
$$\begin{array}{l}{{\mathcal{A}}}_{{\rm{A}}\leftarrow {A}^{{\prime} }}^{(2)}(t)=-\frac{{\rm{i}}{\gamma }_{1}^{2}}{\varDelta }{\int }_{{t}_{0}}^{t}{\rm{d}}{t}^{{\prime} }\{[{e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{A}}{\rm{B}}}\cdot {\bf{A}}({t}^{{\prime} })}-\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{A}}{\rm{B}}}\cdot {\bf{A}}(t)}\rangle ]\\ \,\,\times \,[{e}^{-{\rm{i}}{{\bf{r}}}_{{{\rm{B}}{\rm{A}}}^{{\prime} }}\cdot {\bf{A}}({t}^{{\prime} })}-\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{{\rm{B}}{\rm{A}}}^{{\prime} }}\cdot {\bf{A}}(t)}\rangle ]-{e}^{-{\rm{i}}\varDelta ({t}_{0}-{t}^{{\prime} })}\\ \,\,\times \,[1-\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{{\rm{B}}{\rm{A}}}^{{\prime} }}\cdot {\bf{A}}(t)}\rangle ][{e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{A}}{\rm{B}}}\cdot {\bf{A}}({t}^{{\prime} })}-\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{A}}{\rm{B}}}\cdot {\bf{A}}(t)}\rangle ]\}.\end{array}$$
(17)
The above expression is for the lowest-order transition amplitude between two sites in sublattice A. From it, we can identify the transition matrix element between the next-nearest neighbour sites A′ and A (\({\gamma }_{{{\rm{AA}}}^{{\prime} }}^{\left({\rm{NNN}}\right)}\)) by realizing that the transition amplitude is defined as
$${{\mathcal{A}}}_{i\leftarrow k}\left(t\right)=-{\rm{i}}{\int }_{{t}_{0}}^{t}{\rm{d}}{t}^{{\prime} }{e}^{{\rm{i}}\left({E}_{i}-{E}_{k}\right){t}^{{\prime} }}{\gamma }_{ik}\left({t}^{{\prime} }\right).$$
(18)
Hence,
$$\begin{array}{c}{\gamma }_{{{\rm{A}}{\rm{A}}}^{{\prime} }}^{({\rm{N}}{\rm{N}}{\rm{N}})}=\frac{{\gamma }_{1}^{2}}{\varDelta }\{[{e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{A}}{\rm{B}}}\cdot {\bf{A}}({t}^{{\prime} })}-\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{A}}{\rm{B}}}\cdot {\bf{A}}(t)}\rangle ][{e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{B}}{\rm{A}}}^{{\prime} }\cdot {\bf{A}}({t}^{{\prime} })}-\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{{\rm{B}}{\rm{A}}}^{{\prime} }}\cdot {\bf{A}}(t)}\rangle ]\\ \,\,-{e}^{-i\varDelta ({t}_{0}-{t}^{{\prime} })}[1-\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{{\rm{B}}{\rm{A}}}^{{\prime} }}\cdot {\bf{A}}(t)}\rangle ][{e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{A}}{\rm{B}}}\cdot {\bf{A}}({t}^{{\prime} })}-\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{A}}{\rm{B}}}\cdot {\bf{A}}(t)}\rangle ]\}\\ \,\,=\,\frac{{\gamma }_{1}^{2}}{\varDelta }\{{e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{A}}{\rm{A}}}^{{\prime} }\cdot {\bf{A}}({t}^{{\prime} })}-{e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{A}}{\rm{B}}}\cdot {\bf{A}}({t}^{{\prime} })}\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{B}}{\rm{A}}}^{{\prime} }\cdot {\bf{A}}(t)}\rangle \\ \,\,-\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{A}}{\rm{B}}}\cdot {\bf{A}}(t)}\rangle {e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{B}}{\rm{A}}}^{{\prime} }\cdot {\bf{A}}({t}^{{\prime} })}+\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{A}}{\rm{B}}}\cdot {\bf{A}}(t)}\rangle \langle {e}^{-{\rm{i}}{{\bf{r}}}_{{{\rm{B}}{\rm{A}}}^{{\prime} }}\cdot {\bf{A}}(t)}\rangle \\ \,\,-{e}^{-{\rm{i}}\varDelta ({t}_{0}-{t}^{{\prime} })}[{e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{A}}{\rm{B}}}\cdot {\bf{A}}({t}^{{\prime} })}-\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{A}}{\rm{B}}}\cdot {\bf{A}}(t)}\rangle \\ \,\,-{e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{A}}{\rm{B}}}\cdot {\bf{A}}({t}^{{\prime} })}\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{{\rm{B}}{\rm{A}}}^{{\prime} }}\cdot {\bf{A}}(t)}\rangle +\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{A}}{\rm{B}}}\cdot {\bf{A}}({t}^{{\prime} })}\rangle \langle {e}^{-{\rm{i}}{{\bf{r}}}_{{{\rm{B}}{\rm{A}}}^{{\prime} }}\cdot {\bf{A}}(t)}\rangle ]\}.\end{array}$$
(19)
Averaging now over one cycle,
$$\langle {\gamma }_{{{\rm{AA}}}^{{\prime} }}^{\left({\rm{NNN}}\right)}\rangle =\frac{{\gamma }_{1}^{2}}{\varDelta }\left\{\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{{\rm{AA}}}^{{\prime} }}\cdot {\bf{A}}\left(t\right)}\rangle -\langle {e}^{-{\rm{i}}{{\bf{r}}}_{{\rm{AB}}}\cdot {\bf{A}}\left(t\right)}\rangle \langle {e}^{-{\rm{i}}{{\bf{r}}}_{{{\rm{BA}}}^{{\prime} }}\cdot {\bf{A}}\left(t\right)}\rangle \right\}.$$
(20)
If the transition is between two next-nearest neighbours of the other sublattice (B), then we need to substitute \(\varDelta \to -\varDelta \), and
$$\langle {\gamma }_{{{\rm{BB}}}^{{\prime} }}^{\left({\rm{NNN}}\right)}\rangle =-\langle {\gamma }_{{{\rm{AA}}}^{{\prime} }}^{\left({\rm{NNN}}\right)}\rangle .$$
(21)
Extended Data Fig. 2 shows some representative examples of \(\langle {\gamma }_{{{\rm{BB}}}^{{\prime} }}^{\left({\rm{NNN}}\right)}\rangle \) and \(\langle {\gamma }_{{{\rm{AA}}}^{{\prime} }}^{\left({\rm{NNN}}\right)}\rangle \) for different vector potentials. Note the following. First, the cycle-averaged, laser-induced, next-nearest neighbour, hoppings are complex whenever the vector potential is not pointing towards the M direction. The imaginary component is a maximum when the vector potential points along K or K′, and it switches sign between these two orientations. Second, the hopping depends only on the orientation of the vector potential and not its sense of rotation (compare Extended Data Figs. 2b and 2d). The band structure is modified by rotating the vector potential, which can be achieved on a subcycle timescale by changing the two-colour phase delay \(\varphi \). Finally, also note that even when the vector potential is pointing towards the point M, the bandgap is reduced relative to the field-free case due to a modification of the hopping. In this case, however, as the next-nearest neighbour hopping is real, the gap is reduced equally in both valleys.
Numerical calculations
We performed time-dependent simulations of a tight model of gapped graphene using the code described in ref. 44. The field-free tight-binding parameters are \({r}_{0}=2.73\) a.u., \(\varDelta ={\varDelta }_{{\rm{A}}}-{\varDelta }_{{\rm{B}}}\), where \({\varDelta }_{{\rm{A}}}=5.9/2\) eV and \({\varDelta }_{{\rm{B}}}=-5.9/2\) eV, and \({\gamma }_{1}=0.089\) a.u. The atomic distance and first neighbour hopping are taken to be like those of graphene. Next-nearest neighbour hoppings and higher were neglected. Owing to the uncertainty in the experimental intensity, we simulated several ratios and intensities of the bicircular (trefoil) field. The fields were simulated using a Gaussian envelope with 30 fs of full-width at half-maximum for both fields, which matches the estimated duration of the fields in the experiment. The time-dependent propagation was converged on a Monkhorst–Pack k grid of 300 × 300 points and a time step of 0.1 a.u. The dephasing time was set to 6 fs, but different values did not change our results.
Extended Data Fig. 2 shows the results for two different helicities of the bicircular field. In this case, the field parameters are \({F}_{\omega }=2{F}_{2\omega }=0.0085\) a.u., which correspond to an intensity in the crystal of \({I}_{\omega }=2.5\) TW cm−2 and \({I}_{2\omega }=0.63\) TW cm−2. First, note that valley polarization does not change between Extended Data Fig. 2c and Extended Data Fig. 2f, despite having switched the field helicity. Second, the valley polarization switches in Extended Data Fig. 2c and Extended Data Fig. 2f upon a 60° rotation of the bicircular field. Third, there is a not a perfect interchange of the K and K′ valley populations in Extended Data Fig. 2f,i. This is because the fields (or the associated vector potential) that give rise to the populations in Extended Data Fig. 2f,i are not time-reversal partners. However, the fields in Extended Data Fig. 2c and Extended Data Fig. 2i do switch exactly the K and K′ populations, since the fields in this case are time-reversal partners. Yet, it is clear that, for fixed helicity, the orientation of the field relative to the lattice controls the valley polarization, even if the switching is not fully symmetric. This effect signals band modification by the strong bicircular field and allows for subcycle control.
k-resolved populations
Extended Data Fig. 3 shows the k-resolved populations after the interaction with the bicircular field for different intensities using a fixed ω to 2ω ratio of 2:1 in intensity. We found similar results for intensity ratios of 1:1, 4:1 and 6:1. The polarized valley always corresponds to that in which the model predicts that the bandgap is reduced. Also, the valley polarization increases as a function of intensity, in accordance with the model prediction that the effective bandgap decreases as the intensity is increased.
Probing the valley polarization
To transfer these valley populations into an optical degree of freedom that can be measured in the experiment, we used a linearly polarized probe pulse after the bicircular pulse which, as explained in the main text, allowed us to map them into the helicity of its nonlinear harmonic response (H3 in this case). As in the experiment, we used a probe light field of the same ω frequency as the fundamental field in the bicircular pulse. The field strength of the probe field was \({F}_{\text{probe}}=0.01\) a.u., it was polarized along the Γ–M direction (x direction in the figures) and its duration was 30 fs of full-width at half-maximum. We tested different probe field strengths, but these did not affect our results.
Extended Data Fig. 4 shows the results of the polarimetry analysis. The two helical components that the two photodiodes measure are defined as
$$\begin{array}{l}{\widehat{e}}_{{\rm{l}}}\,=\,\frac{1}{\sqrt{2}}({\widehat{e}}_{x}+{\rm{i}}{\widehat{e}}_{y}),\\ {\widehat{e}}_{{\rm{r}}}\,=\,\frac{1}{\sqrt{2}}({\widehat{e}}_{x}-{\rm{i}}{\widehat{e}}_{y}).\end{array}$$
(22)
Therefore,
$$\begin{array}{l}{\widehat{e}}_{{\rm{l}}}\,=\,{\left|\frac{{A}_{x}}{\sqrt{2}}{e}^{{\rm{i}}{\phi }_{x}}+\frac{{A}_{y}}{\sqrt{2}}{e}^{{\rm{i}}({\phi }_{y}+{\rm{\pi }}/2)}\right|}^{2},\\ {\widehat{e}}_{{\rm{r}}}\,=\,{\left|\frac{{A}_{x}}{\sqrt{2}}{e}^{{\rm{i}}{\phi }_{x}}+\frac{{A}_{y}}{\sqrt{2}}{e}^{{\rm{i}}({\phi }_{y}-{\rm{\pi }}/2)}\right|}^{2},\end{array}$$
(23)
where Ax and Ay (ϕx and ϕy) are the Fourier amplitudes (phases) of the harmonic of interest of the probe (H3 in this case), along directions x and y, respectively. We can rewrite the above as
$$\begin{array}{l}{\widehat{e}}_{{\rm{l}}}\,=\,\frac{1}{2}| {A}_{x}| +\frac{1}{2}| {A}_{y}| +| {A}_{x}| | {A}_{y}| \cos ({\phi }_{x}-{\phi }_{y}-{\rm{\pi }}/2),\\ {\widehat{e}}_{{\rm{r}}}\,=\,\frac{1}{2}| {A}_{x}| +\frac{1}{2}| {A}_{y}| +| {A}_{x}| | {A}_{y}| \cos ({\phi }_{x}-{\phi }_{y}+{\rm{\pi }}/2).\end{array}$$
(24)
The phase ϕy changes by π as the valley polarization changes sign since the y component of the current comes from the anomalous contribution and, thus, is proportional to the Berry curvature, which has the same magnitude but opposite sign in both valleys. As the valley polarization changes upon rotation of the bicircular field, the interference (cosine) term in the expression above oscillates sinusoidally, switching sign with a switch in valley polarization. Additionally, the interference term of the two helicities oscillates out of phase because of the factor ±π/2. Extended Data Fig. 4 shows the interference term of the two helicity components, which indeed oscillate sinusoidally and out of phase. As expected, these interference terms are a maximum or a minimum when there is maximum valley polarization and zero when there is none. The asymmetry of these two interference terms completely characterizes the valley polarization.
Yet, the signal observed in the experiment also includes the amplitude term in the equation above, \(\frac{1}{2}|{A}_{x}|+\frac{1}{2}|{A}_{y}|\), which also oscillates. Its oscillation comes from unequal population injection as a function of rotation and other nonlinear effects occurring during the harmonic generation. However, this term is the same for both helicities, and thus, it is merely a background introducing higher-order Fourier components into the oscillation. We plot this term in Extended Data Fig. 4.
The helicity signal is then the sum of the amplitude term, which is common to both helicities, and the interference term, which is different for each helicity and which contains the information on the valley polarization. The total left and right signals give, respectively, the red and green curves in Extended Data Fig. 4 and also in the main text, which is the curve to be compared with the experiment.
Importantly, regardless of the value of the amplitude term, since it is a background common to both helicities, we can remove its influence. For this, we Fourier-filtered the oscillation to extract only the 120° periodic oscillation. In this way, the asymmetry of the two helicity signals after Fourier filtering characterizes the valley polarization, with the change of sign indicating valley switching.
Note that the helicity-resolved H3 probe signal as a function of the pump trefoil rotation is essentially the same regardless of whether the probe is polarized along Γ–K or Γ–M (Extended Data Fig. 5).
Experimental details
Laser system
The laser pulses used for the experiments were from a mid-infrared laser system based on optical parametric chirped-pulse amplification (OPCPA). The details are in ref. 45. In brief, the laser for the mid-infrared OPCPA was from a 1 µm Innoslab Yb laser system. The pulses were used to generate the seed and the pump for the OPCPA. The OPCPA system produced pulses at 2.03 µm centre wavelength with a 100 kHz repetition rate. During the experiment, the average power obtained from the laser was 5.5–6 W with a pulse duration of approximately 25 fs. The pulse duration was carefully characterized using a frequency-resolved optical gating46 with a set-up developed in-house. The laser pulses were characterized after they had passed through the interferometer. To accurately determine the temporal characteristics of the pulses used in the experiment, they were picked up right before they entered the reflective objective of the microscope arrangement. A second-harmonic generation (SHG) frequency-resolved optical gating was used to characterize the 1 µm pulses in a beta-barium borate crystal, whereas the 2 µm pulses were characterized by a surface third-harmonic generation frequency-resolved optical gating on glass (NBK7), both in a non-collinear geometry. The results of these measurements are shown below.
Interferometer
A large part of the experimental set-up was a three-arm Mach–Zehnder interferometer, as shown in Extended Data Fig. 6. Various aspects of the interferometer are described in detail in the following sections.
Pump
The 2 µm wavelength pulses from the laser source were fed directly into the three-arm interferometer. First, the pulses were split in two by a 50–50 beam splitter (red beam from right to left in Extended Data Fig. 6). The reflected part proceeded to the pump and the transmitted one to the probe arm. The pump beam was sent through a 1.5-mm-thick lithium niobate (LiNbO3) crystal (SHG) cut at 45.6°. This interaction produced a second pulse at half the wavelength (frequency of 2ω) of the fundamental (frequency of ω) by sum frequency generation (SHG)47. LiNbO3 is suited well for a 2 µm fundamental wavelength, given its wide transparency range in the mid-infrared of up to 3.5 µm. The phase-matching conditions (type 1) led to the second-harmonic component being cross-polarized with respect to the fundamental pulse47. The co-propagating two-colour pulses (ω and 2ω) were separated using a custom-designed dichroic beam splitter. The ω arm passed through a variable neutral-density filter, a retro-reflection stage and CaF2 as used in the probe arm, but without the piezo motion rails, as the length of this arm was fixed with respect to the other arms. Additionally, a wire grid polarizer was placed in its path to clean its polarization state before it combined with the 2ω arm at another dichroic beam splitter. The 2ω arm was reflected off the first dichroic beam splitter at 90° and traversed the exact same arrangement as the ω arm. The only difference was a different pair of custom chirped mirrors for 2ω on the delay stage above a closed-loop stick–slip piezo nano-positioning rail with a positioning resolution of 1 nm, a different amount of dispersion-compensating material and a perpendicular axis set on the polarizer. The pump arm (co-propagating ω and 2ω) was further sent through a long-pass spectral filter (F1) to block the parametric optical signals generated at the harmonics of the two-colour pump in the LiNbO3 crystal. Finally, the two-colour pump (ω and 2ω) pulses passed through the fused-silica wedged pair to recombine with the probe pulses (ω), which were reflected off it.
After the three-arm interferometer, the laser beam was expanded using a reflective telescope arrangement to roughly match the input diameter of the tight-focusing reflective objective and adjust the beam divergence. However, right before the focusing objective, there was a broadband quarter-wave plate (λ/4) with its optical axis at 45° with respect to the cross-polarized axis of the two-colour pump. This wave plate arrangement transformed the linear, cross-polarized, pump pulses into circularly polarized, counter-rotating pulses that look like a trefoil or a three-leaf pulse in the X–Y plane, as required in the experiment. The quarter-wave plate was intentionally placed right before the focusing element to prevent any changes in ellipticity due to a phase shift between the s- and p-polarized components upon reflection, especially from beam-folding mirrors and other coated optical elements. For the ω and 2ω pump components, the pulse durations measured just before the focusing objective were about 26 and 48 fs, respectively. Extended Data Fig. 7 depicts the spectral and temporal characteristics of the ω and 2ω pump components. A similar kind of waveform synthesizer was also used earlier in attosecond-controlled strong-field experiments by the group36. The data shown in the manuscript were obtained with the pump intensity in the range 4–7 TW cm−2. For the intensity scaling measurements, the overall power of the pump beam was changed with a variable neutral-density filter keeping the relative power ratio between the components intact.
Probe
In the probe arm, the pulses went through a variable neutral-density filter followed by a delay stage, which, along with a pair of silver retro-reflecting mirrors, hosted two customized chirped mirrors for simultaneous positive dispersion and spectral filtering. Spectral filtering is crucial given the existence of weak optical signals at lower wavelengths arising from the parametric processes in various crystals in the laser system. The delay stage was mounted on a closed-loop stick–slip piezo rail (Smaract SLC-24 series) with a positioning resolution of 1 nm. Further down, the pulses went through a defined thickness of material (CaF2) to compensate for excess positive dispersion. Additionally, another long-pass spectral filter was placed in the beam path to further suppress the unwanted optical signals at lower wavelengths. Finally, the probe pulses went through a half-wave (λ/2) plate and a quarter-wave (λ/4) plate, after which they were reflected off the wedge plate and recombined with the pump beam. This wave-plate combination allowed us to control the polarization state of the probe pulses, and the mechanism is described in greater detail later. This wedge pair arrangement not only acted as a beam recombination element but also as a power attenuator for the probe pulses, as only 4% of the power was reflected. After being recombined, the probe beam followed the collinear path with the pump beam. Just before the focusing objective, a pulse duration of about 26 fs was achieved for the probe pulses.
As shown in Extended Data Fig. 6, a λ/4 plate was the last optical element before the pulses entered the reflective objective. This led to a major problem in which the third arm (or the linearly polarized probe as in this experiment) cannot remain linear once it has passed through the λ/4 plate unless it is along the optical axis at 45°. Also, a linear polarization launched at 45° with respect to the s- or p-polarized states would lose its linear contrast, as it would become elliptical on acquiring a different phase shift in the s and p components on every reflective optic in its beam path. To have full flexibility over the polarization state of the probe pulses after the λ/4 plate, a scheme was implemented such that a combination of λ/2 and λ/4 wave plates were additionally placed in the probe arm, as depicted in Extended Data Fig. 6. Intuitively, one can think of these additional plates as inducing perfectly opposing elliptical polarization, which cancels out in the final λ/4 before the reflective objective to produce linearly polarized light with a high extinction ratio. This scheme was numerically tested using a Jones matrix approach. It was observed that any arbitrary shifts in phase between the s and p components in the beam path between the two λ/4 wave plates can be compensated. However, changes in magnitude between the s and p components lead to a deviation from linearity and cannot be compensated for by this scheme.
A movable silver mirror was placed at 45° right before the reflective objective intercepting the probe pulses to optimize and characterize the polarization extinction ratio or linearity. The pulses were then guided to an InGaAs photodiode with a polarizer attached to it at a fixed angle. The fixed angle was such that the probe polarization was aligned along s or p to prevent any additional ellipticity induced by the intercepting mirror, which would not be present otherwise during the experiment.
The data shown in the main manuscript were obtained at a pump–probe delay of about 60–110 fs.
Trefoil pump rotation
When the bicircular ω and 2ω components of the pump were combined in the interferometer such that the E field ratio at the focus was 2:1, the coherent sum of their electric fields in the X–Y plane (the plane perpendicular to propagation) was transformed to that of a trefoil waveform. The rotation of the trefoil waveform was controlled by the subcycle phase delay between the ω and 2ω components. Experimentally, this was achieved by introducing an optical path difference between the ω and 2ω arms. When the central wavelength (λ0) of the ω arm was 2 µm, a rotation of 360° was induced by delaying the piezo stage by 3 µm. This information was used to convert the stage delay to angular rotation, which was recorded as the experimental data.
Interferometric stability
During the experiment, it was critical that the angle of the three-leaf or trefoil pump remained stable, as this was directly linked to the delay between the two colours in the pump arm. To characterize the delay stability, an additional temperature-stabilized continuous-wave diode laser (Thorlabs L785P090 with LDMT9) was sent through the two pump arms to interferometrically measure its path difference over time.
Using the above-mentioned scheme, the interferometric stability was found to be highest around 2 h after switching on the driving laser system. Additionally, tests were carried out to measure the stability when the laser was going through the interferometer and the piezo delay stage being scanned, as during the experiment, as illustrated in Extended Data Fig. 8. Over a period of 10 min, the standard deviation of the position generated by the closed-loop piezo stages from the position extracted from the continuous-wave interferometer was close to 38.1 nm, which roughly translates to about 4.6° in rotation error of the ω and 2ω bicircular trefoil structure. A similar scheme with active stabilization was used earlier and achieved few tens of attosecond interferometric stability48.
Microscope
The core segment of the microscope (Extended Data Fig. 9) applied tight focusing of laser pulses from all three interferometer arms, followed by a collection and collimation configuration before the light proceeded to the detection apparatus. A 500-µm-thick fused-silica substrate with the hBN monolayer on its front surface was placed between these components.
The focusing was achieved without adding any dispersion or chromatic aberration by using a reflective objective with a numerical aperture of 0.4. The beam profile in the interaction region was characterized by placing a CMOS sensor in the focal plane. A cross section of the resulting two-dimensional signal on the sensor is shown in Extended Data Fig. 10 for all three interferometer arms. The beam waists extracted from the fit were about 8.0, 10.0 and 9.7 µm for the 2ω pump, ω pump and ω probe arms, respectively. The beam was largely circular at the focus, and the residual ellipticity was below 5% for each arm. The small focus spot sizes not only allowed us to spatially restrict the laser interaction to a single hBN monocrystalline grain but also gave a Rayleigh length (201 µm for 2ω and 157 µm for ω) that was significantly less than the thickness of the fused-silica substrate (Extended Data Fig. 9). By observing the surface-enhanced perturbative harmonics, we could individually optimize the substrate position (back or front of the hBN-coated surfaces) at the focus.
Once the fused-silica substrate was at the focus, the diverging beam along with other nonlinear signals was recollimated using a long-working-distance transmissive objective with a numerical aperture of 0.45. Using a closely matched numerical aperture allowed the divergent beams to be entirely collected.
Further down, the laser pulses were polarization resolved using a quartz Wollaston prism. The prism introduced an angular separation of 10° between the p- and s-polarized components, which were then loosely focused by a CaF2 lens onto the two separate large-area photodiodes. A λ/4 wave plate was placed between the transmission objective and the Wollaston prism with its optical axis at 45° with respect to the prism to convert the photodiodes into helicity detectors.
Spatio-temporal overlap
As described earlier, the implementation of a three-arm interferometer leads to the laser pulses travelling through three different non-collinear optical paths. Upon recombination, they require careful alignment such that they spatially and temporally overlap in the interaction region. To ensure this, the second-order and third-order nonlinearity in fused silica is utilized such that a two- and three-photon transition of ω and 2ω leads to the production of 3ω and 4ω photons, respectively. To utilize this effect, the 2ω pump is chosen as a common reference and the other pulses are aligned onto this particular arm. An initial spatial alignment of all three arms using irises is enough to ensure a coarse spatial overlap between the pulses at the focus on the fused-silica substrate. Soon after, the 2ω pump arm is delayed with respect to the ω pump arm while observing the 4ω light with an appropriate band-pass filter before the photodiodes. As the pump arms are counter-rotating bicircularly in their polarization state with respect to each other, the 3ω channel is disallowed. After the 2ω pump delay is fixed, its delay with respect to the probe arm is determined by scanning the probe delay stage and observing the 3ω or 4ω channels, as both these transitions are allowed. After the temporal alignment, the spatial overlap is checked and optimized based on the respective parametric signals.
Lock-in polarization detection
In this experiment, the signal of interest is the 3ω harmonic generated by the probe and modulated by the pump. However, several sources of light are co-propagating along with the signal and impinging on the diodes. Hence, a simple amplified detection of the diode current is not enough to distinguish the relevant 3ω signal (λ ≈ 667 nm). Signals that are spectrally different are already filtered using high- or low-pass spectral filters, which block the fundamental 2ω, components above λ = 750 nm and other components below λ = 600 nm. However, 3ω light can be produced by the pump (ω) and probe (ω) pulses individually by the third-harmonic upconversion process and by the combination of pump (2ω) and probe (ω) through a two-photon process on the silica substrate, thereby further polluting the experimental signal of interest.
We overcome this major problem by modulating either the ω arm or the 2ω arm of the interferometer using a mechanical chopper wheel with a 50% duty cycle at a known frequency. Hence, the 3ω signal of interest, which was probed by a modulated bicircular pump, was also modulated at the chopping frequency. To chop the ω arm, an unwanted 3ω signal, generated by the third-harmonic upconversion of ω pump, falls entirely on a single diode (D1) because its helicity remains the same as that of the ω pump. Therefore, the other diode (D2) remains background-free. Moreover, when the 2ω arm is chopped, an unwanted 3ω signal, generated through a pump (2ω) and probe (ω) two-photon process, falls on the D2 diode, as its helicity is the same as that of the 2ω signal. In that case, the D1 diode remains background-free. Therefore, a combination of two independent measurements while keeping all the parameters unchanged, first chopping the ω arm and subsequently the 2ω arm, can fulfil the requirement for detecting both the helicities of the 3ω signal in a background-free condition. Note that the continuous service of the continuous-wave reference interferometer ensures the interferometric stability between these two measurements. This modulated signal impinging on the photodiodes is amplified by two high-gain, low-noise, trans-impedance amplifiers placed right behind the respective photodiodes to minimize the noise picked up by the cables. The amplified voltage signals are then fed into two time-synchronized lock-in amplifiers (Zurich Instruments).
Extended Data Table 1 lists the signals detected along with their origin. It also depicts the method by which the background signals are rejected while keeping the relevant signal detectable.
hBN monolayer sample
The hBN monolayer samples used in this experiment were obtained from a commercial supplier (ACS Materials). These were grown using chemical vapour deposition and subsequently transferred onto a 500-µm-thick, large-area (5 cm diameter), fused-silica substrate. Although the monolayer patches, of characteristic size of about 10 µm as typical for samples grown by chemical vapour deposition, had their large areas covered, their orientations were random. These individual hBN patches were characterized in situ by microscopy.
SHG polarization-resolved microscopy
To characterize the quality of the hBN patches individually, we used in situ polarization-resolved SHG microscopy. This was repeated before every experimental run.
SHG polarization-resolved microscopy was done with the 1 µm wavelength arm of the set-up. The perturbative second harmonic was detected by a photodiode or spectrometer after being collected by the transmission objective and passed through band-pass filters and eventually a polarizer. The SHG is not only the strongest component but being an even harmonic, it is a pure indicator of broken inversion symmetry in a crystal. Hence, the emission was also fully absent from the fused-silica surface, leaving only odd layers of hBN (1L, 3L, 5L, …, (2n + 1)L), which break inversion symmetry to produce SHG. As n increases, the magnitude of SHG drops49,50. Hence, a spatial scan over the substrate with the SHG strength as an optimization parameter can locate the thinner odd layers of hBN.
Another important advantage of this technique is that it can be performed in situ during the final band engineering experiment to make sure that all three probe and bicircular pump pulses also excite the same hBN grain as independently observed by this characterization technique.
Further, not only can this technique determine the inversion-symmetry breaking but it can also detect the orientation of the hBN crystals by observing the polarization response of the material. The two unique light-induced oscillation directions (armchair and zigzag) make the outgoing SHG polarization rotate in the opposite direction and at twice the rate, as a function of the rotation of the incoming (pump) light polarization. This has also been widely investigated and demonstrated in recent works49,50. Two results from two different positions on the hBN sample used are displayed in Extended Data Fig. 11. The characteristic four-lobe polarization-dependent signal implies that the interaction region of the laser irradiation is small enough to fit into the monocrystalline hBN patch. Thus, the single crystalline orientation is dominant and is not averaged out over several grains with different crystalline orientations. Further, Extended Data Fig. 11a,b shows that the four-lobe structure exhibits an angular offset with the four-lobe structure obtained from a neighbouring hBN patch. This manifests from the hBN crystalline patches, which are oriented at arbitrary angles at different laser excitation sites. Usually, a six-lobe feature in the SHG signals can be used to indicate the orientation of hBN (ref. 49). This difference is attributed to a lab-to-incident polarization frame conversion. In the present case, the SHG polarization microscopy was done in transmission mode while keeping the sample fixed for a given hBN patch. The sensitivity of this technique to the presence of hBN on the fused-silica substrate allows it to be used as a marker for laser-induced damage. This confirms that the suitable intensity range before damage occurs is close to 1013 W cm−2 .
