# Melting of unidirectional charge density waves across twin domain boundaries
in GdTe3

Sanghun Lee  These authors contributed equally to this work  Department of
Physics, Yonsei University, Seoul 03722, Republic of Korea  Eunseo Kim  These
authors contributed equally to this work  Department of Physics, Yonsei
University, Seoul 03722, Republic of Korea  Junho Bang  Department of Physics,
Yonsei University, Seoul 03722, Republic of Korea  Jongho Park  Center for
Correlated Electron Systems, Institute for Basic Science, Seoul 08826,
Republic of Korea  Department of Physics and Astronomy, Seoul National
University, Seoul 08826, Republic of Korea  Changyoung Kim  Center for
Correlated Electron Systems, Institute for Basic Science, Seoul 08826,
Republic of Korea  Department of Physics and Astronomy, Seoul National
University, Seoul 08826, Republic of Korea  Dirk Wulferding
<EMAIL_ADDRESS>Center for Correlated Electron Systems, Institute for Basic
Science, Seoul 08826, Republic of Korea  Department of Physics and Astronomy,
Seoul National University, Seoul 08826, Republic of Korea  Doohee Cho
<EMAIL_ADDRESS>Department of Physics, Yonsei University, Seoul 03722,
Republic of Korea

###### Abstract

Solids undergoing a transition from order to disorder experience the
proliferation of topological defects. The melting process generates transient
quantum states. However, their dynamical nature with femtosecond lifetime
hinders exploration with atomic precision. Here, we suggest an alternative
approach to the dynamical melting process by focusing on the interface created
by competing degenerate quantum states. We use a scanning tunneling microscope
(STM) to visualize the unidirectional charge density wave (CDW) and its
spatial progression (”static melting”) across a twin domain boundary (TDB) in
the layered material GdTe3. Combining STM with a spatial lock-in technique, we
reveal that the order parameter amplitude attenuates with the formation of
dislocations and thus two different unidirectional CDWs coexist near the TDB,
reducing the CDW anisotropy. Notably, we discover a correlation between this
anisotropy and the CDW gap. Our study provides valuable insight into the
behavior of topological defects and transient quantum states.

Topological defects are defined as singularities that appear in the order
parameter space[1]. They have been of importance in condensed matter physics
due to their key role in understanding the emergence of exotic quantum phases
in materials, such as those involving superconductivity[2, 3, 4, 5],
ferroicity[6, 7, 8], charge density wave (CDW)[9, 10, 11, 12, 13], and so on.
In general, they appear when a system undergoes phase transitions that change
its charge(or spin) order or symmetry[14, 15]. Their proliferation often gives
rise to phenomena different from those inherent in the pristine lattice[16].
These characteristics make topological defects highly relevant to non-
equilibrium states driven by external control parameters[17, 18, 19].

The family of RTe3 (R = rare earth elements) is notable for its potential to
serve as an ideal platform for exploring the mechanism of phase transitions
and proliferation of topological defects, as it exhibits various broken
symmetry ground states including CDW[20], antiferromagnetic spin ordering[21],
and superconductivity[22]. In particular, some of them are known for their
unidirectional CDW, while hydrostatic pressure can generate a bidirectional
CDW with a superconducting ground state[23, 24]. Previous studies on RTe3 have
reported that light pulses can also induce melting of the equilibrium CDW
state and the temporal evolution of an out-of-equilibrium CDW[25, 26, 27, 28].
These transient states host the bidirectional CDW demonstrated by the pump-
and-probe diffraction[26, 27, 28]. It has been proposed that the phase
transitions occurring in such non-equilibrium regions are closely related to
the proliferation of topological defects[29]. However, since it is challenging
to acquire microscopic details of the ultra-fast transient states, information
about the spatial distribution of the topological defects, its relation with
the order parameter amplitude, and the electronic structure of the excited
states has been considerably limited so far.

Here, we present the spatial variations of CDW across twin domain boundaries
(TDB), where two orthogonally aligned unidirectional CDW domains compete. We
utilized low-temperature scanning tunneling microscopy and spectroscopy
(STM/S) to visualize the charge modulations and electronic structure across
the domain boundary. Our analysis employs a spatial lock-in technique to
demonstrate the static CDW melting with a spatial gradient of its order
parameter across the TDB and the proliferation of topological defects. We also
find that two unidirectional CDWs coexist at the TDB, and that the normalized
difference in their amplitude anti-correlates with the density of states (DOS)
close to the Fermi level ($E_{\rm F}$). These findings offer a microscopic
understanding of the formation of topological defects and their accompanying
melting process, while also suggesting a connection to the superconducting
phase of this system.

Figure 1: Unidirectional CDW in GdTe3. (a) Crystal structure of GdTe3. The
semiconducting GdTe slab is encapsulated between conducting Te layers. The
dashed line between Te layers indicates the cleavage plane. (b) An STM
topographic image of the cleaved surface of GdTe3 acquired with $V_{\rm
set}=-100$ mV and $I_{\rm set}=100$ pA. $Scale$ $bar$, 5 nm. (c) Zoomed-in
view of the area marked by the black square in (b). The Gd lattice is resolved
with a unidirectional CDW propagating along the $c-$axis. (d) FFT of the STM
topographic image shown in (b). Orange and light-blue circles indicate the
Bragg peaks of Gd lattice along $a-$ and $c-$axes. Dark blue arrows are
associated with the wavevectors of the unidirectional CDW along the $c-$axis.
$Scale$ $bar$, 1/Å. (e) The $dI/dV$ tunneling spectrum of GdTe3 taken with
$V_{\rm set}=-400$ mV, $V_{\rm mod}=10$ mV, and $I_{\rm set}=1$ nA. Two dark
blue arrows mark the band edges of the partially opened CDW gap. One in the
empty state is determined by the spectral peak, and the other in the filled
state by an assumption of particle-hole symmetry. To make a comparison, the
gaps identified by STM [30, 31] and optical spectroscopy [32] are marked with
red and orange triangles, respectively. All color bars are linear scales.

GdTe3 is a layered material comprising two conducting square Te nets separated
by an insulating GdTe slab, as illustrated in Figure 1(a). The Te layers are
weakly coupled by van der Waals interactions, suggesting that a Te layer
terminates each cleavage plane. This metallic surface layer is associated with
the unidirectional CDW observed in previous STM studies. An STM topographic
image of the cleaved surface, shown in Figure 1(b), reveals a square lattice
whose unit cell is consistent with the Gd lattice. This can be rationalized by
its protrusion closer to the surface and its higher charge density compared to
Te atoms in the insulating GdTe layer, resulting in a higher STM current [33].
The image also exhibits regular stripe patterns with wavefronts propagating
along the $c-$axis. This unidirectional charge modulation has been known as an
incommensurate CDW phase with a wavelength of about $7/2c$, where $c$ is the
lattice constant along the direction perpendicular to the wavefronts.
Superimposing the CDW superstructure with the Gd lattice results in
unidirectional patterns with a periodicity of about $7c$, as indicated by the
black arrow in Figure 1(c).

The Bragg peaks for the Gd lattice and the unidirectional CDW are clearly
visible in the fast Fourier transform (FFT) of the STM topographic image. The
Bragg peaks for the pristine lattice are highlighted by orange and light-blue
circles in Figure 1(d). Additionally, the satellite peaks marked by blue
arrows correspond to the unidirectional CDW. As shown in Figure 1(a), the bulk
GdTe3 has a glide symmetry along the $c-$axis rather than the $a-$axis, with
the mirror plane between the Te layers. Due to the anisotropy in the crystal
symmetry, the lattice constant $c$ is slightly larger than $a$, and thus the
unidirectional CDW forms preferably along the $c-$axis[22]. Therefore, the
orientation of the CDW may serve as a reference for determining the direction
of the lattice vectors. By fitting the Bragg peaks in the FFT, we determine
that the lattice constants is $a(c)\sim 4.3$ Å, and the incommensurate CDW
wavevector $\textbf{q}_{\rm CDW}\sim 2/7$qc. These results are consistent with
previous studies[34].

The differential conductance ($dI/dV$) spectrum (Figure 1(e)) allows an
estimate of the CDW gap of about $2\Delta\approx 420$ meV[30, 31, 32], which
is determined by assigning the band edges marked with dark blue arrows. In
addition, the derivative of the spectrum highlights a subtle change around the
assigned edge of the valence band (See Figure 5). In 2D metallic systems, a
CDW gap opens when the states are connected with the nesting vector
corresponding to the CDW wavevector (qCDW). Therefore, a CDW state in 2D
metallic systems is likely to have a partially gapped Fermi surface, resulting
in the finite DOS at $E_{\rm F}$. The Fermi surface of GdTe3 exhibits a
diamond-like structure, which is the result of the combination of Fermi
surfaces originating from two perpendicular 1D metallic chains. Those are
attributed to the inter-site hopping of $p_{\rm x}$ and $p_{\rm y}$ orbitals
and the weak inter-orbital interaction between them in the Te square net [35,
36]. This Fermi surface topology satisfies the unidirectional nesting
condition with band-folding driven by the subsurface Gd-lattice. In accordance
with our data, recent studies using angle-resolved photoemission spectroscopy
(ARPES) have revealed an anisotropic gap opening and the remnant Fermi surface
at low temperatures[30, 37]. In addition, charge modulations between the
conduction and the valence band edge exhibit an out-of-phase relation, which
has been commonly observed in Peierls-type CDW systems[38, 39, 40, 41] (See
Figure 6).

Figure 2: A TDB with perpendicularly oriented unidirectional CDWs. (a) An STM
topographic image of two different unidirectional CDW domains and a TDB
acquired with $V_{\rm set}=400$ mV and $I_{\rm set}=300$ pA. $Scale$ $bar$, 10
nm. (b) FFT of the STM topographic image shown in (a). Orange and light-blue
circles highlight the Bragg peaks of the Gd lattice. Red and blue circles are
associated with the Bragg peaks of the unidirectional CDW propagating along
the $a-$ and $c-$axes, respectively. $Scale$ $bar$, 1/Å. (c) Line profiles of
the FFT intensity in (b) along $a-$ (red line) and $c-$axes (blue line). The
wavevectors are normalized with the reciprocal lattice vectors for the same
direction. (d-f) Zoomed-in view of the area marked by the white squares in
(a). $Scale$ $bar$, 2 nm. All maps in (d-f) were spatially differentiated. All
color bars are linear scale. Figure 3: Visualizing amplitude and phase of the
unidirectional CDW. (a) and (b) Phase maps for the unidirectional CDWs along
the $c-$ and $a-$axes, respectively, obtained by the spatial lock-in
technique. (c) and (d) The CDWs reconstructed by combining the amplitude and
phase acquired with reference signals corresponding to the unidirectional CDW
along $c-$ and $a-$ axes. (e) and (f) The zoomed-in images of topological
defect pairs in the area marked with white squares in (a(c)) and (b(d)),
respectively. The singular points where phase winding occurs in $\pm 2\pi$,
highlighted with black (clockwise) and white (counterclockwise) dots, are
characterized by the CDW dislocations. All maps in (a-d) were acquired within
the same field of view as the STM topographic image in Figure 2(a). $Scale$
$bar$, 10 nm. All color bars are linear scale.

X-ray diffraction[42], transport[20], and STM[30, 31] studies have reported
that GdTe3 has a unidirectional CDW below $T_{\rm c}\sim 340$ K. This
anisotropy has been attributed to the broken $C_{4}$ rotation symmetry
(—$\mathbf{q}_{\rm a}$— ¿ —$\mathbf{q}_{\rm c}$—) inherent in the crystal
structure. We confirm this symmetry breaking by additional polarization-
resolved Raman spectroscopy measurements[43], which can directly probe the
orientation of unidirectional CDW domains on a larger micrometer-scale (see
Figure 7). However, external control parameters, such as strain[44] or laser
pulses[27], may enable the coexistence of two perpendicularly oriented
unidirectional CDW states domain by domain. These degenerate CDWs share the
same crystallographic structure but propagate along orthogonal directions. The
CDW domains are separated by the interface, which we have previously referred
to as twin domain boundaries (TDB). In the following, we outline our
observation and characterization of such TDBs.

Our STM topographic image (Figure 2(a)) illustrates the presence of the TDB,
where two CDW domains with perpendicular orientations intersect. The arrows
with ticks in Figure 2(a) indicate the direction of propagation and wavefronts
of the CDW for each domain, respectively. Notably, the TDB in the center line
of the STM topographic image exhibits a bidirectional CDW resulting from the
superposition of two unidirectional CDWs. The differentiated STM images
(Figure 2(d)-(f)) provide a clearer representation of two distinct CDWs and
their coexistence at the TDB. These unidirectional CDWs are not sharply
separated but spread across the TDB with amplitude attenuation.

In Figure 2(b), the FFT of our STM topographic image shows two sets of
distinct CDW peaks near the $\Gamma$ point, each representing unidirectional
CDWs along the $a-$axis (red circles) and $c-$axis (blue circles),
respectively. Further measurements demonstrate that the CDW along the $c-$axis
is dominant throughout the entire surface (see Figure 8), while the CDW along
the $a-$axis is locally formed. It has been well established that the
direction of the CDW is determined by the anisotropy of the pristine
lattice[22]. However, as shown in Figure 2(b), the Bragg peaks of the Gd
lattice, highlighted by orange and light-blue circles, are sharply defined
without peak splitting since the structural anisotropy is beyond our
experimental resolution.

To compare two perpendicular CDWs, we acquired line profiles of the FFT
intensity along $a-$ and $c-$axes in Figure 2(c). The wavevectors are
normalized with the reciprocal lattice vectors of GdTe3. They commonly show
the strongest peak and four additional peaks at 1, 2/7, 3/7, 4/7, and 5/7,
corresponding to $\mathbf{q}_{\rm a(c)}$, $\mathbf{q}_{\rm CDW}$,
1-2$\mathbf{q}_{\rm CDW}$, 2$\mathbf{q}_{\rm CDW}$, and 1-qCDW, respectively.
The normalized two different CDW wavevectors are identical within the
resolution of our data.

When examining the spatial variation across the TDB, we can observe a gradual
diminishing of the CDW in each domain as it crosses into the other domain. The
process, which is termed ”static melting” in this paper, can be accompanied by
a reduction of the CDW amplitude and a twist of its phase. To acquire the
spatial variations of the amplitude($A_{\rm CDW}(\mathbf{r})$) and the
phase(${\phi_{\rm CDW}(\mathbf{r})}$) of the unidirectional CDW, we employed a
spatial lock-in technique (see details in Figures 9 and 10). Figure 3(a) and
(b) display changes in the phase of the individual unidirectional CDWs, which
propagate along the $c-$ and $a-$axes, respectively. We note that the CDWs are
almost coherent in their domains, while several phase windings of $\pm 2\pi$,
so-called topological defects, are observable across the TDB.

Topological defects in the form of $2\pi$ phase winding are known to induce
dislocations in the CDW wavefront [1]. We can construct the individual
unidirectional CDWs with $A_{\rm CDW}(\mathbf{r})$ and ${\phi_{\rm
CDW}(\mathbf{r})}$ acquired by the lock-in technique. Figure 3(c) and (d)
clearly show the attenuation of the CDW amplitude and the proliferation of
dislocations across the TDB. These results suggest that the phase stiffness of
the unidirectional CDW can be reduced when its order parameter amplitude
decreases[45]. In addition, by comparing the spatial relationship between
their phase and amplitude, it is evident that phase singularities occur where
the CDW amplitude is strongly reduced (Figure 3(e) and (f)). The local
excitation is unstable against the CDW reconstruction, either restoring the
initial state or creating the $90^{\circ}$ rotated unidirectional CDW. In that
case, the system exhibits the bidirectional CDW[27].

Figure 4: Electronic structure of the TDB. (a) The spatial derivative of the
STM topographic image of a TDB region acquired simultaneously with the two-
dimensional spectroscopic mapping. $Scale$ $bar$, 5 nm. (b) The line cut of
$dI/dV$ map along a yellow arrow in (a). (c) Averaged spectra acquired in the
TDB (red line) and in the CDW domain (black line). Spectral weight transfer
leads to the in-gap like feature close to $E_{F}$ and also suppression of the
DOS around +300 meV. In addition, no gap opening is apparent at higher binding
energy. (d) In-gap spectral weight map of the area marked by the red dashed
square in (a). In-gap spectral weight is obtained by integrating $dI/dV$ from
$-100$ to $+100$ meV. Tunneling conditions for $dI/dV$ mapping are $V_{\rm
set}=-400$ mV and $I_{\rm set}=1$ nA. $Scale$ $bar$, 5 nm. (e) Anti-
correlation between the anisotropy parameter and the in-gap spectral weight
along the line marked by a yellow arrow in (a). All color bars are linear
scale.

The gradual attenuation of each unidirectional CDW amplitude across the TDB
allows us to investigate the difference between unidirectional and
bidirectional CDWs shown in Figure 4(a). We performed spectroscopic mapping in
the vicinity of the TDB. Figure 4(b) shows differential conductance ($dI/dV$)
spectra acquired along the yellow arrow in Figure 4(a). The $dI/dV$ spectra
averaged over the TDB and the domains are compared in Figure 4(c). The direct
comparison demonstrates two notable distinctions: a reduced CDW gap and an
enhanced DOS near $E_{\rm F}$ at the TDB. Similar observations have been
reported in other RTe3 samples, in cases where checkerboard CDW is formed
through phase transitions[32], or when a bidirectional CDW is induced by laser
pulses[26].

We further analyze the spatial evolution of DOS near $E_{\rm F}$ by
calculating the in-gap spectral weight ($\Sigma dI/dV$). For a direct
comparison with previous time-resolved ARPES studies[26], we integrate the
differential conductance from $-100$ meV to $+100$ meV, which is within the
CDW gap. As shown in Figure 4(d), the in-gap spectral weight exhibits a local
enhancement at the TDB, while the bidirectional CDW states exist across a
spatially extended region and topological defects are scattered in the
proximity of the TDB. To correlate the localized spectral feature with the CDW
order parameter, we determined the amplitude of the CDW, $A_{\rm a}$ and
$A_{\rm c}$, as a function of position (see details in Figure 11).
Subsequently, we define an anisotropy parameter $\rm{AP}$ as follows:
$\rm{AP}=\left|\frac{A_{\rm c}-A_{\rm a}}{A_{\rm c}+A_{\rm a}}\right|$. The
anisotropy parameter is close to 1 for the unidirectional CDW and 0 for the
bidirectional (checkerboard-like) CDW. As shown in Figure 4(e), the anisotropy
parameter strongly anti-correlates with the in-gap spectral weight in Figure
4(d).

Note that such bidirectional CDWs in $R$Te3 have been commonly observed in
excited states driven by external stimuli with topological defects[26, 27].
They have been reported to have higher LDOS than their ground states around
$E_{F}$. However, these quantum states are neither static nor spatially
homogeneous, which hinders an accurate determination of their atomic and
electronic properties. In this regard, we highlight the TDB, which can be
locally pinned by strain gradients and exist at the interface between two
perpendicularly aligned unidirectional CDWs. Examining the TDB allows us to
observe the decay behavior statically of the order parameter amplitude and
phase, which is likely to exhibit a charge order similar to that of the
excited state. Remarkably, we find that the electronic structure is primarily
associated with the relative difference of the CDW amplitude between
degenerate ground states, rather than with their phase. This is consistent
with the decay behavior of the laser-induced excited states in previous time-
resolved ARPES experiments [26], highlighting that the restoration of the CDW
gap occurs much faster than its phase coherence. Therefore, this approach is
an ideal tool for investigating the electronic properties of local
excitations, such as bidirectional CDW or topological defects, to better
understand short-lived quantum states out of equilibrium.

In conclusion, our STM results reveal the spatial progression of the
unidirectional CDW close to a TDB in GdTe3. Our spatial lock-in technique
enables us to find the local excitations and their distribution which so far
have been extensively discussed for understanding the dynamical melting
process of the unidirectional CDW without microscopic experimental evidence.
We note that the amplitude of the CDW order parameter decreases with winding
its phases, thus constructing the unidirectional CDW wavefronts’ dislocations
and the bidirectional CDW. This intermediate region where two unidirectional
CDWs coexist exhibits the enhanced DOS close to $E_{\rm F}$, which was
observed in pump-probe photoelectron spectroscopy studies [26]. The higher
electron density near $E_{\rm F}$, accompanied by a bidirectional CDW, might
be connected to the emergence of superconductivity in the bidirectional (or
checkerboard) CDW in RTe3 systems [20, 46, 47]. Lastly, our approach will be
instructive to study competing phenomena between charge(or spin)-ordered
phases and superconductivity in low-dimensional correlated systems. To gain a
deeper understanding of this complex relationship, it is required to conduct
additional STM experiments on the systems that exhibit the charge order with
superconductivity below the critical temperatures. These atomically resolved
experiments should focus on investigating the connection between the CDW and
the superconducting order parameter on the nanoscale [48].

###### Acknowledgements.

The authors thank M. P. Allan, T. Benschop, J. -F. Ge, G. Kim, S, Kim, J.
Rhim, S. Park, and Y. Son for valuable discussions. S. L., E. K., J. B., and
D. C. were supported by the National Research Foundation of Korea (NRF) grant
funded by the Korea government (MSIT) (No. 2020R1C1C1007895, 2017R1A5A1014862,
and RS-2023-00251265) and the Yonsei University Research Fund of 2019-22-0209.
J. P., C. K., and D. W. acknowledge support from the Institute for Basic
Science (IBS) (Grants No. IBS-R009-G2, IBS-R009-Y3).

## References

  * [1] Mermin, N. D. The topological theory of defects in ordered media. Rev. Mod. Phys. 1979, 51, 591–648.
  * [2] Agterberg, D.; Tsunetsugu, H. Dislocations and vortices in pair-density-wave superconductors. Nat. Phys. 2008, 4, 639–642.
  * [3] Du, Z.; Li, H.; Joo, S. H.; Donoway, E. P.; Lee, J.; Davis, J. S.; Gu, G.; Johnson, P. D.; Fujita, K. Imaging the energy gap modulations of the cuprate pair-density-wave state. Nature 2020, 580, 65–70.
  * [4] Chen, W.; Ren, W.; Kennedy, N.; Hamidian, M.; Uchida, S.-i.; Eisaki, H.; Johnson, P. D.; O’Mahony, S. M.; Davis, J. S. Identification of a nematic pair density wave state in Bi2Sr2CaCu2O8+x. Proc. Natl. Acad. Sci. U.S.A. 2022, 119, e2206481119.
  * [5] Liu, Y.; Wei, T.; He, G.; Zhang, Y.; Wang, Z.; Wang, J. Pair density wave state in a monolayer high-Tc iron-based superconductor. Nature 2023, 618, 934–939.
  * [6] Chae, S.; Lee, N.; Horibe, Y.; Tanimura, M.; Mori, S.; Gao, B.; Carr, S.; Cheong, S. W. Direct observation of the proliferation of ferroelectric loop domains and vortex-antivortex pairs. Phys. Rev. Lett. 2012, 108, 167603.
  * [7] Artyukhin, S.; Delaney, K. T.; Spaldin, N. A.; Mostovoy, M. Landau theory of topological defects in multiferroic hexagonal manganites. Nat. Mater. 2014, 13, 42–49.
  * [8] Liu, J.; Balents, L. Anomalous Hall effect and topological defects in antiferromagnetic Weyl semimetals: Mn3Sn/Ge. Phys. Rev. Lett. 2017, 119, 087202.
  * [9] Le Bolloc’h, D.; Ravy, S.; Dumas, J.; Marcus, J.; Livet, F.; Detlefs, C.; Yakhou, F.; Paolasini, L. Charge density wave dislocation as revealed by coherent x-ray diffraction. Phys. Rev. Lett. 2005, 95, 116401.
  * [10] Pásztor, Á.; Scarfato, A.; Spera, M.; Barreteau, C.; Giannini, E.; Renner, C. Holographic imaging of the complex charge density wave order parameter. Phys. Rev. Res. 2019, 1, 033114.
  * [11] Lim, S.; Kim, J.; Won, C.; Cheong, S.-W. Atomic-scale observation of topological vortices in the incommensurate charge density wave of 2H-TaSe2. Nano Lett. 2020, 20, 4801–4808.
  * [12] Aishwarya, A.; May-Mann, J.; Almoalem, A.; Ran, S.; Saha, S. R.; Paglione, J.; Butch, N. P.; Fradkin, E.; Madhavan, V. Visualizing the melting of the charge density wave in UTe2 by generation of pairs of topological defects with opposite winding. 2023, arXiv:2306.09423. [cond-mat.supr-con]. https://doi.org/10.48550/arXiv.2306.09423 (accessed Nov. 06, 2023).
  * [13] LaFleur, A.; Li, H.; Frank, C. E.; Xu, M.; Cheng, S.; Wang, Z.; Butch, N. P.; Zeljkovic, I. Inhomogeneous high temperature melting and decoupling of charge density waves in spin-triplet superconductor UTe2. 2023, arXiv:2308.03721. [cond-mat.supr-con]. https://doi.org/10.48550/arXiv.2308.03721 (accessed Nov. 06, 2023).
  * [14] Saito, H.; Kawaguchi, Y.; Ueda, M. Topological defect formation in a quenched ferromagnetic Bose-Einstein condensates. Phys. Rev. A 2007, 75, 013621.
  * [15] Du, K.; Fang, X.; Won, C.; De, C.; Huang, F.-T.; Xu, W.; You, H.; Gómez-Ruiz, F. J.; del Campo, A.; Cheong, S.-W. Kibble–Zurek mechanism of Ising domains. Nat. Phys. 2023, 19, 1495-1501.
  * [16] Mishra, S.; Beyer, D.; Berger, R.; Liu, J.; Groning, O.; Urgel, J. I.; Mullen, K.; Ruffieux, P.; Feng, X.; Fasel, R. Topological defect-induced magnetism in a nanographene. J. Am. Chem. Soc. 2020, 142, 1147–1152.
  * [17] Vogelgesang, S.; Storeck, G.; Horstmann, J.; Diekmann, T.; Sivis, M.; Schramm, S.; Rossnagel, K.; Schäfer, S.; Ropers, C. Phase ordering of charge density waves traced by ultrafast low-energy electron diffraction. Nat. Phys. 2018, 14, 184–190.
  * [18] Cheng, Y.; Zong, A.; Wu, L.; Meng, Q.; Xia, W.; Qi, F.; Zhu, P.; Zou, X.; Jiang, T.; Guo, Y., et al. Ultrafast formation of topological defects in a 2D charge density wave. 2022, arXiv:2211.05748. [cond-mat.mtrl-sci]. https://doi.org/10.48550/arXiv.2211.05748 (accessed Nov. 06, 2023).
  * [19] Wandel, S.; Boschini, F.; da Silva Neto, E.; Shen, L.; Na, M.; Zohar, S.; Wang, Y.; Welch, S.; Seaberg, M.; Koralek, J., et al. Enhanced charge density wave coherence in a light-quenched, high-temperature superconductor. Science 2022, 376, 860–864.
  * [20] Ru, N.; Condron, C.; Margulis, G.; Shin, K.; Laverock, J.; Dugdale, S.; Toney, M.; Fisher, I. Effect of chemical pressure on the charge density wave transition in rare-earth tritellurides RTe3. Phys. Rev. B 2008, 77, 035114.
  * [21] Ru, N.; Chu, J.-H.; Fisher, I. Magnetic properties of the charge density wave compounds RTe3 (R = Y, La, Ce, Pr, Nd, Sm, Gd, Tb, Dy, Ho, Er, and Tm). Phys. Rev. B 2008, 78, 012410.
  * [22] Yumigeta, K.; Qin, Y.; Li, H.; Blei, M.; Attarde, Y.; Kopas, C.; Tongay, S. Advances in rare-earth tritelluride quantum materials: structure, properties, and synthesis. Adv. Sci. 2021, 8, 2004762.
  * [23] Hamlin, J.; Zocco, D.; Sayles, T.; Maple, M.; Chu, J.-H.; Fisher, I. Pressure-induced superconducting phase in the charge-density-wave compound terbium tritelluride. Phys. Rev. Lett. 2009, 102, 177002.
  * [24] Zocco, D.; Hamlin, J.; Grube, K.; Chu, J.-H.; Kuo, H.-H.; Fisher, I.; Maple, M. Pressure dependence of the charge-density-wave and superconducting states in GdTe3, TbTe3, and DyTe3. Phys. Rev. B 2015, 91, 205114.
  * [25] Schmitt, F.; Kirchmann, P. S.; Bovensiepen, U.; Moore, R. G.; Rettig, L.; Krenz, M.; Chu, J.-H.; Ru, N.; Perfetti, L.; Lu, D., et al. Transient electronic structure and melting of a charge density wave in TbTe3. Science 2008, 321, 1649–1652.
  * [26] Zong, A.; Kogar, A.; Bie, Y.-Q.; Rohwer, T.; Lee, C.; Baldini, E.; Ergeçen, E.; Yilmaz, M. B.; Freelon, B.; Sie, E. J., et al. Evidence for topological defects in a photoinduced phase transition. Nat. Phys. 2019, 15, 27–31.
  * [27] Kogar, A.; Zong, A.; Dolgirev, P. E.; Shen, X.; Straquadine, J.; Bie, Y.-Q.; Wang, X.; Rohwer, T.; Tung, I.-C.; Yang, Y., et al. Light-induced charge density wave in LaTe3. Nat. Phys. 2020, 16, 159–163.
  * [28] Zhou, F.; Williams, J.; Sun, S.; Malliakas, C. D.; Kanatzidis, M. G.; Kemper, A. F.; Ruan, C.-Y. Nonequilibrium dynamics of spontaneous symmetry breaking into a hidden state of charge-density wave. Nat. Commun. 2021, 12, 566.
  * [29] Kosterlitz, J. M.; Thouless, D. J. Ordering, metastability and phase transitions in two-dimensional systems. J. Phys. C: Solid State Phys. 1973, 6, 1181.
  * [30] Lei, S.; Lin, J.; Jia, Y.; Gray, M.; Topp, A.; Farahi, G.; Klemenz, S.; Gao, T.; Rodolakis, F.; McChesney, J. L., et al. High mobility in a van der Waals layered antiferromagnetic metal. Sci. Adv. 2020, 6, eaay6407.
  * [31] Xu, Z.; Ji, S.-H.; Tang, L.; Wu, J.; Li, N.; Cai, X.; Chen, X. Molecular beam epitaxy growth and electronic structures of monolayer GdTe3. Chin. Phys. Lett. 2021, 38, 077102.
  * [32] Hu, B.; Cheng, B.; Yuan, R.; Dong, T.; Wang, N. Coexistence and competition of multiple charge-density-wave orders in rare-earth tritellurides. Phys. Rev. B 2014, 90, 085105.
  * [33] Tomic, A.; Rak, Z.; Veazey, J.; Malliakas, C.; Mahanti, S.; Kanatzidis, M.; Tessmer, S. Scanning tunneling microscopy study of the CeTe3 charge density wave. Phys. Rev. B 2009, 79, 085422.
  * [34] Chen, Y.; Wang, P.; Wu, M.; Ma, J.; Wen, S.; Wu, X.; Li, G.; Zhao, Y.; Wang, K.; Zhang, L., et al. Raman spectra and dimensional effect on the charge density wave transition in GdTe3. Appl. Phys. Lett. 2019, 115, 151905.
  * [35] Laverock, J.; Dugdale, S.; Major, Z.; Alam, M.; Ru, N.; Fisher, I.; Santi, G.; Bruno, E. Fermi surface nesting and charge-density wave formation in rare-earth tritellurides. Phys. Rev. B 2005, 71, 085114.
  * [36] Yao, H.; Robertson, J. A.; Kim, E.-A.; Kivelson, S. A. Theory of stripes in quasi-two-dimensional rare-earth tellurides. Phys. Rev. B 2006, 74, 245126.
  * [37] Liu, J.; Huan, S.; Liu, Z.; Liu, W.; Liu, Z.; Lu, X.; Huang, Z.; Jiang, Z.; Wang, X.; Yu, N., et al. Electronic structure of the high-mobility two-dimensional antiferromagnetic metal GdTe3. Phys. Rev. Mater. 2020, 4, 114005.
  * [38] Dai, J.; Calleja, E.; Alldredge, J.; Zhu, X.; Li, L.; Lu, W.; Sun, Y.; Wolf, T.; Berger, H.; McElroy, K. Microscopic evidence for strong periodic lattice distortion in two-dimensional charge-density wave systems. Phys. Rev. B 2014, 89, 165140.
  * [39] Spera, M.; Scarfato, A.; Pasztor, A.; Giannini, E.; Bowler, D. R.; Renner, C. Insight into the charge density wave gap from contrast inversion in topographic STM images. Phys. Rev. Lett. 2020, 125, 267603.
  * [40] Pásztor, Á.; Scarfato, A.; Spera, M.; Flicker, F.; Barreteau, C.; Giannini, E.; Wezel, J. v.; Renner, C. Multiband charge density wave exposed in a transition metal dichalcogenide. Nat. Commun. 2021, 12, 6037.
  * [41] Jang, W.-J.; Sim, J.; Heo, J. E.; Noh, H.; Kim, R.; Chang, S. H.; Jeon, S.; Kim, H. W. Direct observation of multiband charge density waves in NbTe2. Phys. Rev. B 2022, 106, 125110.
  * [42] Banerjee, A.; Feng, Y.; Silevitch, D.; Wang, J.; Lang, J.; Kuo, H.-H.; Fisher, I.; Rosenbaum, T. Charge transfer and multiple density waves in the rare earth tellurides. Phys. Rev. B 2013, 87, 155131.
  * [43] Wang, Y.; Petrides, I.; McNamara, G.; Hosen, M. M.; Lei, S.; Wu, Y.-C.; Hart, J. L.; Lv, H.; Yan, J.; Xiao, D., et al. Axial Higgs mode detected by quantum pathway interference in RTe3. Nature 2022, 606, 896–901.
  * [44] Straquadine, J.; Ikeda, M.; Fisher, I. Evidence for realignment of the charge density wave state in ErTe3 and TmTe3 under uniaxial stress via elastocaloric and elastoresistivity measurements. Phys. Rev. X 2022, 12, 021046.
  * [45] Emery, V.; Kivelson, S. Importance of phase fluctuations in superconductors with small superfluid density. Nature 1995, 374, 434–437.
  * [46] Fang, A.; Straquadine, J. A.; Fisher, I. R.; Kivelson, S. A.; Kapitulnik, A. Disorder-induced suppression of charge density wave order: STM study of Pd-intercalated ErTe3. Phys. Rev. B 2019, 100, 235446.
  * [47] Fang, A.; Singh, A. G.; Straquadine, J. A.; Fisher, I. R.; Kivelson, S. A.; Kapitulnik, A. Robust superconductivity intertwined with charge density wave and disorder in Pd-intercalated ErTe3. Phys. Rev. Res. 2020, 2, 043221.
  * [48] Liu, X.; Chong, Y. X.; Sharma, R.; Davis, J. S. Discovery of a Cooper-pair density wave state in a transition-metal dichalcogenide. Science 2021, 372, 1447–1452.

##  I Methods

Sample growth. We used stoichiometric single-crystalline GdTe3 flakes to grow
crystals through the self-flux method. We mixed high-purity Gd metal (99.9%)
and Te chips (99.999%) in a 1:30 molar ratio for the Te-rich self-flux
condition. The mixed precursor of GdTe3 was loaded into quartz tubes and
sealed at $\sim$10-5 Torr with a turbo pump to avoid oxygen contamination. We
heated the mixture to a maximum temperature of 900∘C and maintained the
temperature for 24 hours to achieve a homogeneous melt. Afterward, we cooled
the sample down to 500∘C at a rate of $-$2∘C/hour. We decanted the melted
mixture in an ampule and waited for it to solidify at room temperature. The
crystal structure was determined using a high-resolution X-ray diffractometer
with an airtight holder to prevent oxygen-caused degradation.

STM/S measurements. STM/S measurements were carried out using a commercial
low-temperature STM (UNISOKU) at liquid helium temperature ($T=4.2$ K) in
ultra-high vacuum (UHV) environment down to $1\times 10^{-10}$ Torr. The
single crystal samples of GdTe3 were exfoliated in situ at room temperature
and loaded to the precooled STM head. Mechanically sharpened Pt-Ir (90/10)
wires were used for STM tips. STM topographic images were obtained in constant
current mode by applying bias voltages ($V_{\rm set}$) to the sample. The
differential conductance ($dI/dV$) images were acquired using a standard lock-
in technique, with an a.c. voltage (a peak-to-peak voltage amplitude $V_{\rm
pp}=10$ mV and a frequency $f_{\rm mod}=613$ Hz added to the d.c. sample bias.

Raman scattering measurements. Raman scattering measurements shown in Figure 7
were carried out at room temperature (295 K) with a single mode laser (Oxxius,
561 nm). To avoid degradation of the surface, the sample was freshly cleaved,
installed, and sealed into the sample chamber while inside a glove box. After
removing the sealed sample chamber from the glove box, we evacuated it with a
turbo pump, thereby avoiding any air contamination. The polarization of the
incident laser light was rotated with a superachromatic $\lambda/2$ waveplate
within the crystallographic $ac$-plane. For the polarization-resolved Raman
maps, the polarization of the incident light and the polarization of the
scattered light entering the spectrometer are at a fixed parallel relation to
each other (”parallel polarization”). The inelastically Raman-scattered light
was dispersed through a Princeton Instruments TriVista spectrometer and
detected via a PyLoN eXcelon charge-coupled device.

##  II Supporting Information

Figure 5: Assignment of the CDW gap in the tunneling spectrum. Due to partial
gap opening, the CDW gap does not show the zero conductivity in the
differential tunneling spectrum. In the main text, we show the spatially
averaged tunneling spectrum in Figure 1(e). We assign the small peak (marked
by a dark blue arrow) in the empty state as the edge of the conduction band.
Due to the lack of a distinct shoulder-like feature in the filled state, we
determine the band edge at the energy where the rate of increase in the DOS
changes, shown in its derivative (red curve) and highlighted by the dashed
black lines. Figure 6: Phase relation of the energy dependent charge
modulations in the unidirectional CDW. a, A $dI/dV$ map of the unidirectional
CDW at the energy of -200 meV. It shows the charge modulations and the
quasiparticle interference from the defects. b, The FFT spectrum of the
$dI/dV$ map shown in a. This spectrum highlights the wavevectors of the
lattice and the unidirectional CDW, with colored circles indicating their
positions. c, The energy-dependent $dI/dV$ line profiles. They are acquired
along the perpendicular direction (red arrow in a) of the CDW wavefronts. All
profiles have been normalized for clarity. d, Phase relation of the charge
modulations in valence and conduction bands. Notably, the profiles at the
edges of the valence and conduction bands exhibit an out-of-phase relation, a
common feature observed in Peierls-type CDW systems. To gain a deeper
understanding of energy-dependent charge modulations, theoretical calculations
with the multiband structure, partially opened energy gap, and incommensurate
charge modulations superimposed on atomic corrugations are required. Figure 7:
Unidirectional CDW and its rotation. a, The prominent CDW amplitude mode
(marked by a green arrow) observed with Raman spectroscopy and its distinct
2-fold symmetry in the polarization-resolved Raman maps. A four-fold symmetric
phonon mode (marked by an orange arrow) used to assign crystallographic $a$
and $c$ axes, indicated by dashed black lines. Two different measurement
positions, a few tens $\mu$m apart from each other, show similar behavior,
although the amplitude mode is offset by a 90∘ phase between these two
positions ”1” and ”2”. b, Angle-dependent integrated Raman scattering
intensity of the CDW amplitude mode at position ”1” (blue squares) and
position ”2” (red circles). Note that the diameter of the laser spot (i.e.,
the surface area from which the signal was collected) was about 2 $\mu$m. The
observed two-fold symmetry of the CDW amplitude mode can be direct evidence
that larger, $\mu$m-sized domains exist in GdTe3. Figure 8: Coexistence of two
different unidirectional CDWs. a, Large-scale STM topographic image including
twin domains that have perpendicularly oriented unidirectional CDWs. One
domain with the CDW propagating along the $c-$axis is larger than the other.
b, A zoomed-in STM topographic image of the TDB. Figure 9: CDW phase analysis
with spatial lock-in technique. a, An STM image of the TDB shown in Figure
2(a): $T_{\rm raw}=\sum_{\rm j}\frac{A_{\rm j}(\mathbf{r})}{2}({\rm
e}^{\mathrm{i}[\mathbf{q}_{\rm j}\cdot\mathbf{r}+\phi_{\rm
j}(\mathbf{r})]}+\mathrm{e}^{-\mathrm{i}[\mathbf{q}_{\rm
j}\cdot\mathbf{r}+\phi_{\rm j}(\mathbf{r})]})$, where $\phi_{\rm
j}(\mathbf{r})$ is the position-dependent phase. b, Block diagram of the
spatial lock-in algorithm composed of a signal mixer and a low-pass filter. c
and d, The reference signals corresponding to the unidirectional CDWs: $T_{\rm
ref}={\rm e}^{-\mathrm{i}\mathbf{q}_{\rm CDW}\cdot\mathbf{r}}$. As shown in b,
$T_{\rm raw}$ is multiplied by $T_{\rm ref}$ with the mixer. The subsequent
low-pass filtering can extract the DC component from the mixing signal. In the
case of equal wavevector $\mathbf{q}_{\rm j}=\mathbf{q}_{\rm CDW}$, the lock-
in output given by $\operatorname{Re}(T_{\rm out})+{\rm
i}\operatorname{Im}(T_{\rm out})$. e and f, The amplitude ($A(\mathbf{r})$) of
the lock-in output. g and h, The phase ($\phi(\mathbf{r})$) of the lock-in
output can be obtained by calculating the arctangent of
$\operatorname{Im}(T_{\rm out})/\operatorname{Re}(T_{\rm out})$. i and j, The
real part ($A(\mathbf{r}){\rm cos}(\mathbf{q}_{\rm
CDW}\cdot\mathbf{r}+\phi(\mathbf{r}))$) of the lock-in output. Figure 10:
Proof of validity of the spatial lock-in technique. a, Phase maps for the
unidirectional CDWs along the $c-$(left panel) and $a-$(right panel)axes,
obtained by the spatial lock-in technique. The set of phase maps in the top
left corner is shown in Figures. 3(a) and 3(b). The other sets are acquired by
applying the lock-in technique to the STM images for the same field of view
with different scan directions (denoted in a). Note that all of them show
topological defects at the same locations, which are marked by white dashed
ellipses. To demonstrate the validity of our analysis, we construct artificial
images, including two different unidirectional CDWs, TDB, and topological
defects b, The phase information implemented in the artificial images. We
assume that the phase windings (topological defects) of the unidirectional CDW
in the left domain only appear on the right side of the image. c, The
artifical STM image is constructed with atomic corrugations, two different
unidirectional CDWs, and their phase information. Each CDW is gradually
attenuated across the TDB marked by a black arrow. The insets are zoomed in
images for the regions marked by colored boxes, including topological defects.
d, The phase map is acquired by analyzing the artificial STM image shown in c
with the unidirectional CDW in the left domain as a reference signal. Note
that the acquired phase map clearly identifies the topological defects that
can be obscured in the STM image as shown in c. Figure 11: Determination of
the anisotropy parameter. a, The first derivative of the STM topographic image
of the TDB shown in Figure 4(a). b and c, The amplitude ($A_{\rm
c(a)}(\mathbf{r})$) maps for each unidirectional CDW. d, The anisotropy
parameter determined by calculating the normalized difference,
$\rm{AP}=\left|\frac{A_{\rm c}-A_{\rm a}}{A_{\rm c}+A_{\rm a}}\right|$, of two
amplitude maps shown in b and c. e, In-gap spectral weight map shown in Figure
4(d). f, Correlation between the anisotropy parameter and the in-gap spectral
weight yielding the correlation factor of -0.77. g, We note that the energy
ranges for the determination of the in-gap spectral weight do not result in a
significant difference.


# Trapped mode control in metasurfaces composed of particles with the form
birefringence property

Anton S. Kupriianov1 Volodymyr I. Fesenko2 Andrey B. Evlyukhin3
<EMAIL_ADDRESS>Wei Han1 Vladimir R. Tuz4,5<EMAIL_ADDRESS>1College of Physics, International Center of Future Science, Jilin University,
2699 Qianjin St., Changchun 130012, China  2Institute of Radio Astronomy of
National Academy of Sciences of Ukraine, 4, Mystetstv St., Kharkiv 61002,
Ukraine  3Institute of Quantum Optics, Leibniz Universität Hannover, 30167
Hannover, Germany  4State Key Laboratory of Integrated Optoelectronics,
College of Electronic Science and Engineering, International Center of Future
Science, Jilin University, 2699 Qianjin Street, Changchun, 130012, China
5School of Radiophysics, Biomedical Electronics and Computer Systems, V. N.
Karazin Kharkiv National University, 4, Svobody Square, Kharkiv 61022, Ukraine

###### Abstract

Progress in developing advanced photonic devices relies on introducing new
materials, discovered physical principles, and optimal designs when
constructing their components. Optical systems operating on the principles of
excitation of extremely high-quality factor trapped modes (also known as the
bound states in the continuum, BICs) are of great interest since they allow
the implementation of laser and sensor devices with outstanding
characteristics. In this paper, we discuss how one can utilize the anisotropic
properties of novel materials (transition metal dichalcogenides, TMDs),
particularly, the bulk molybdenum disulfide (MoS2), to realize the excitation
of trapped modes in dielectric metasurfaces. The bulk MoS2 is a thin-film
structure in which the light wave behaves the same way as that in the uniaxial
anisotropic material with the form birefringence property. Our metasurface is
composed of an array of disk-shaped nanoparticles (resonators) made of the
MoS2 material under the assumption that the anisotropy axis of MoS2 can be
tilted to the rotation axis of the disks. We perform a detailed analysis of
eigenwaves and scattering properties of such anisotropic resonators as well as
the spectral features of the metasurface revealing dependence of the
excitation conditions of the trapped mode on the anisotropy axis orientation
of the MoS2 material used.

††preprint: APS/123-QED

##  I Introduction

Thanks to astonishing advances in technology, electrically thin layers of
materials can be fabricated with subwavelength periodic patterns to create
planar metamaterials (also known as metasurfaces). Such two-dimensional
structures make it possible to control the phase, amplitude, frequency,
angular momentum, and polarization state of electromagnetic waves in a
predetermined way [1, 2]. Due to abrupt phase changes on device surfaces
rather than phases of light propagating inside the systems, this type of
metamaterials can be ultra-thin and flat, which is highly desirable in
integrated optics and photonics. The main task in this field is to design
metasurfaces that provide the most efficient interaction of light with
nanostructured matter. With the use of novel materials, like transition metal
dichalcogenides (TMDs) [3, 4, 5], one can also expect to make such
metasurfaces manageable [6].

Atomically thin TMDs of the form MX2 (M = Mo, W; X = S, Se, Te) are a new
class of semiconductor materials proposed for use in both fundamental physics
exploration in two-dimensional systems and device applications. Among the
TMDs, MoS2 is an indirect bandgap semiconductor with large bandgap energy
which can be transformed into a direct semiconductor in the monolayer limit
[7, 8]. A bulk MoS2 is considered as a layer-by-layer stack of MoS2 monolayers
(thin films), where strong interactions only exist within the basal plane of
covalent bonds. Therefore, the bulk MoS2 is an anisotropic crystal, which is
macroscopically characterized with a tensor-valued permittivity possessing in-
plane $\varepsilon_{\parallel}$ and out-of-plane $\varepsilon_{\bot}$ terms
[9, 10], resembling the form birefringence property [11, 12]. In particular,
these extraordinary properties of MoS2 have been utilized in constructing
anisotropic resonators for lasing systems [13], plasmon-assisted [14, 15] and
dielectric metasurfaces [16, 17, 18, 19, 20].

The optical features of dielectric metasurfaces depend on the excitation
conditions of corresponding Mie-type modes [21, 22] arising in resonators. In
particular, in the spectra of dielectric metasurfaces, there are two basic
resonant states which correspond to the excitation of the magnetic dipole and
electric dipole modes whose position on the frequency scale depends on the
aspect ratio of the particles [23]. Moreover, when metasurfaces containing
anisotropic particles are considered, an additional mode splitting may occur
[24]. Performing tuning mode conditions by changing the aspect ratio of
dielectric nanoparticles allows one to realize flexible manipulation by
resonant states [23, 25] of the nanostructure, and implement, for example,
Huygens metasurfaces [26] where destructive interference between the electric
and magnetic dipole modes of comparable strength suppresses backward
scattering making such metasurfaces to be fully transparent.

Among different designs of dielectric metasurfaces, we further distinguish a
particular class of structures that allow obtaining the strongest resonant
response due to the excitation of so-called trapped modes [27] (recently such
modes are referred to as bound states in the continuum, BICs [28]). In
metasurfaces, such modes are related to non-radiating (dark) eigenwaves which
cannot couple with the field of incident radiation (continuum). The origin and
feature of a particular trapped mode depend on the symmetry of a system and
its spatial, geometric, and material parameters. In particular, a trapped mode
can be excited in metasurfaces if their constitutive particles (resonators)
are in some way perturbed [29]. In this case, a perturbation transforms
inherently dark mode to a weakly radiative one when the spatial symmetry of
the translation unit cell of the metasurface is broken. The cluster-based unit
cell designs [30, 31] as well as the use of specific irradiation conditions
(oblique incidence [32], near-field sources [33], etc.) are other possible
ways to realize the excitation of trapped modes in metasurfaces.

From the point of view of electromagnetic theory, breaking the spatial
symmetry of resonators forming a metasurface leads to the appearance of
artificial anisotropy of their averaged constitutive parameters resulting in
their anisotropic or bianisotropic electromagnetic response [34]. The symmetry
breaking enables the trapped mode coupling with the continuum transforming the
dark mode into a quasi-trapped one. Despite the leakage of energy from the
quasi-trapped mode due to such artificial anisotropy, the concentration of
electromagnetic energy in the resonators can be maintained at a very high
level, if the introduced perturbation is small [28, 35]. Nevertheless,
technologically, producing nanoscale resonators with the required small
perturbation is a rather complicated task. Evidently, this issue can be solved
with the use of anisotropic substances, in particular, the MoS2 material.

Although, generally, the mechanism of excitation of trapped modes in
dielectric metasurfaces is well described [36, 37, 34], the effect of
anisotropy in nanoparticles on such states, as far as we know, has not yet
been thoroughly studied. Therefore, in this paper, we combine the anisotropic
properties of bulk MoS2 considered for particle fabrication and the geometric
parameters of these particles to construct a metasurface and study the
influence of the birefringence property on the resonant conditions of the
trapped mode excitation. In particular, in our metasurface, we utilize disk-
shaped resonators due to their amenability to modern lithography techniques
used for producing metasurfaces operating in the near-infrared range. All our
numerical simulations are targeted on the third telecommunication window [38].
In contrast to the results published earlier for the BIC-assisted metasurfaces
composed of resonators made of TMDs [39, 40, 41, 42], here we consider
resonators without any perturbation of their geometric symmetry when the
metasurface is under the normal irradiation conditions. Our approach is based
on the analysis of eigenwaves of such resonators and their electric and
magnetic dipole responses revealing the dependence of the excitation
conditions of the trapped mode on the anisotropy axis orientation in the bulk
MoS2 material.

##  II Eigenwaves analysis

Figure 1: Schematic illustration of the bulk MoS2 material (thin-film
structure) and sketch view of the metasurface made of MoS2 anisotropic disk-
shaped nanoparticles, the incident wave conditions, and geometric parameters
of the metasurface unit cell. The $c$-axis tilt angle $\varphi$ of the MoS2
anisotropic material is measured from the $z$-axis.

In what follows, we consider a metasurface composed of MoS2 anisotropic disk-
shaped nanoparticles arranged in a two-dimensional array. The translation unit
cell of this array is a square with the side size $p$. The radius and height
of the disks are $a$ and $h$, respectively. The bulk MoS2 material is a thin-
film structure, which is characterized by tensor-valued permittivity
$\hat{\epsilon}$. The light wave in such a material behaves the same way as
that in the uniaxial anisotropic crystal with the form birefringence property
[11, 12]. In particular, it is a nonmagnetic ($\mu$ = $\mu_{0}$) uniaxial
anisotropic crystal whose anisotropy axis ($c$-axis) is directed perpendicular
to the plane of layers forming the MoS2 thin-film structure [9, 43]. For the
sake of definiteness, we assume that the $c$-axis of the crystal is tilted in
the $y$-$z$ plane and can change its orientation on the angle $\varphi$ with
respect to the $z$-axis of the chosen Cartesian coordinate frame, where the
$z$-axis is also the rotation (axial symmetry) axis of the disks, as shown in
Fig. 1. In the operating frequency range under study, material losses in the
bulk MoS2 material are negligibly small (see Appendix A) and, therefore, we
exclude them from our consideration. Without loss of generality, we assume
that the particles are located in a homogeneous surrounding space (air). The
trapped mode manifestation is related to the in-plane symmetry breaking of the
metasurface unit cell. The substrate presence does not break the in-plane
symmetry and has no significant influence on the effects discussed. Therefore,
in order not to overload our study, we consider that the substrate is made of
an air-like material $\varepsilon_{s}\approx 1$. Accounting for the substrate
can be performed later in the engineering implementation of the metasurface by
adjusting the geometric parameters of the resonators used.

Figure 2: (a) Patterns of the electric and magnetic fields plotted for the
$\textsc{TE}_{01\delta}$, $\textsc{EH}_{11\delta}^{x,y}$, and
$\textsc{HE}_{11\delta}^{x,y}$ eigenwaves of the MoS2 metasurface at
$\varphi=0^{\circ}$, where the red and dark blue arrows represent the flow of
electric and magnetic fields while brown and gray arrows represent moments of
the electric dipole $\bf p$ and magnetic dipole $\bf m$, respectively. (b)
Eigenfrequencies ($\omega^{\prime}/2\pi$) and (c) quality factors
($\omega^{\prime}/2\omega^{\prime\prime}$) of eigenwaves as functions of the
$c$-axis tilt angle $\varphi$. Parameters of the metasurface are: $a=258$ nm,
$h=275$ nm, and $p$ = 1500 nm.

Initially, we perform an eigenwave analysis of the given metasurface. For our
simulations, we use the RF module of the commercial COMSOL Multiphysics
finite-element electromagnetic solver. In the solver, we construct one square
unit cell of the metasurface imposing the Floquet-periodic boundary conditions
on its sides to simulate the infinitely expanded two-dimensional array of
nanoparticles. From the solutions found, we select those that correspond to
the excitation of the lowest-order eigenwaves only. We classify these
eigenwaves by considering the electromagnetic field configuration inside the
resonator when the $c$-axis tilt angle $\varphi$ is zero and permittivity of
the disk acquires the diagonal form
$\hat{\epsilon}|_{\varphi=0^{\circ}}=\\{\varepsilon_{\parallel},\varepsilon_{\parallel},\varepsilon_{\bot}\\}$.
For the chosen geometric parameters of the resonators, the eigenwaves arise as
the transverse electric TE01δ, and hybrid EH11δ and HE11δ modes, sequentially
arranged with increasing frequency (this classification is based on the mode
nomenclature of cylindrical dielectric waveguides, see Ref. [44]). In the
subscripts of the mode abbreviations, the first index denotes the azimuthal
variation of the fields, the second index denotes the order of variation of
the field along the radial direction, and the third index denotes the order of
variation of fields along the $z$-direction. Since we do not fix the exact
number of field variations along the $z$-axis, we substitute the third index
with the letter $\delta$. In general, the characteristics of eigenwaves under
study inherit the properties of the corresponding modes of an anisotropic
nanowire derived in Appendix B.

Moreover, as known from the theory of dielectric resonators [45], the TE01δ
mode of an isolated cylindrical resonator radiates like a magnetic dipole $\bf
m$ oriented along its rotation axis (out-of-plane), whereas the EH11δ and
HE11δ modes radiate like an electric dipole $\bf p$ and magnetic dipole $\bf
m$, respectively, oriented along the transverse direction (in-plane). Thus,
each hybrid mode is complemented by its own degenerated state, depending on
the direction in which the corresponding dipole moment is oriented, along
either the $x$-axis or the $y$-axis in the chosen Cartesian coordinate frame.
We distinguish these degenerated hybrid modes by adding the corresponding
superscript $x$ or $y$ to the mode abbreviations. The appearance of all
eigenwaves of our interest at the initial stage $\varphi=0$ are collected in
Fig. 2(a).

Among these eigenwaves, the field for the TE01δ mode is axisymmetric and thus
has no azimuthal variation, whereas the fields of the hybrid modes
EH${}_{11\delta}^{x,y}$ and HE${}_{11\delta}^{x,y}$ are azimuthally dependent.
This property of the TE01δ mode leads to the appearance of conditions
corresponding to the existence of the trapped mode in our lossless metasurface
signifying that the corresponding eigenfrequency ($\omega=\omega^{\prime}+{\rm
i}\omega^{\prime\prime}$) is a purely real quantity
($\omega^{\prime\prime}=0$).

Further in our study, we vary the value of the tilt angle $\varphi$. When
$\varphi\in(0^{\circ},90^{\circ})$, there are off-diagonal components in the
tensor $\hat{\epsilon}$, whereas when $\varphi=90^{\circ}$, the $c$-axis is
directed along the $y$-axis, and permittivity of the disk becomes diagonal
$\hat{\epsilon}|_{\varphi=90^{\circ}}=\\{\varepsilon_{\parallel},\varepsilon_{\bot},\varepsilon_{\parallel}\\}$
(see Appendix B). We track the change in the values of the real part of the
eigenfrequency ($\omega^{\prime}$) and quality factor
($\omega^{\prime}/2\omega^{\prime\prime}$) of the selected eigenwaves. The
corresponding results of our calculations are presented in Figs. 2(b) and
2(c), respectively. One can see that the hybrid EH${}_{11\delta}^{x,y}$ and
HE${}_{11\delta}^{x,y}$ modes lose their degeneracy and split apart with a
change in the tilt angle $\varphi$ of the anisotropy axis, while the TE01δ
mode does not. Moreover, the quality factor of the TE01δ mode changes from
infinity to finite values undergoing a transformation from the trapped to a
radiative state. Hereinafter, our task is to reveal the manifestation of these
eigenwaves in the transmitted spectra of our metasurface.

##  III Spectral features

Figure 3: Transmission coefficient magnitude as a function of the frequency
$f=\omega/2\pi$ and $c$-axis tilt angle $\varphi$ for the (a) $x$-polarized
and (b) $y$-polarized incident waves. (c) Patterns of the electric and
magnetic fields plotted at the corresponding resonant frequencies for the
$\textsc{TE}_{01\delta}$, $\textsc{EH}_{11\delta}^{x}$, and
$\textsc{HE}_{11\delta}^{y}$ modes of the metasurface at three different
values of $\varphi$ and the $x$-polarized wave. The red and dark blue arrows
represent the flow of electric and magnetic fields while the brown and gray
arrows represent the electric $\bf p$ and magnetic $\bf m$ dipole moments,
respectively. Parameters of the metasurface are the same as in Fig. 2.

In this section, we reveal specific spectral features of the MoS2 metasurface
caused by the material anisotropy. For this study we consider that the
metasurface is illuminated by a plane electromagnetic wave under the normal
incidence conditions ($\textbf{k}=\\{0,0,-k_{z}\\}$) with the electric field
polarized along either the $x$-axis ($\textbf{E}=\\{E_{x},0,0\\}$,
$x$-polarization) or the $y$-axis ($\textbf{E}=\\{0,E_{y},0\\}$,
$y$-polarization). In the framework of the RF module of the COMSOL
Multiphysics solver, we dispose of the radiating and receiving ports above and
below the metasurface, respectively. The port boundary conditions are placed
on the interior boundaries of the perfectly matched layers, adjacent to the
air domains. Perfectly matched layers on the top and bottom of the unit cell
absorb the excited wave from the radiating port and prevent unwanted re-
reflection inside the computational domain. The port boundary conditions
automatically determine the reflection and transmission characteristics of the
metasurface in terms of S-parameters. After simulation, we retrieve values of
the transmission coefficient ($|T|^{2}=|S_{21}|^{2}$) as a function of the
frequency and $c$-axis tilt angle $\varphi$ for two orthogonal polarizations
of the incident wave. The results of our calculations are summarized in Figs.
3(a) and 3(b) for the $x$-polarized and $y$-polarized waves, respectively.

When $\varphi=0^{\circ}$, the metasurface is polarization-independent. In the
frequency band of interest, there are two resonant states related to the
excitation of the degenerated hybrid EH11δ and HE11δ modes. The identified
resonances appear due to an existing electromagnetic coupling between the
linearly polarized incident wave and transversely located electric and
magnetic dipoles inherent to the corresponding eigenwaves. Since the TE01δ
mode is axially symmetric with the corresponding magnetic dipole oriented
orthogonally to the metasurface plane, the field of the incident wave does not
interact with this mode and there is no corresponding resonance in the spectra
of the transmitted wave.

Variation of the $c$-axis tilt angle $\varphi$ breaks the in-plane symmetry of
the structure, resulting that the metasurface becomes to be polarization-
dependent. Although the degeneracy for hybrid modes is lifted, their resonant
positions on the frequency scale change a little compared to those of the
degenerated modes that existed in the $\varphi=0^{\circ}$ case. However, since
the rotation of the anisotropy $c$-axis is maintained in the $y$-$z$ plane,
noticeable changes are observed in the spectra of the $x$-polarized wave. At a
certain value of the $c$-axis tilt angle ($\varphi\approx 50^{\circ}$), the
resonant frequencies of the EH${}_{11\delta}^{x}$ and HE${}_{11\delta}^{y}$
modes coincide, and with a further increase in $\varphi$, the resonances
rearrange their order on the frequency scale. This mode-crossing effect is
widely discussed in the context of the implementation of Huygens metasurfaces
[26]. Moreover, the quality factor of the EH${}_{11\delta}^{y}$ and
HE${}_{11\delta}^{y}$ resonances increases, since the electromagnetic coupling
of the corresponding either electric $\bf p$ or magnetic $\bf m$ dipole with
the field of the incident wave decreases.

Aside from the resonances related to the above-considered hybrid EH11δ and
HE11δ modes which possess a good electromagnetic coupling with the field of
the normally incident linearly polarized wave, in the frequency band of
interest, an additional resonance arises in the transmitted spectra of the
$x$-polarized waves as soon as $\varphi$ is nonzero. This resonance is related
to the manifestation of the TE01δ mode. It appears as an alone-standing lowest
frequency resonance whose quality factor decreases with increasing $\varphi$.
This resonance acquires a peak-and-trough (Fano) profile as is typical for the
trapped modes excitation [27, 36]. While the dip in curves corresponds to the
maximum of reflection, the peak corresponds to the maximum of transmission.
These extremes approach 0 and 1, respectively, since the material losses in
the given metasurface are considered to be negligibly small.

To reveal the origin of changes in the emerging resonances, we plot the
patterns of the electric and magnetic fields within the unit cell in three
basic projections and show the orientation of the dipole moments $\bf m$ and
$\bf p$ at the resonator center. These calculations are made for three
different values of the $c$-axis tilt angle $\varphi$. The results of our
calculations are collected in Fig. 3(c). One can see that variation in the
material anisotropy leads to a change in the slope of the dipole moment of the
corresponding mode, thereby changing the degree of its electromagnetic
coupling with the field of the incident wave. This mechanism is fully
consistent with the behavior of trapped modes in resonators with perturbed
unit cells [36, 46, 47] and the corresponding theory [37, 34, 40] can be
adapted in this case. It is based on the scattering characteristics of a
single anisotropic nanoparticle which we discussed in Appendix C.

##  IV Polarization features

Figure 4: Polarization features of the transmitted field: (a) transmission
coefficient magnitude of the circularly polarized waves, ellipticity, and
rotation angle (polarization azimuth) as functions of the frequency
$f=\omega/2\pi$ for two particular values of the $c$-axis tilt angle
$\varphi$; (b) polarization ellipse and (c) transmission coefficient magnitude
mapped on the surface of the Poincaré sphere related to the polarization
states of the incident wave at the resonant frequency of the trapped mode
excitation. The point at the north (south) pole of the Poincaré sphere
represents the left-handed (right-handed) circularly polarized field, and the
point on the equator represents the $45^{\circ}$ linearly polarized field.
Parameters of the metasurface are the same as in Fig. 2.

Since the metasurface is composed of anisotropic particles, the polarization
features of the metasurface at the resonant frequency of the trapped mode
excitation should be elucidated. To this end, we study the transmission
characteristics of our metasurface with respect to all polarization states of
the incident wave for several values of the $c$-axis tilt angle $\varphi$.

Figure 4(a) shows the corresponding frequency dependences of the transmission
coefficient magnitude of circularly polarized waves, ellipticity angle
($\eta$), and polarization azimuth ($\theta$) for the transmitted field.
According to the definition of the Stokes parameters, we introduce the
ellipticity $\eta$ so that the field is linearly polarized when $\eta=0$, and
$\eta=-\pi/4$ for the left-handed circularly polarized (LCP) wave and $+\pi/4$
for the right-handed circularly polarized (RCP) wave (note that in the latter
cases the preferential azimuthal angle of the polarization ellipse $\theta$
becomes undefined). In all other cases ($0<|\eta|<\pi/4$), the field is
elliptically polarized. In the considered frequency band, the transmission
coefficient magnitude of the given metasurface is the same for the LCP and RCP
incident waves, whereas the transmitted field undergoes the rotation of its
polarization ellipse at the resonant frequencies where it can change its
polarization state. However, at the resonant frequency of the trapped mode
excitation, the transmitted field remains to be almost linearly polarized
regardless of the value of the $c$-axis tilt angle $\varphi$. The
corresponding polarization ellipse of the transmitted field for the trapped
mode resonance is presented in Fig. 4(b), whose frequency is depicted in Fig.
4(a) by a vertical dashed line.

This polarization feature of the trapped mode is further confirmed by our
mapping of the transmission coefficient magnitude of the metasurface for
different polarization states onto the surface of the full Poincaré sphere as
shown in Fig. 4(c). In our calculations, the electric field vector ${\bf E}$
of the incident wave is defined by components $E_{x}=[0.5(P+1)]^{1/2}$ and
$E_{y}=[0.5(P-1)]^{1/2}\exp(\rm i\beta)$, where $P\in[-1,1]$ and
$\beta\in[-\pi,\pi]$ [48, 49, 50]. The axes of the Poincaré sphere are labeled
in terms of the Stokes parameters. Linearly polarized states are arranged
along the sphere’s equator, while right-handed (RCP) and left-handed (LCP)
circularly polarized waves are located at their north and south poles,
respectively. All other points of the sphere depict elliptically polarized
states of the incident wave. From the uniform color gradient of the sphere, it
is easy to conclude that at the resonant frequency of the trapped mode
excitation, there are the total reflection and total transmission of $x$\- and
$y$-polarized incident waves, respectively, while circularly polarized waves
pass through the metasurface with half the amplitude of the incident wave.

##  V Conclusions

We have revealed the role of material anisotropy in the mechanism of the
trapped mode excitation in a metasurface composed of MoS2 disk-shaped
nanoparticles. Our study is based on the analysis of the properties of the
three lowest-order modes of cylindrical resonators forming the metasurface.
They are the axially symmetric TE01δ mode as well as the non-symmetric hybrid
EH11δ and HE11δ modes.

When the anisotropy axis of the MoS2 material used for the disk-shaped
particles fabrication coincides with the axis of their rotational symmetry,
the hybrid EH11δ and HE11δ modes degenerate, whereas the TE01δ mode behaves
like a trapped mode in the metasurface, possessing a purely real
eigenfrequency. When the anisotropy axis is tilted, the hybrid EH11δ and HE11δ
modes are split, and the TE01δ mode transforms into a radiative one, whose
quality factor decreases with increasing the anisotropy tilt angle. It is
shown how a change in the anisotropy tilt angle affects the orientation of
vectors of the electric and magnetic dipole moments of the corresponding
modes.

While we have considered a specific MoS2 material, the trapped mode excitation
mechanism in metasurfaces composed of dielectric nanoparticles that we have
demonstrated is general for both natural and artificial materials with
uniaxial anisotropy. By performing control ability over the high-quality
factor trapped modes in MoS2 metasurfaces by manipulating material anisotropy,
we propose a platform for many important applications in light-matter
interaction, nonlinear photonics, quantum optics, and spectroscopy.

###### Acknowledgements.

V.R.T. is grateful for the hospitality and support from Jilin University,
China. A.B.E. thanks funding support from the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) under Germany’s Excellence Strategy within
the Cluster of Excellence PhoenixD (EXC 2122, Project ID 390833453).

##  Appendix A MoS2 relative permittivity

Figure 5: Relative permittivity of MoS2 extracted from experimental data of
Ref. [10].

Dispersion characteristics of components of the tensor of the relative
permittivity ($\hat{\epsilon}=\epsilon^{\prime}+{\rm
i}\epsilon^{\prime\prime}$) of the MoS2 material are shown in Fig. 5. Since in
our work, we focus on the spectral range within the third telecommunication
window ($\lambda\in[1500,1600]~{}\rm{nm}$), where negligibly weak dispersion
exists (see the spectral region marked light blue), we adopt
$\varepsilon_{\parallel}^{\prime}=16.6$, $\varepsilon_{\perp}^{\prime}=6.3$,
$\varepsilon_{\parallel}^{\prime\prime}=\varepsilon_{\perp}^{\prime\prime}=0.0$
in all considered region.

##  Appendix B Lowest-order modes of an anisotropic nanowire

For reference, here we present dispersion features of three lowest-order modes
(TE01, EH11, and HE11) of an anisotropic nanowire. The nanowire is made of a
nonmagnetic material (MoS2) with the tensor-valued relative permittivity
$\hat{\epsilon}$. The nanowire is disposed of in a surrounding medium with
permittivity $\varepsilon_{0}$. The radius of the nanowire is $a$. The
symmetry axis of the nanowire coincides with the $z$-axis of the Cartesian
coordinate system [see the inset in Fig. 6(a)].

Figure 6: (a)-(c) Dispersion curves of the axially symmetric
$\textsc{TE}_{01}$ and non-symmetric hybrid $\textsc{HE}_{11}^{x,y}$,
$\textsc{EH}_{11}^{x,y}$ modes of the cylindrical MoS2 nanowire at different
values of the $c$-axis tilt angle $\varphi$: (a) $\varphi=0^{\circ}$, (b)
$\varphi=45^{\circ}$, and (c) $\varphi=90^{\circ}$. The circles in panel (a)
correspond to the analytical solution given by Eq. (3). (d) The effective mode
index $\tilde{\beta}$ of the $\textsc{TE}_{01}$, $\textsc{HE}^{x,y}_{11}$ and
$\textsc{EH}^{x,y}_{11}$ modes versus the $c$-axis tilt angle $\varphi$ at the
fixed frequency $f=200.0$ THz. (e) Simulated electric field intensities
$|\textbf{E}|$ of the nanowire modes corresponding to panels (a), (b), and
(c). A schematic of the MoS2 nanowire is shown in the inset of panel (a). The
red arrow represents the $c$-axis direction in the nanowire core. The nanowire
radius is $a=258$ nm, and components of the tensor of relative permittivity
are: $\varepsilon_{\parallel}=16.6$ and $\varepsilon_{\bot}=6.24$.

The optic $c$-axis of the uniaxial crystalline MoS2 core of the nanowire is
tilted in the $y$-$z$-plane and makes an angle $\varphi$ with the $z$-axis.
The full relative permittivity tensor of the nanowire core
$\hat{\epsilon}(\varphi)$ can be evaluated by using the rotation matrix

$\hat{R}=\left({\begin{matrix}1&0&0\cr 0&{\cos{\varphi}}&{-\sin{\varphi}}\cr
0&{\sin{\varphi}}&{\cos{\varphi}}\cr\end{matrix}}\right)$ (1)

as follows

$\hat{\epsilon}(\varphi)=\hat{R}\cdot\hat{\epsilon}|_{\varphi=0}\cdot\hat{R}^{T}=\left({\begin{matrix}\varepsilon_{xx}&0&0\cr
0&{\varepsilon_{yy}\cos^{2}{\varphi}+\varepsilon_{zz}\sin^{2}{\varphi}}&{\varepsilon_{yy}\sin{\varphi}\cos{\varphi}-\varepsilon_{zz}\sin{\varphi}\cos{\varphi}}\cr
0&{\varepsilon_{yy}\sin{\varphi}\cos{\varphi}-\varepsilon_{zz}\sin{\varphi}\cos{\varphi}}&{\varepsilon_{zz}\cos^{2}{\varphi}+\varepsilon_{yy}\sin^{2}{\varphi}}\cr\end{matrix}}\right),$
(2)

where
$\hat{\epsilon}|_{\varphi=0^{\circ}}=\\{\varepsilon_{xx},\varepsilon_{yy},\varepsilon_{zz}\\}=\\{\varepsilon_{\parallel},\varepsilon_{\parallel},\varepsilon_{\bot}\\}$
and superscript $T$ denotes the matrix transposition operation.

As follows from Eq. (2), the permittivity of the nanowire can be defined in
the diagonal form as
$\hat{\epsilon}|_{\varphi=0^{\circ}}=\\{\varepsilon_{\parallel},\varepsilon_{\parallel},\varepsilon_{\bot}\\}$
and
$\hat{\epsilon}|_{\varphi=90^{\circ}}=\\{\varepsilon_{\parallel},\varepsilon_{\bot},\varepsilon_{\parallel}\\}$
in two special cases when the $c$-axis is directed along the $z$-axis and the
$y$-axis, respectively. For the rest range
$\varphi\in\left(0^{\circ},90^{\circ}\right)$, the tensor
$\hat{\epsilon}(\varphi)$ has nonzero off-diagonal components.

When the directions of the $c$-axis and the $z$-axis coincide, there is an
analytical solution for the nanowire modes [51]. Generally, such an
anisotropic nanowire supports an infinite set of the axially symmetric TEmn
and TMmn modes as well as the non-symmetric hybrid EHmn and HEmn modes, whose
dispersion equation can be expressed as [52]:

$\begin{split}&\left[\frac{\varepsilon_{0}}{\gamma
a}\frac{K^{\prime}_{m}(\gamma a)}{K_{m}(\gamma
a)}+\frac{\sqrt{\varepsilon_{\parallel}\varepsilon_{\bot}}}{\kappa
a}\frac{J^{\prime}_{m}(p\kappa a)}{J_{m}(p\kappa a)}\right]\\\
&\times\left[\frac{1}{\gamma a}\frac{K^{\prime}_{m}(\gamma a)}{K_{m}(\gamma
a)}+\frac{1}{\kappa a}\frac{J^{\prime}_{m}(\kappa a)}{J_{m}(\kappa
a)}\right]\\\ &=\left(\frac{m\beta}{k_{0}}\right)^{2}\left[\frac{1}{(\gamma
a)^{2}}+\frac{1}{(\kappa a)^{2}}\right]^{2},\end{split}$ (3)

where, $k_{0}$ is the wavenumber in free space, $m$ and $n$ are the integer
values known as the azimuthal and radial indices, respectively, $\beta$ is the
longitudinal propagation constant,
$p=\sqrt{\varepsilon_{\bot}/\varepsilon_{\parallel}}$ is the anisotropy
parameter, $\kappa=\sqrt{k_{0}^{2}\varepsilon_{\parallel}-\beta^{2}}$ and
$\gamma=\sqrt{\beta^{2}-k_{0}^{2}\varepsilon_{0}}$ are the transverse
propagation constants inside and outside the nanowire, respectively,
$J_{m}(\cdot)$ and $K_{m}(\cdot)$ are the Bessel function of the first kind
and the modified Bessel function of the second kind, and
$J^{\prime}_{m}(\cdot)$ and $K^{\prime}_{m}(\cdot)$ are their derivatives with
respect to the function argument.

Figure 7: Contributions of the electric dipole (ED), magnetic dipole (MD), and
total scattering cross-section of a single anisotropic MoS2 nanoparticle
(disk) for three different orientations of the $c$-axis of anisotropy with
respect to the rotation axis of the disk, where (a) $\varphi=0^{\circ}$, (b)
$\varphi=45^{\circ}$, and (c) $\varphi=90^{\circ}$, and the solid and dashed
lines correspond to the frontal disk irradiation by the $x$-polarized and
$y$-polarized wave, respectively. Geometric parameters of the disk are:
$a=258$ nm, $h=275$ nm, and components of the tensor of relative permittivity
are: $\varepsilon_{\parallel}=16.6$ and $\varepsilon_{\bot}=6.24$.

A formal substitution of $m=0$ into Eq. (3) yields us two separate dispersion
equations related to the transverse magnetic TM0n and transverse electric TE0n
modes [see, the first and second multiplier on the left side of the dispersion
equation (3), respectively], similarly to the case of a nanowire with an
isotropic core [53].

In the general case of a uniaxial crystal core when the $c$-axis does not
coincide with the $z$-axis, the propagation constants and field components can
be determined numerically [54, 55, 56]. The calculated effective mode index
$\tilde{\beta}=\beta/k_{0}$ and the electric field patterns of the lowest-
order $\textsc{TE}_{01}$, $\textsc{EH}^{x,y}_{11}$, and
$\textsc{HE}_{11}^{x,y}$ modes propagated in the MoS2 nanowire are presented
in Fig. 6 for different values of the $c$-axis tilt angle $\varphi$. In
particular, when the $c$-axis is directed along the $z$-axis, the results of
our simulation are checked against those obtained from the analytical solution
given by Eq. (3).

From Figs. 6(a)-6(c), one can conclude that when
$\varphi\in\left[0^{\circ},90^{\circ}\right]$, the dispersion characteristic
of the $\textsc{TE}_{01}$ mode does not depend on the transverse polarization,
that is, this mode does not split into the $x$-polarized
($\textsc{TE}_{01}^{x}$) and $y$-polarized ($\textsc{TE}_{01}^{y}$) components
(see, also Refs. [57, 58]). Nevertheless, when $\varphi\neq 0^{\circ}$, the
axial symmetry of the electric field of the $\textsc{TE}_{01}$ mode becomes
violated, as shown in Fig. 6(e).

In turn, the dispersion characteristics of hybrid modes are more complex.
Figure 6(a) shows that at $\varphi=0^{\circ}$, no modal birefringence
phenomenon exists, it means that the hybrid $\textsc{HE}^{x,y}_{11}$ and
$\textsc{EH}^{x,y}_{11}$ modes are degenerated as shown in Fig. 6(e), which is
typical for uniaxial nanowires [57, 58].

In the range $\varphi\in(0^{\circ},90^{\circ}]$, the dispersion
characteristics of hybrid $\textsc{EH}^{x,y}_{11}$ and
$\textsc{HE}^{x,y}_{11}$ modes depend on the transverse polarization
orientations, which leads to degeneracy lifting as presented in Fig. 6(e).
Thus, the different polarizations of $\textsc{HE}_{11}^{x}$
($\textsc{EH}_{11}^{x}$) and $\textsc{HE}_{11}^{y}$ ($\textsc{EH}_{11}^{y}$)
modes have different dispersion relations and propagate like different modes,
as shown in Figs. 6(b) and 6(c). For the chosen parameters of anisotropy, the
effect of the $c$-axis orientation is much stronger for the $y$-polarized
modes compared to the $x$-polarized ones due to the fact that
$\varepsilon_{yy}(\varphi)\in\left[\varepsilon_{\bot},\varepsilon_{\parallel}\right]$
while $\varepsilon_{xx}(\varphi)=\varepsilon_{\parallel}$ [see, for clarity,
Eq. (2)]. A similar feature has been reported previously for an anisotropic
rectangular waveguide [55]. Besides, Fig. 6(d) indicates that difference
between the $x$-polarized and $y$-polarized modes increases as the $c$-axis
tilt angle $\varphi$ rises, and the difference between $\textsc{EH}_{11}^{x}$
and $\textsc{EH}_{11}^{y}$ modes is smaller than the difference between the
$\textsc{HE}_{11}^{x}$ and $\textsc{HE}_{11}^{y}$ modes.

##  Appendix C Scattering cross-section of an anisotropic nanoparticle

Figure 8: Variation of the magnitude of the $z$ component of the electric $\bf
p$ and magnetic $\bf m$ dipole moments for several values of the $c$-axis tilt
angle $\varphi$ for the frontal irradiation of the disk by the (a)
$x$-polarized and (b) $y$-polarized wave. Parameters of the disk are the same
as in Fig. 7.

For numerical simulations of the single particle scattering cross-section and
the particle electric and magnetic dipole moments, we use the discrete dipole
approximation (DDA) and the decomposed discrete dipole approximation (DDDA)
methods [23]. The main idea of the DDA method consists in the replacement of
the scattering object by a cubic lattice of electric point dipoles with known
polarizability tensor $\hat{\alpha}_{p}$. The corresponding dipole moment
${\bf d}_{i}$ induced in each lattice point $i$ (with the radius vector ${\bf
r}_{i}$) is found by solving coupled dipole equations:

${\bf d}^{i}=\hat{\alpha}_{p}{\bf
E}^{i}_{0}+\hat{\alpha}_{p}\frac{k_{0}^{2}}{\varepsilon_{0}}\sum_{j\neq
i}^{N}\hat{G}_{ij}{\bf d}^{j},\quad i=1...N,$ (4)

where ${\bf E}^{i}_{0}$ is the electric field of the incident wave at the
point ${\bf r}_{i}$, $k_{0}$ is the wavenumber in a vacuum, $\varepsilon_{0}$
is the vacuum dielectric constant, $\hat{G}_{ij}=\hat{G}({\bf r}_{i},{\bf
r}_{j})$ is the Green tensor of the system without the particle, $N$ is the
number of dipoles replacing the particle after the discretization procedure.
For a MoS2 particle (disk), when the $c$-axis of anisotropy is directed along
the $z$-axis in the chosen Cartesian coordinate frame ($\varphi=0$), the
permittivity tensor is written as

$\hat{\epsilon}=\left(\begin{array}[]{ccc}\varepsilon_{\parallel}&0&0\\\
0&\varepsilon_{\parallel}&0\\\ 0&0&\varepsilon_{\bot}\end{array}\right),$ (5)

and the polarizability tensor can be taken as

$\hat{\alpha}_{p}\simeq
3\varepsilon_{0}\varepsilon_{s}V_{p}\left(\begin{array}[]{ccc}\frac{\varepsilon_{\parallel}-\varepsilon_{s}}{\varepsilon_{\parallel}+2\varepsilon_{s}}&0&0\\\
0&\frac{\varepsilon_{\parallel}-\varepsilon_{s}}{\varepsilon_{\parallel}+2\varepsilon_{s}}&0\\\
0&0&\frac{\varepsilon_{\bot}-\varepsilon_{s}}{\varepsilon_{\bot}+2\varepsilon_{s}}\end{array}\right),$
(6)

where $V_{p}$ is the volume of the discretization cell, $\varepsilon_{s}$ is
the relative permittivity of the surrounding medium. Here we assume that the
particle is localized in a vacuum or in the air where $\varepsilon_{s}=1$.
After the solution of Eqs. (4) and choosing the origin of the Cartesian
coordinate system in the center of mass of the particle, the scattered
electric field ${\bf E}^{\textrm{scat}}$ in the far-field zone and the
spherical coordinates are written as

$\displaystyle{E}^{\textrm{scat}}_{\phi}(r,\phi,\theta)$ $\displaystyle=$
$\displaystyle\frac{k_{0}^{2}e^{{\rm
i}k_{s}r}}{4\pi\varepsilon_{0}r}\sum_{j=1}^{N}e^{-{\rm i}k_{s}({\bf
n}\cdot{\bf r}_{j})}\left(d^{j}_{y}\cos\phi\right.$
$\displaystyle\left.-d_{x}^{j}\sin\phi\right),$
$\displaystyle{E}^{\textrm{scat}}_{\theta}(r,\phi,\theta)$ $\displaystyle=$
$\displaystyle\frac{k_{0}^{2}e^{{\rm
i}k_{s}r}}{4\pi\varepsilon_{0}r}\sum_{j=1}^{N}e^{-{\rm i}k_{s}({\bf
n}\cdot{\bf r}_{j})}\left(d^{j}_{x}\cos\phi\cos\theta\right.$ (7)
$\displaystyle\left.+d_{y}^{j}\sin\phi\cos\theta-d_{z}^{j}\sin\theta\right),$

where $k_{s}=k_{0}\sqrt{\varepsilon_{s}}$, ${\bf
n}=(\sin\theta\cos\phi,~{}\sin\theta\sin\phi,~{}\cos\theta)$ is the unit
vector directed to the observation point with spherical coordinate
$(r,~{}\phi,~{}\theta)$, $r$ is the distance to the observation point, $\phi$
and $\theta$ are the azimuthal and polar angles, respectively.

The scattering cross-section $\sigma_{\textrm{scat}}$ is calculated as

$\sigma_{\textrm{scat}}=\frac{1}{|{\bf
E}_{0}|^{2}}\int_{\Omega}\left\\{|{E}^{\textrm{scat}}_{\phi}|^{2}+|{E}^{\textrm{scat}}_{\theta}|^{2}\right\\}r^{2}d\Omega,$
(8)

where $\Omega$ and $d\Omega=\sin\theta d\phi d\theta$ is the full solid angle
and its differential.

In the DDDA, a discrete representation of the induced polarization inside the
particle is obtained as

${\bf P}({\bf r})=\sum_{j=1}^{N}{\bf d}^{j}\delta({\bf r}-{\bf r}_{j}),$ (9)

where $\delta({\bf r}-{\bf r}_{j})$ is the Dirac delta function.

Using the relation ${\bf j}=-{\rm i}\omega{\bf P}$, where $\omega$ is the
field angular frequency, between the induced polarization $\bf P$ and the
induced electric current density $\bf j$ and Eq. (9), the multipole moments of
the particle can be calculated by applying the expressions from Ref. [59]. In
this approach the multipole decomposition of the scattering cross-section is

$\displaystyle\sigma_{\textrm{scat}}$
$\displaystyle\simeq\frac{k_{0}^{4}}{6\pi\varepsilon_{0}^{2}|{\bf
E}_{0}|^{2}}|{\bf
p}|^{2}+\frac{k_{0}^{4}\varepsilon_{s}\mu_{0}}{6\pi\varepsilon_{0}|{\bf
E}_{0}|^{2}}|{\bf m}|^{2}$
$\displaystyle+\frac{k_{0}^{6}\varepsilon_{s}}{720\pi\varepsilon_{0}^{2}|{\bf
E}_{0}|^{2}}\sum_{\alpha\beta}|Q_{\alpha\beta}|^{2}+\frac{k_{0}^{6}\varepsilon_{s}^{2}\mu_{0}}{80\pi\varepsilon_{0}|{\bf
E}_{0}|^{2}}\sum_{\alpha\beta}|M_{\alpha\beta}|^{2},$

where $\mu_{0}$ is vacuum permeability, $\bf p$ and $\bf m$ are the electric
and magnetic dipole moments, $\hat{Q}$ and $\hat{M}$ are the electric and
magnetic quadrupole moments of the particle.

Note that from the practical point of view, it is convenient to rotate the
disk under the fixed polarizability tensor (6) in the numerical simulation in
order to change the anisotropic properties of the disk.

The results of our modeling of the cross-section characteristics of a single
particle for several different values of the $c$-axis tilt angle $\varphi$ are
shown in Fig. 7. First, one can see that the resonant response of the disk
depends on both orientations of the MoS2 layers and irradiation conditions.
Second, in the chosen frequency range, the resonant features are associated
only with the resonant excitation of either the electric (ED) or magnetic (MD)
dipole moments of the particle. Importantly, in the range,
$\varphi\in(0,90^{\circ})$, the material anisotropy of the disk leads to the
resonant excitation of the longitudinal (directed along the wave vector $\bf
k$ of the incident wave) components of the ED and MD dipole moments (Fig. 8).
The excitation of $p_{z}$ ($m_{z}$) component is realized for the
$y$-polarized ($x$-polarized) waves, respectively, and spectral positions of
their extremes are completely correlated with the resonant features of the
transmission spectra of the MoS2 metasurface presented in Figs. 3(a) and 3(b),
respectively. This demonstrates that the discussed resonant states of the
anisotropic metasurface are basically influenced by the anisotropy and modes
of an individual dielectric resonator rather than by the electromagnetic
coupling between them.

## References

  * Chen  _et al._ [2016] H.-T. Chen, A. J. Taylor, and N. Yu, A review of metasurfaces: physics and applications, Rep. Prog. Phys. 79, 076401 (2016).
  * Li  _et al._ [2018] A. Li, S. Singh, and D. Sievenpiper, Metasurfaces and their applications, Nanophotonics 7, 989 (2018).
  * Wang  _et al._ [2018] G. Wang, A. Chernikov, M. M. Glazov, T. F. Heinz, X. Marie, T. Amand, and B. Urbaszek, Colloquium: Excitons in atomically thin transition metal dichalcogenides, Rev. Mod. Phys. 90, 021001 (2018).
  * Verre  _et al._ [2019] R. Verre, D. G. Baranov, B. Munkhbat, J. Cuadra, M. Käll, and T. Shegai, Transition metal dichalcogenide nanodisks as high-index dielectric Mie nanoresonators, Nat. Nanotechnol. 14, 679 (2019).
  * Munkhbat  _et al._ [2022] B. Munkhbat, B. Küçüköz, D. G. Baranov, T. J. Antosiewicz, and T. O. Shegai, Nanostructured transition metal dichalcogenide multilayers for advanced nanophotonics, Laser Photonics Rev. , 2200057 (2022).
  * Xiao  _et al._ [2020] S. Xiao, T. Wang, T. Liu, C. Zhou, X. Jiang, and J. Zhang, Active metamaterials and metadevices: a review, J. Phys. D: Appl. Phys. 53, 503002 (2020).
  * Yan  _et al._ [2015] X. Yan, L. Zhu, Y. Zhou, Y. E, L. Wang, and X. Xu, Dielectric property of MoS2 crystal in terahertz and visible regions, Appl. Opt. 54, 6732 (2015).
  * Ermolaev  _et al._ [2020] G. A. Ermolaev, Y. V. Stebunov, A. A. Vyshnevyy, D. E. Tatarkin, D. I. Yakubovsky, S. M. Novikov, D. G. Baranov, T. Shegai, A. Y. Nikitin, A. V. Arsenin, and V. S. Volkov, Broadband optical properties of monolayer and bulk MoS2, npj 2D Mater. Appl. 4, 21 (2020).
  * Islam  _et al._ [2021] K. M. Islam, R. Synowicki, T. Ismael, I. Oguntoye, N. Grinalds, and M. D. Escarra, In-plane and out-of-plane optical properties of monolayer, few-layer, and thin-film MoS2 from 190 to 1700 nm and their application in photonic device design, Adv. Photonics Res. 2, 2000180 (2021).
  * Ermolaev  _et al._ [2021] G. A. Ermolaev, D. V. Grudinin, Y. V. Stebunov, K. V. Voronin, V. G. Kravets, J. Duan, A. B. Mazitov, G. I. Tselikov, A. Bylinkin, D. I. Yakubovsky, S. M. Novikov, D. G. Baranov, A. Y. Nikitin, I. A. Kruglov, T. Shegai, P. Alonso-González, A. N. Grigorenko, A. V. Arsenin, K. S. Novoselov, and V. S. Volkov, Giant optical anisotropy in transition metal dichalcogenides for next-generation photonics, Nat. Commun. 12, 854 (2021).
  * Born and Wolf [1959] M. Born and E. Wolf,  _Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light_ (Pergamon Press Ltd., 1959) Chap. 14.
  * Yariv and Yeh [1984] A. Yariv and P. Yeh,  _Optical Waves in Crystals: Propagation and Control of Laser Radiation_ (Wiley, 1984) Chap. 6.
  * Jayasekara  _et al._ [2016] C. Jayasekara, M. Premaratne, S. D. Gunapala, and M. I. Stockman, MoS2 spaser, J. Appl. Phys. 119, 133101 (2016).
  * Ren  _et al._ [2021] Q. Ren, F. Feng, X. Yao, Q. Xu, M. Xin, Z. Lan, J. You, X. Xiao, and W. E. I. Sha, Multiplexing-oriented plasmon-MoS2 hybrid metasurfaces driven by nonlinear quasi bound states in the continuum, Opt. Express 29, 5384 (2021).
  * Wang  _et al._ [2022] X. Wang, J. Xin, Q. Ren, H. Cai, J. Han, C. Tian, P. Zhang, L. Jiang, Z. Lan, J. You, and W. E. I. Sha, Plasmon hybridization induced by quasi bound state in the continuum of graphene metasurfaces oriented for high-accuracy polarization-insensitive two-dimensional sensors, Chin. Opt. Lett. 20, 042201 (2022).
  * Srivastava  _et al._ [2017] Y. K. Srivastava, A. Chaturvedi, M. Manjappa, A. Kumar, G. Dayal, C. Kloc, and R. Singh, MoS2 for ultrafast all-optical switching and modulation of THz Fano metaphotonic devices, Adv. Opt. Mater. 5, 1700762 (2017).
  * Bucher  _et al._ [2019] T. Bucher, A. Vaskin, R. Mupparapu, F. J. F. Löchner, A. George, K. E. Chong, S. Fasold, C. Neumann, D.-Y. Choi, F. Eilenberger, F. Setzpfandt, Y. S. Kivshar, T. Pertsch, A. Turchanin, and I. Staude, Tailoring photoluminescence from MoS2 monolayers by Mie-resonant metasurfaces, ACS Photonics 6, 1002 (2019).
  * Muhammad  _et al._ [2021] N. Muhammad, Y. Chen, C.-W. Qiu, and G. P. Wang, Optical bound states in continuum in MoS2-based metasurface for directional light emission, Nano Lett. 21, 967 (2021).
  * Li  _et al._ [2022] H. Li, G. Wei, H. Zhou, H. Xiao, M. Qin, S. Xia, and F. Wu, Polarization-independent near-infrared superabsorption in transition metal dichalcogenide Huygens metasurfaces by degenerate critical coupling, Phys. Rev. B 105, 165305 (2022).
  * Zhang  _et al._ [2022] Y. Zhang, D. Chen, W. Ma, S. You, J. Zhang, M. Fan, and C. Zhou, Active optical modulation of quasi-BICs in Si–VO2 hybrid metasurfaces, Opt. Lett. 47, 5517 (2022).
  * Evlyukhin  _et al._ [2010] A. B. Evlyukhin, C. Reinhardt, A. Seidel, B. S. Luk’yanchuk, and B. N. Chichkov, Optical response features of Si-nanoparticle arrays, Phys. Rev. B 82, 045404 (2010).
  * Kuznetsov  _et al._ [2016] A. I. Kuznetsov, A. E. Miroshnichenko, M. L. Brongersma, Y. S. Kivshar, and B. Luk’yanchuk, Optically resonant dielectric nanostructures, Science 354, aag2472 (2016).
  * Evlyukhin  _et al._ [2011] A. B. Evlyukhin, C. Reinhardt, and B. N. Chichkov, Multipole light scattering by nonspherical nanoparticles in the discrete dipole approximation, Phys. Rev. B 84, 235429 (2011).
  * Khromova  _et al._ [2016] I. Khromova, P. Kužel, I. Brener, J. L. Reno, U.-C. Chung Seu, C. Elissalde, M. Maglione, P. Mounaix, and O. Mitrofanov, Splitting of magnetic dipole modes in anisotropic TiO2 micro-spheres, Laser Photon. Rev. 10, 681 (2016).
  * Razmjooei  _et al._ [2021] N. Razmjooei, Y. H. Ko, F. A. Simlan, and R. Magnusson, Resonant reflection by microsphere arrays with AR-quenched Mie scattering, Opt. Express 29, 19183 (2021).
  * Decker  _et al._ [2015] M. Decker, I. Staude, M. Falkner, J. Dominguez, D. N. Neshev, I. Brener, T. Pertsch, and Y. S. Kivshar, High-efficiency dielectric Huygens’ surfaces, Adv. Opt. Mater. 3, 813 (2015).
  * Fedotov  _et al._ [2007] V. A. Fedotov, M. Rose, S. L. Prosvirnin, N. Papasimakis, and N. I. Zheludev, Sharp trapped-mode resonances in planar metamaterials with a broken structural symmetry, Phys. Rev. Lett. 99, 147401 (2007).
  * Koshelev  _et al._ [2018] K. Koshelev, S. Lepeshov, M. Liu, A. Bogdanov, and Y. Kivshar, Asymmetric metasurfaces with high-$Q$ resonances governed by bound states in the continuum, Phys. Rev. Lett. 121, 193903 (2018).
  * Khardikov  _et al._ [2010] V. V. Khardikov, E. O. Iarko, and S. L. Prosvirnin, Trapping of light by metal arrays, J. Opt. 12, 045102 (2010).
  * Kupriianov  _et al._ [2019] A. S. Kupriianov, Y. Xu, A. Sayanskiy, V. Dmitriev, Y. S. Kivshar, and V. R. Tuz, Metasurface engineering through bound states in the continuum, Phys. Rev. Applied 12, 014024 (2019).
  * Xiao  _et al._ [2022] S. Xiao, M. Qin, J. Duan, F. Wu, and T. Liu, Polarization-controlled dynamically switchable high-harmonic generation from all-dielectric metasurfaces governed by dual bound states in the continuum, Phys. Rev. B 105, 195440 (2022).
  * He  _et al._ [2018] Y. He, G. Guo, T. Feng, Y. Xu, and A. E. Miroshnichenko, Toroidal dipole bound states in the continuum, Phys. Rev. B 98, 161112 (2018).
  * van Hoof  _et al._ [2021] N. J. van Hoof, D. R. Abujetas, S. E. Ter Huurne, F. Verdelli, G. C. Timmermans, J. A. Sánchez-Gil, and J. G. Rivas, Unveiling the symmetry protection of bound states in the continuum with terahertz near-field imaging, ACS Photonics 8, 3010 (2021).
  * Evlyukhin  _et al._ [2020] A. B. Evlyukhin, V. R. Tuz, V. S. Volkov, and B. N. Chichkov, Bianisotropy for light trapping in all-dielectric metasurfaces, Phys. Rev. B 101, 205415 (2020).
  * Tsilipakos  _et al._ [2021] O. Tsilipakos, L. Maiolo, F. Maita, R. Beccherelli, M. Kafesaki, E. E. Kriezis, T. V. Yioultsis, and D. C. Zografopoulos, Experimental demonstration of ultrathin broken-symmetry metasurfaces with controllably sharp resonant response, Appl. Phys. Lett. 119, 231601 (2021).
  * Tuz  _et al._ [2018] V. R. Tuz, V. V. Khardikov, A. S. Kupriianov, K. L. Domina, S. Xu, H. Wang, and H.-B. Sun, High-quality trapped modes in all-dielectric metamaterials, Opt. Express 26, 2905 (2018).
  * Yu  _et al._ [2019] P. Yu, A. S. Kupriianov, V. Dmitriev, and V. R. Tuz, All-dielectric metasurfaces with trapped modes: Group-theoretical description, J. Appl. Phys. 125, 143101 (2019).
  * Bi  _et al._ [2021] K. Bi, Q. Wang, J. Xu, L. Chen, C. Lan, and M. Lei, All-dielectric metamaterial fabrication techniques, Adv. Opt. Mater. 9, 2001474 (2021).
  * Evlyukhin  _et al._ [2021] A. B. Evlyukhin, M. A. Poleva, A. V. Prokhorov, K. V. Baryshnikova, A. E. Miroshnichenko, and B. N. Chichkov, Polarization switching between electric and magnetic quasi-trapped modes in bianisotropic all-dielectric metasurfaces, Laser Photonics Rev. 15, 2100206 (2021).
  * Prokhorov  _et al._ [2022a] A. V. Prokhorov, A. V. Shesterikov, M. Y. Gubin, V. S. Volkov, and A. B. Evlyukhin, Quasitrapped modes in metasurfaces of anisotropic ${\mathrm{MoS}}_{2}$ nanoparticles for absorption and polarization control in the telecom wavelength range, Phys. Rev. B 106, 035412 (2022a).
  * Prokhorov  _et al._ [2022b] A. V. Prokhorov, P. D. Terekhov, M. Y. Gubin, A. V. Shesterikov, X. Ni, V. R. Tuz, and A. B. Evlyukhin, Resonant light trapping via lattice-induced multipole coupling in symmetrical metasurfaces, ACS Photonics 9, 3869 (2022b).
  * Qin  _et al._ [2022] M. Qin, J. Duan, S. Xiao, W. Liu, T. Yu, T. Wang, and Q. Liao, Strong coupling between excitons and quasi-bound states in the continuum in the bulk transition metal dichalcogenides, arXiv preprint arXiv:2209.00416 (2022).
  * Ushkov  _et al._ [2022] A. A. Ushkov, G. A. Ermolaev, A. A. Vyshnevyy, D. G. Baranov, A. V. Arsenin, and V. S. Volkov, Anapole states and scattering deflection effects in anisotropic van der Waals nanoparticles, Phys. Rev. B 106, 195302 (2022).
  * Snitzer [1961] E. Snitzer, Cylindrical dielectric waveguide modes, J. Opt. Soc. Am. 51, 491 (1961).
  * Mongia and Bhartia [1994] R. K. Mongia and P. Bhartia, Dielectric resonator antennas—a review and general design relations for resonant frequency and bandwidth, Int. J. RF Microw. Comput. Aided Eng. 4, 230 (1994).
  * Sayanskiy  _et al._ [2019] A. Sayanskiy, A. S. Kupriianov, S. Xu, P. Kapitanova, V. Dmitriev, V. V. Khardikov, and V. R. Tuz, Controlling high-$Q$ trapped modes in polarization-insensitive all-dielectric metasurfaces, Phys. Rev. B 99, 085306 (2019).
  * Vaity  _et al._ [2022a] P. Vaity, H. Gupta, A. Kala, S. Dutta Gupta, Y. S. Kivshar, V. R. Tuz, and V. G. Achanta, Polarization-independent quasibound states in the continuum, Adv. Photonics Res. 3, 2100144 (2022a).
  * Tuz  _et al._ [2020] V. R. Tuz, P. Yu, V. Dmitriev, and Y. S. Kivshar, Magnetic dipole ordering in resonant dielectric metasurfaces, Phys. Rev. Applied 13, 044003 (2020).
  * Tuz and Evlyukhin [2021] V. R. Tuz and A. B. Evlyukhin, Polarization-independent anapole response of a trimer-based dielectric metasurface, Nanophotonics 10, 4373 (2021).
  * Vaity  _et al._ [2022b] P. Vaity, H. Gupta, A. Kala, S. Dutta Gupta, Y. S. Kivshar, V. R. Tuz, and V. Gopal Achanta, Polarization-independent quasi-bound states in the continuum, Adv. Photonics Res. 3, 2100144 (2022b).
  * Tonning [1982] A. Tonning, Circularly symmetric optical waveguide with strong anisotropy, IEEE Trans. Microwave Theory Techn. 30, 790 (1982).
  * Paul and Shevgaonkar [1981] D. K. Paul and R. K. Shevgaonkar, Multimode propagation in anisotropic optical waveguides, Radio Sci. 16, 525 (1981).
  * Yu  _et al._ [2018] P. Yu, V. I. Fesenko, and V. R. Tuz, Dispersion features of complex waves in a graphene-coated semiconductor nanowire, Nanophotonics 7, 925 (2018).
  * Vandenbulcke and Lagasse [1976] P. Vandenbulcke and P. E. Lagasse, Eigenmode analysis of anisotropic optical fibers or integrated optical waveguides, Electron. Lett. 12, 120 (1976).
  * Koshiba  _et al._ [1986] M. Koshiba, K. Hayata, and M. Suzuki, Finite-element solution of anisotropic waveguides with arbitrary tensor permittivity, J. Lightwave Technol. 4, 121 (1986).
  * Bagatskaya  _et al._ [2005] O. V. Bagatskaya, V. I. Fesenko, and S. N. Shulga, Electric field lines in a rectangular waveguide with an inhomogeneous anisotropic insert, Telecommun. Radio Eng. 63, 1 (2005).
  * Dai and Jen [1991] J. D. Dai and C. K. Jen, Analysis of cladded uniaxial single-crystal fibers, J. Opt. Soc. Am. A 8, 2021 (1991).
  * Dai  _et al._ [1992] J. D. Dai, C. A. S. de Oliveira, and C. K. Jen, Dispersion characteristics of cladded biaxial optical fibers, J. Opt. Soc. Am. A 9, 1564 (1992).
  * Evlyukhin and Chichkov [2019] A. B. Evlyukhin and B. N. Chichkov, Multipole decompositions for directional light scattering, Phys. Rev. B 100, 125415 (2019).