Full 2π tunable phase modulation using avoided crossing of resonances

Active metasurfaces have been proposed as one attractive means of achieving high-resolution spatiotemporal control of optical wavefronts, having applications such as LIDAR and dynamic holography. However, achieving full, dynamic phase control has been elusive in metasurfaces. In this paper, we unveil an electrically tunable metasurface design strategy that operates near the avoided crossing of two resonances, one a spectrally narrow, over-coupled resonance and the other with a high resonance frequency tunability. This strategy displays an unprecedented upper limit of 4π range of dynamic phase modulation with no significant variations in optical amplitude, by enhancing the phase tunability through utilizing two coupled resonances. A proof-of-concept metasurface is justified analytically and verified numerically in an experimentally accessible platform using quasi-bound states in the continuum and graphene plasmon resonances, with results showing a 3π phase modulation capacity with a uniform reflection amplitude of ~0.65.


Supplementary Note 1. Radiative and dissipative loss rates of the metasurface
Following the formalism of the TCMT (temporal coupled mode theory) section below, we can extract the values of the qBIC (quasi-bound states in the continuum) radiative loss rate and its dissipative loss rate respectively, from full wave electromagnetic simulation results of a setup with a graphene sheet deposited through the entire unit cell underneath the Si pillars. In this setup graphene plasmons will be off-resonance and too weak to play a role, therefore satisfying the condition of a singular resonance (the qBIC) system and simplifying the equations down to Supplementary Eq. 1, discussed in Supplementary Note 3. The extracted radiative loss rates and the dissipative loss rates of the qBIC depending on the metasurface structural changes are shown in Supplementary Fig. 1. We can note that the structural changes do not largely affect dissipative loss rates, whereas they heavily influence radiative loss rates.
Supplementary Figure 1. The dependence of the qBIC's radiative and dissipative loss rates with respect to structural parameters. The dependence with respect to a, the Si pillar post width b, the post height c, the and d, the substrate thickness. The dissipative loss rates are not largely affected by the changes in the structural geometry, whereas the radiative coupling is heavily influenced by the structural changes due to the corresponding changes in the strengths of symmetry-breaking-induced dipole moments.

Radiative loss rates
The radiative loss is largely affected because it is dependent on the strength of the dipole moment (allowing coupling to free-space radiation) that is formed due to the dimerization (symmetrybreaking) of the metasurface, and this is greatly influenced by the structural changes. For example, raising the height of the Si pillars would increase the amount of charge accumulation that leads to a dipole formation, resulting in a larger radiative coupling. Increasing post width will decrease the central gap, leading to a shorter dipole length/smaller dipole moment and therefore a lower radiative loss rate. The dimerizing perturbation is defined to express the percentage of the shift of the center of the Si pillars in the direction of dimerization, and the magnitude of the shift is defined by /2. The dimerizing perturbation affects the degree of the perturbation (symmetry breaking/dimerization) and a larger delta would correspond to a larger dipole moment, resulting in a larger radiative coupling. This relation is shown to be linear in the limit of smaller perturbation and non-linear as the perturbation gets larger (33). The substrate thickness changes the degree of the Fabry-Perot effect between the back-reflector and the Si pillars, and the radiative coupling draws a rough sinusoidal shape as the dipole moment is affected by the constructive/destructive interference with the light reflected from the back reflector.

Dissipative loss rates
The dissipative loss rate, on the other hand, is not majorly affected by structural changes but instead greatly on the changes in the graphene Fermi level. This is shown in Supplementary Fig. 2a, where the initial hump of the dissipative loss rate near 0 eV is due to the interband transition losses and the steady increase from 0.2 eV is due to the changes in the intraband transition losses. The resonance frequency of the qBIC is roughly 41 THz, which would mean that at this frequency interband transition losses would be cut off from around roughly 0.08-0.09 eV due to Pauli blocking and the corresponding change in the graphene's optical conductivity (27). This explains the drop in the dissipative loss rate from 0 eV to 0.2 eV. Then the intraband scattering in graphene slowly increases as the Fermi level increases due to the rise in the carrier density and an increase in the optical conductivity (27), which explains the steady increase in the dissipative loss rate.
We can employ graphene plasmons as well as the qBIC by depositing graphene ribbons between the Si pillars, as opposed to a whole graphene sheet underneath the Si pillars throughout the unit cell. In this case, the coupling between the qBIC and the graphene plasmon affects the radiative and the dissipative loss rates and therefore the extracted loss rate values from using a single resonance TCMT model becomes unreliable. Using the two-resonance TCMT model explained below would entail too many parameters at once and the extraction of the loss rate values would be unviable as well. This is shown in Supplementary Fig. 2b, where the high coupling zone region from 0.5-0.8 eV is denoted by a grey box. The loss rate values for the qBIC shows a similar trend to that of Supplementary Fig. 2a, with slight differences in the high Fermi level range due to the presence of graphene plasmons occurring in the side graphene ribbons. This can be seen in Fig. 3e of the main paper, where GP2 (the secondary graphene plasmon) approaches the qBIC line and further broadens the linewidth. This would technically introduce errors to our loss rate extraction method based on the single resonance TCMT model explained below. The loss rates for the graphene plasmon are shown as well, with the dissipative loss rate nicely following the DC scattering rate 1/ = 2 / where is the electron charge, is the mobility currently set at 1000 cm 2 /(V⋅s), and ≈ /300 is the Fermi velocity (42). The graphene plasmon is a dipole radiator as can be seen Fig. 4c of the main paper, and radiation from a dipole radiator increases as the frequency of the radiation increases (43). Because the resonance frequency of the graphene plasmon drastically increases as the Fermi level is increased (Fig. 3e of the main paper), its radiative loss rate also increases until the high coupling zone.
Supplementary Figure 2. The dependence of the radiative and dissipative loss rates with respect to the graphene Fermi level. a, The qBIC's radiative loss rate is not largely affected by the changes in the graphene fermi level, whereas it heavily influences the dissipative loss rate due to the interband/intraband scattering in graphene. b, Dissipative loss rates for both the qBIC and graphene plasmons (employed through graphene ribbons as opposed to a graphene sheet in a). The qBIC's radiative and dissipative loss rates follow a similar trend to that in a with the coupling between the qBIC and the graphene plasmon, as well as the effect of a secondary graphene plasmon from the side graphene ribbons, playing a significant role in the higher fermi levels close to 1 eV.

Supplementary Note 2. Structure parameter optimization
Since the optical response of our graphene -Si structure involves complicated dynamics of both qBIC and graphene plasmon resonances, the operating wavelength and the corresponding geometric parameters have to be numerically optimized and fine-tuned. For this purpose, we employed Reticolo RCWA (rigorous coupled-wave analysis) software (44,45) and the gradientfree algorithm BOBYQA (bound optimization by quadratic approximation) (46,47) to obtain the exact geometric parameters for maximal phase modulation.
First, we defined a figure of merit (FoM) in order to quantify the phase modulation performance.
The FoM was defined as the area enclosed by the complex reflectivity curve drawn on the complex plane as the graphene Fermi level being tuned at a single target frequency. This definition facilitated with the simultaneous optimization of both the reflection amplitude and the phase coverage.
Second, we utilized the gradient-free algorithm BOBYQA for the structural parameter optimization. Gradient-free algorithms are suitable for problems in which calculating the objective function gradient over the design variable space is infeasible or impossible. Our structure exhibiting complicated dynamics with respect to changes in geometry that is analytically intractable, making gradient-free algorithm a preferred optimization algorithm over gradient-based algorithms.
To complete the optimization program we had to find a suitable simulation tool for our structure.
We found that RCWA method was proper for this task, since its simulation times are much faster than that of FDTD (finite-difference time-domain) or FEM (finite element method) methods, while providing satisfactory accuracy. For the graphene sheet case, the Fourier order needed to reproduce similar reflectance spectra to FEM results was only 100, for which evaluating 51 Fermi energy points took less than 5 seconds on our server computer (Intel E5-2680v4). However, for the graphene ribbon case, the required Fourier order was about 600, which took roughly 200 times longer than that of the case with 100 Fourier orders. The reason for such a long computation time stems from the extremely narrow width of the graphene nanoribbons compared to the structure period. In order to bypass this problem, we first optimize geometric structural parameters for a maximum FOM at low Fourier orders, and then set the obtained geometric parameters as an initial point for the next optimization procedure with higher Fourier orders. This gradual increase of the Fourier orders rendered the simulation more accurate.
Finally, the simulation results generated by the RCWA was compared to that of the commercial FEM tool (COMSOL Multiphysics), which demonstrated excellent agreement.

Supplementary Note 4. Comparison of the theoretical TCMT model and simulation results
With the analytical two-resonance model for the complex reflection amplitude , we can compare the theoretical reflectance spectra for an avoided crossing with two resonances and that of the simulation obtained from full wave simulations (Fig. 3 of the main paper). To model the resonance frequency of the qBIC as a function of fermi level, we assumed that the qBIC resonance follows a linear trend due to the perturbative change of the graphene permittivity with respect to a much larger qBIC mode volume spanned around the Si pillars. The fitting parameters were chosen as

Supplementary Note 6. Phase modulation performance dependence on the TCMT parameters
We investigated whether it was better performance-wise for the spectrally narrow resonance to have the smallest FWHM as possible and for the second resonance to have the greatest ∆ as possible. In general, the decreasing FWHM of the spectrally narrow resonance (denoted as Resonance 1 hereafter) with respect to ∆ of the second resonance (denoted as Resonance 2 hereafter) leads to a better device performance (See Fig. 1b main text). However, there exist a few additional requirements that must be satisfied simultaneously: • Resonance 1 should be highly over-coupled ( ≫ ): The FWHM of Resonance 1 can be decomposed to FWHM1 = 2( 1 + 1 ) , where 1 and 1 are the radiative and dissipative coupling rate of Resonance 1. Even if the FWHM1 is very small, if Resonance 1 is under-coupled ( 1 < 1 ) the device will fail to cover the range 0-2π in phase modulation as shown in Fig.1a and Supplementary Fig. 4b. A decrease in 1 while 1 is kept constant, on the other hand, is beneficial not only because it decreases the FWHM 1 but also increases the ratio = / , thereby increasing the reflection amplitude through = | | 2 = 1 − 4 /( + 1) 2 and making it more uniform as shown in Supplementary Fig. 4c. It is evident that the complex reflection amplitude draws a larger phase modulation circle while having a more uniform radius in Supplementary Fig. 4c than the reference case in Supplementary Fig.   4a.
• should not be too small compared to : Even if it's ensured that Resonance 1 stays highly over-coupled ( 1 ≫ 1 ), decreasing 1 does not always bring performance enhancement. If 1 is much smaller than 2 of Resonance 2 in terms of their order of magnitudes ( 1 ≪ 2 ), the phase modulation will fail to cover the range 0-2π as shown in Supplementary Fig. 4d, a case where both 1 and 1 are decreased while keeping 2 the same as the reference. This is because the hybridized mode (with the characteristics of having both 1 and 2 ) will be too under-coupled near the avoided crossing region and the phase modulation will fail to circle around the origin.
• Resonance 2 should be highly tunable (∆ > FWHM1): Having a greater ∆ 2 will lead to a better phase modulation performance. This is because the hybridized modes will cross the operating frequency in greater speed as the control parameter (in our case the Fermi level) is varied. This is illustrated in Supplementary Fig. 4e.
From the discussions above, we can combine the different cases and plot the ideal case of having large ∆ 2 and small 1 , 1 , 2 with 1 ≫ 1 & 1 > 2 , as shown in Supplementary Fig.   4f.  Table S1 below. c, The case where 1 = 1 0 and

Supplementary Note 7. Relocation of the avoided crossing for limited Fermi level tuning range
The device presented in the main paper can only cover 282° for EF = 0-0.6eV as shown in Fig. 3d.
The reason for this sub-2π phase coverage is that the device utilizes only one resonance in this EF range since the avoided crossing between the graphene plasmon resonance and qBIC mode occurs at the upper bound of the EF range (EF = 0.6 eV). As mentioned before, achieving full 2π coverage requires having a double crossover between the operating frequency line and the two hybrid modes.
Therefore, to have a better phase coverage with a limited EF tuning range, the device parameters should be altered to have the avoided crossing at a lower Fermi level. This can be done by either narrowing the center graphene width to obtain the GP resonance at a lower Fermi energy or by increasing the meta-atom period to lower the frequency of the qBIC resonance. The other design parameters, which determines the radiative and dissipative decay rates of the resonances, also need to be fine-tuned to achieve a nearly uniform reflection amplitude.
To test this possibility, we conducted an optimization for 0-0.6eV Fermi level range as shown in Supplementary Fig. 5. Compared to the structure presented in the main paper (blue line in Supplementary Fig. 5a), the low-EF-optimized device (red line in Supplementary Fig. 5a) has a 15.6% longer period Λ but the gap between the Si pillars is nearly the same for the two structures.
As shown in Supplementary Fig. 5a, the reflection curve of the low-EF-optimized device can cover full 2π phase change while its magnitude is slightly lower than the curve in Fig. 3d. Note that, as shown in Supplementary Fig. 5b, c, the avoided crossing occurs at around EF = 0.5 eV, which is included in the 0-0.6 eV tuning range.