Tailoring accidental double bound states in the continuum in all-dielectric metasurfaces

Bound states in the continuum (BICs) have been thoroughly investigated due to their formally divergent Q-factor, especially those emerging in all-dielectric, nanostructured metasurfaces from symmetry protection at the $\Gamma$ point (in-plane wavevector $k_{||}=0$). Less attention has been paid to accidental BICs that may appear at any other $k_{||}\not =0$ in the band structure of supported modes, being in turn difficult to predict. Here we make use of a coupled electric/magnetic dipole model to determine analytical conditions for the emergence of accidental BICs, valid for any planar array of meta-atoms that can be described by dipolar resonances, which is the case of many nanostructures in the optical domain. This is explored for all-dielectric nanospheres through explicit analytical conditions that allow us in turn to predict accidental BIC positions in the parameter space $(\omega,\bf{k_{||}}$). Finally, such conditions are exploited to determine not only single, but also double (for both linear polarizations) accidental BICs occurring at the same position in the dispersion relation $\omega-\bf{k_{||}}$ for realistic semiconductor nanodisk meta-atoms. This might pave the way to a variety of BIC-enhanced light-matter interaction phenomena at the nanoscale such as lasing or non-linear conversion, that benefit from emerging at wavevectors away from the $\Gamma$ point (off-normal incidence) overlapping for both linear polarizations.


Introduction
Planar arrays of particles, often called metasurfaces, are versatile thin platforms designed to control the properties of light at BICs have been largely explored and they can only be found at the Γ point, where there is a mismatch between the symmetry of the field mode and the available radiation channel imposed by the in-phase oscillation of the particles in the array. By contrast, accidental BICs emerge when different modes are coupled and they interfere destructively at a specific angle, typically arising in high-order resonant bands in photonic crystal slabs. [19,34,36,41,[49][50][51] Despite both photonic crystals and metasurfaces can be considered analogous in the nondiffracting regime as quasi-planar periodic arrangements, there are subtle theoretical differences in the manner that connected domains (slab with holes in photonic crystals) polarize light, as compared with isolated domains (meta-atoms in the case of metasurfaces) that allow for resonant electromagnetic confinement inside, with nontrivial consequences in the resulting phenomenology; and Babinet-like principles are not at all evident as long as holes/meta-atoms are thick and penetrable. [52] In the case of metasurfaces, since particles in the array suffer from depolarization effects from the own lattice, the conditions for the formation of accidental BICs are not evident, and there are no guidelines to design them in the wavevector space; in fact, it has been assumed that multipolar contributions are needed to yield accidental BICs in metasurfaces, [49] condition that we will demonstrate not strictly necessary, as in the case of 1D gratings, [23] with periodicity only along one in-plane direction, invariant along the other.
In this work, the emergence of accidental BICs in a rectangular metasurface of dielectric particles is investigated on the basis of a coupled electric and magnetic dipole (CEMD) formulation. [53,54] First, generic conditions are established, demonstrating that, for arrays of axially symmetric particles, accidental BIC can only arise along the symmetry lines ΓX (and ΓY for rectangular arrays). Analytic expressions are then derived to determine the condition for accidental BICs along these symmetry lines, showing that the individual polarizability of the particles must be related to the lattice depolarization Green function. These conditions only depend on the geometry of the lattice and meta-atom polarizabilities, enabling us to establish a general guideline for the formation of accidental BICs. Later, the accidental BIC condition for rectangular arrays of alldielectric spheres is studied, giving examples of systems supporting TE and TM accidental BICs. Interestingly, it is shown that rectangular arrays can support more than one accidental BIC for the same metasurface at a given polarization. Finally, a feasible example with a square array of semiconductor nanodisks is presented, where the aspect ratio of the disk offers us a new degree of freedom for engineering the polarizability of the meta-atom. In addition, rectangular metasurfaces also allow to tailor the emergence of, not only single, but also double (for both linear polarizations) accidental BICs at the same position in ω, k-space (ω being the frequency and k the in-plane wavevector), at any off-normal angle of incidence along ΓX or ΓY for a given frequency. This is in turn confirmed through full numerical simulations, in good agreement with our quasianalytical CEMD calculations. Although all the phenomenology is analyzed for rectangular arrays and one meta-atom per unit cell, the principles used in this study can be extended to metasurfaces with any geometry as long as the CEMD approach remains valid.

CEMD Formulation
First, let us consider an infinite rectangular array of identical particles, labeled as (n, m) and placed atˆˆ= where a and b are the lattice constants along the x and y axis, respectively. For the homogeneous problem in the absence of external illumination, each particle in the array is excited by the waves emitted from the rest of the array. The self-consistent incident field on the (n, m) = (0, 0) particle ( inc Ψ , is then given by the solution of where ∑ nm ′ means that the sum runs for all indices except for (n, m) = (0, 0). ( ) G r and α are matrices representing the dyadic Green function and the dipolar polarizability of the particles, respectively, and their representations depend on the basis chosen to describe the electromagnetic fields, ( ) r Ψ Ψ . For the fields, the time dependence exp (−iωt) is assumed; ω is the angular frequency, related to the modulus of the wavevector through k = ω/c, c being the speed of light. The dyadic Green function is obtained from the scalar Green function, ( ) g r , by applying a linear differential operator, L, that also depends on the chosen basis. [53] For periodic arrays the Bloch's theorem holds, where k x and k y are the in-plane components of the wavevectors. Thus, the self-consistent incident field can be written as We have defined G b , the lattice depolarization dyadic (or return Green function), as and we have explicitly shown the dependence on both the modulus wavevector and the in-plane components of the wavevector, k k ( , ) || = k k x y . G b tells us about the coupling strength between particles, and is crucial to determine all the lattice properties and the nature of the modes supported by the metasurface. Finally, we rearrange Equation (3) as follows: where I is the unit dyadic. In order to find the resonant states of the metasurfaces we need to find a solution to the homogeneous linear system of equations Equation (5), appearing only when the determinant is equal to zero: The complex frequencies at which Equation (6) is satisfied are the eigenfrequencies, denoted by ν = ν′ − iν″, where ν′ and ν″ are real numbers that denote the real and imaginary parts, respectively. In the latter eigenmode equation, it is more convenient to use the second expression rather than the first one because the imaginary part is well defined. Since G b can be expressed in the reciprocal space, [53] Equation (6) can thus be employed to determine the dispersion relation of resonant modes in metasurfaces. For the sake of simplicity, we will assume in what follows that the array is embedded in a uniform medium (vacuum, indeed), so that G in Equation (3) is the free space Green dyadic. This is equivalent in an experimental configuration with a metasurface on a substrate to adding an index-matching layer. [27] Nonetheless, the impact of a substrate could be straightforwardly incorporated in our CEMD formulation by including the Green dyadic for two semi-infinite media.

Generic Condition for Accidental BICs
In the more generic case (with any restriction for the values of k||), the determinant can be only factorized into two terms (for diagonal polarizabilities), one for each polarization: x yz xyz with 1 1 1 In these expressions, G bij are the matrix elements of G b , ( ),( ) α i e m are the diagonal matrix elements of the electric (e) and magnetic (m) polarizability, and TM, TE stand for transverse electric and magnetic modes, respectively. To simplify, the upper index (EM) used in ref. [ 53 ] (indicating that they connect electric and magnetic dipoles), is omitted in G byz and G bzx , understood for physical reasons. Note that G bxy couples either electric only or magnetic only dipoles, while G byz and G bzx connect electric with magnetic dipoles. Also, for these expressions, the sign of G byz and G bzx is taken in such a way that: since their signs must be consistent with their definitions. The resonant surface modes arising from the zeroes of (TE) η xyz (respectively, (TM) η xyz ) represent hybrid modes with different polarization, where the electric (respectively, magnetic) in-plane dipoles are coupled among them (by G bxy ) and with the magnetic (respectively, electric) dipoles along the z axis (by G byz and G bzx ). The interference between the different dipolar modes leads to interesting phenomenology, with the open question yet on whether accidental BICs can be supported: recall that the imaginary part of the above mentioned resonant modes must vanish (ν″ = 0) to be considered BICs.
To this end, let us analyze Equations (8) in a general manner. First, note that it is possible to solve for the zeroes of each of Equations (8) as a system of two equations for the real and imaginary parts (separately) with three unknowns: the real parts of the inverse of the polarizabilities along the three axis. Then the accidental BIC conditions can be solved for either mode as a function of a parameter (the real part of the inverse of one of the polarizabilities), leading to plausible solutions (zeroes) of Equations (8) along any direction in k-space by tuning the remaining polarizabilities. This is true for non axially symmetric meta-atoms (as an ellipsoidal disk or a regular disk with its axis in the xy plane) with three different polarizabilities along the main axes.
Nonetheless, if the particles possess axial symmetry around the z axis (as in a disk), in-plane polarizabilities are equal, α x = α y , and less degrees of freedom are available. In fact, it can be shown that, if both 0 ≠ k x and 0 ≠ k y , no real values for the inverse of the polarizability can be found for which Equation (7) is fulfilled. Therefore, it is not possible to design accidental BIC away from the symmetry lines ΓX and ΓY (for rectangular arrays) for axially symmetric arrays of particles (in the dipolar regime). Moreover, no accidental BIC can be formed either along the symmetry lines ΓM and ΓS.
By contrast, G bxy and G bzx (or G byz ) become zero along the symmetry line ΓX (respectively, ΓY), and only the interference between the electric dipole along the y axis (respectively, x axis) and the magnetic dipole along the z axis is relevant (and vice versa exchanging electric/magnetic). Therefore, since the degrees of freedom remain the same, polarizabilities satisfying Equation (7) can be found, so that accidental BICs along the symmetry line ΓX (or ΓY) can be supported by metasurfaces with axially symmetric meta-atoms, as we reveal in the following sections.

Accidental BICs along ΓX
If we thus consider modes propagating along the x axis ( 0 ≠ k x and k y = 0) and diagonal polarizabilities, Equation (6) can be factorized in four terms [53] with 1 , For waves propagating along the y axis (k x = 0 and 0 ≠ k y ), the phenomenology is the same by replacing the index y by x in Equation (11). Among the different resonant surface modes represented by each term in Equation (10), we exclude (TE) η x and (TM) η x since they have been shown to yield no BIC whatsoever (all their poles inside the continuum of radiation have non-negligible widths, 0 ν ′′ ≠ ). We thus focus on the hybrid modes, (TE) η yz and (TM) η yz , emerging from electric (magnetic) dipoles along the y axis coupled to magnetic (electric) dipoles along the z axis. At the Γ point, k x = k y = 0, well known symmetry-protected BICs arise. [53] For the sake of completeness, we explicitly show the conditions, not shown in ref. [ 53 ], yielding symmetry-protected BICs given by the in-phase oscillation of lossless dipoles along the z axis: Therefore, for a metasurface to support symmetry-protected BICs, the real part of the inverse of the polarizability of the meta-atoms must be consistent with the lattice properties. Searching for accidental BICs at 0 ≠ k x (with k y = 0), also for lossless meta-atoms, the condition η yz = 0 yields a system of two equations (real and imaginary parts) with two unknowns: the real parts of the inverse of the polarizabilities along the y and z axes. After a straightforward derivation, the accidental BIC conditions are obtained note that two conditions, for 1/α y and 1/α z separately, must be satisfied simultaneously.
Equations (12) and (13) give us an interesting perspective about the differences between symmetry-protected and accidental BICs. Symmetry-protected BICs only need to fulfil one condition at a specific frequency for the polarizability (with in-plane modulus wavevector k || = 0), whereas accidental BICs need to fulfil two conditions at the same frequency and longitudinal component of the wavevector. Then, accidental BICs obey more restrictive conditions that are not always reachable, depending on the precise cancellation of radiation from different dipolar modes. These characteristics are illustrated in Figure 1: a symmetry-protected BIC may simply arise from vertical dipolar resonances, whereas an accidental BIC requires the interference between different in-plane and out-of-plane dipolar modes to cancel radiation to the far-field at specific offnormal angles. Quite importantly, Equation (13) is a guideline to engineer accidental BICs for a wide variety of all-dielectric (and plasmonic) photonic metasurfaces with a common (rectangular) lattice symmetry.

Accidental BICs on All-Dielectric Nanosphere Metasurfaces
To shed light on the phenomenology resulting from the previous theoretical formulation, we represent in Figure 2 the accidental BIC position as a function both of the normalized frequency (a/λ, with λ = 2πc/ω) and of the angle of incidence (θ) for rectangular arrays (lattice parameters a, b) of dielectric spheres with different values of their refractive index. The polarizability of the spheres lated from Mie theory, taking into account only the first two (dipolar) terms of the harmonic expansion. The color of each line is correlated with the normalized sphere radius r/a that the metasurface must fulfil in order to support an accidental BIC, where r/a = 0.5 represents a metasurface of touching spheres. The results for TE accidental BICs in square arrays (a, b) are analyzed in Figure 2a. It can be seen that the general tendency is that a higher r/a ratio is needed to support accidental BICs as the refractive index of the particles becomes smaller and as the longitudinal component of the wavevector increases. Also, note that the curves corresponding to n s = 2.8 and n s = 3 end before achieving the condition k = k x (being k x related to θ through the expression k x = ksin θ). Indeed, it can be found that no accidental BIC can be formed for approximately n s < 2.5. On the other hand, as the value of n s increases, accidental BICs become only accessible at  higher values of k x , reaching the diffraction limit, as illustrated in Figure 2a for n s = 4.
By contrast, we have found no accidental BICs in TM polarization for square arrays of all-dielectric spheres. For dielectric spheres, the dipolar mode ( ) α y m resonates at frequencies lower than that of ( ) α z e , precluding its formation according to Equations (13). Alternatively, we study TM accidental BICs in rectangular arrays that indeed support them, as revealed in Figure 2d for a/b = 0.66. As in the case of TE polarization, larger spheres are needed for lower refractive indices. Interestingly, these rectangular metasurfaces support accidental BICs for dielectric spheres with relative small values of n s ; namely, high-refractive indices are not strictly necessary, which widens the range of dielectric materials that could be exploited in this respect. For example, accidental BICs can be engineered around θ = 20°-50° for spheres with n s = 1.8 and meaningful radii. For higher values of n s , accidental BICs arise at higher frequencies, being possible to engineer them for n s ≈ 5, but they are not accessible for all values of θ. Also, it is remarkable that two accidental BICs can be found in parameter space for the same n s and r/a parameters (i.e., for the same metasurface). For example, there are accidental BICs at (a/λ, θ) ≈ (0.68, 19) and (a/λ, θ) ≈ (0.59, 42) for n s = 3.5 and r/a = 0.2.
By way of example, we show in Figure 2 the dispersion relations and associated quality factors for modes approaching the accidental BIC condition in four specific metasurfaces. These magnitudes are calculated by solving Equation (6) and using the formalism described in ref. [ 53 ]. The dispersion relations ν′(k x ) are shown in the parameter space given by angle θ (recall that k x = ksin θ) and normalized frequency a/λ (with λ = c/ν′); the quality factor is defined as Q f = ν′/ν″ and calculated along the mode dispersion. TE modes in square arrays are shown in Figure 2b,c for two different metasurfaces: high-refractive index particles n s = 4 and r/a = 0.2 (blue curves); and moderate high-refractive index particles n s = 2.8 and r/a = 0.3 (red curves). As expected, the quality factor of the modes diverges at the predicted position in parameter space of the accidental BICs, around θ = 43° and θ = 18° for the two cases shown in Figure 2b,c. The dispersion relations of TM modes in rectangular arrays (a/b = 0.66) are also represented in Figure 2e,f. The blue curves correspond to spheres with low index of refraction, n s = 1.8 and r/a = 0.37, whereas red curves are for high refraction index particles, n s = 3.5 and r/a = 0.2. Interestingly, both metasurfaces support two accidental BICs for the same polarization, and their relative distance in parameter space can be tailored by the r/a ratio. Incidentally, if r/a is decreased, a threshold value is found at which both accidental BICs collapse into a single one at the same angle and frequency, below which no accidental BIC is found (although a narrow resonance is still present).
For the sake of illustration, the reflectance is shown in Figure 3 for different configurations supporting accidental BICs, marked as black symbols in Figure 2b,e. The reflectance maps have been calculated using the CEMD formulation. [53] For TE polarization, two square arrays have been considered in Figure 3a  index n s = 4 and radius r/a = 0.2 is examined in Figure 3a. The expected behavior in parameter space (formally equivalent to the ω, k) is observed, associated with a narrow (accessible) quasi-BIC band with high Q-factors turning into an (inaccessible) BIC, with a vanishing reflectance corresponding to the BIC at the predicted angle around θ = 43°. The other example is shown in Figure 3b for a lower (but moderately high) refractive index of n s = 2.8 at r/a = 0.3, for which the accidental BIC appears around θ = 18° over a broad high reflectance background.
For TM polarization, we have considered the two examples as marked as black symbols in Figure 2e. First, the interesting case mentioned above for a low index of refraction (n s = 1.8) and r/a = 0.37: as shown in Figure 3c, two accidental BICs are found in a rectangular array (a/b = 0.66) around θ ≈ 25° and 38° . Note that there are very narrow high-reflectance bands between both BIC conditions, barely visible. Indeed, the accidental BIC occurring at large angles is zoomed in in the inset of Figure 3c.
Finally, Figure 3d presents the reflectance for the rectangular array of high refractive index spheres with n s = 3.5 and r/a = 0.2. This metasurface also supports two accidental BICs with an interesting phenomenology. The accidental BIC around θ ≈ 42° shows the typical vanishing feature, that is, there is no background. Conversely, the BIC around θ ≈ 19° occurs over a broad reflectance background with a complex pattern stemming from the combination of a broad resonance; then the quasi-BIC to BIC transition appears as a narrow dip within the broad band, stemming from the interference that fades away within the high reflectance background, as a BIC-induced transparency band. [42] Although all preceding cases have been discussed in normalized frequencies (which could be easily extrapolated to any spectral regime), we would like to emphasize that all of them can be associated to realistic nanosphere metasurfaces supporting accidental BICs in the optical domain, for example, with lattice parameters of the order of a = 300 nm, and dielectric spheres made of low (oxides) or high (semiconductors) refractive indices. To further stress the ability of Equation (13) to design all-dielectric metasurfaces with accidental BICs, we apply it next to a realistic scenario similar to that already used for symmetry-protected BICs. [39]

Accidental BICs on Semiconductor Nanodisk Metasurfaces
In order to give a broader view of accidental BICs in all-dielectric metasurfaces, the phenomenology is also investigated for a more practical case: a square array of dielectric nanocylinders of diameter, D, and height, H. Moreover, the additional degree of freedom, given by the aspect ratio of the particle, offers a new parameter to engineer the properties of accidental BICs. Unlike for spherical particles, the relative spectral position of the resonances can be tuned with the aspect ratio of cylindrical particles. To this end, the polarizabilities of the disks are numerically calculated through SCUFF, [55,56] rigorously accounting for the different components; [57] a refractive index of n = 3.5 is assumed, as a typical value for dielectric semiconductors at optical and near-IR frequencies. [39] Specifically, the polarizabilities of disks with diam- eters D = 120, 150, 180, and 210 nm, and heights H = 80, 110, 140  The accidental BIC position for an array of semiconductor disks with lattice constant a = b = 300 nm (square array) is presented in Figure 4. Therein the position of the symbols marks the point in the parameter space at which the metasurface supports an accidental BIC for the different disks calculated numerically through Equation (13). In addition, the dashed lines indicate the evolution of the spectral position as a function of the diameter, where the polarizability is calculated by interpolating those for different heights at constant diameter. As a general trend, the angular position of the accidental BIC increases with the diameter in both polarizations, up to a point for larger disks that the conditions are no longer fulfilled. In addition, it should be emphasized that accidental BICs can be supported in square arrays of nanodisks for TM polarization, in contrast to the case of a square array of nanospheres.
The appearance of accidental BICs for TM polarization is possible due to the tunability of the relative spectral position of the dipolar resonances of the individual nanodisks with different aspect ratios. Moreover, exploiting such versatility, it is also possible to find configurations at which the metasurface supports accidental BICs for both polarizations at the same position in the parameter space. The requirements are: i) equal electric and magnetic in-plane polarizabilities, ii) equal electric and magnetic out-of-plane polarizabilities; and iii) fulfillment of Equation (13). We have checked that a cylinder of D = 210 nm and H = 200 nm closely meets (i) and (ii) around the wavelength λ ≈ 910 nm, where the parameters of the lattice can be tuned to fulfil the third requirement. Therefore, based on the properties of this specific disk and fixing the operation wavelength at λ = 910 nm, the double accidental BIC condition is found as a function of the lattice parameter along the y axis, b, and the in-plane wavevector along the x axis, k x (where k x = ksin θ), by changing the lattice parameter along the x axis, a. The results are illustrated in Figure 5a. As the lattice parameter a becomes smaller, the angle at which the BIC is found moves to higher angles of incidence, up to the point below a ≈ 420 nm such that no accidental BIC emerges (not shown in the graph). Also, it can be seen that the lattice parameter b must be always bigger than a. Now let us calculate the dispersion relation of the resonant modes in two specific configurations, proceeding as in the case of sphere arrays shown in Figure 2b,e. Note that, for a proper calculation of the dispersion relation, the polarizability at complex frequencies is needed, trivial for the analytical polarizabilities of spheres, but not for those of disks. Nonetheless, since we are interested in the dispersion relation around the accidental BIC, the imaginary part of the mode frequency, ν″, is vanishingly small. Therefore, the polarizability at complex frequencies can be safely approximated to the polarizability at the real part of the complex frequency. This approximation is very accurate for modes with Q f > 10 3 .
First, Figure 5b ′ : note that δ < 2 × 10 −5 in the range θ = 71°-78°. Also, it is remarkable that the quality factor of both modes diverges around θ = 72°, as predicted. A similar phenomenology is observed in a rectangular array with lattice constant a = 550 nm and b = 985 nm. As illustrated in Figure 5d,e, both modes lie down again close to each other in parameter space, although in this configuration the detuning is higher in the vicinity of the minimum. Nonetheless, around θ = 37°, where the quality factor becomes very large, the detuning is near zero. Significantly, the quality factor of the TM mode exhibits two divergences, associated to two accidental BICs, a phenomenology that was also shown above for dielectric sphere arrays in Figure 2f.
By choosing a = 450 nm, the specular reflectance is represented in Figures 6a,b at the double accidental BIC condition. For both polarizations, the accidental BIC is found at θ ≈ 72° (and at λ ≈ 910 nm) when the lattice constant along the y axis is b = 845 nm. Moreover, Figure 6c,d shows the reflectance at both polarizations for a rectangular array with lattice constants a = 550 nm and b = 985 nm. As predicted by Figure 5a, a double accidental BIC is supported at θ ≈ 37°.
Finally, we verify the accuracy of our theoretical model by carrying out numerical calculations through SCUFF [55] (free software implementation based on the method of moments), for a specific case corresponding to Figure 6c Figure 7 confirming that accidental BICs are observed in both polarizations with similar dispersion relations and Q factors. Nevertheless, there is a slight spectral/angular shift with respect to the predictions of our CEMD model, especially in TM polarization: the TE-BIC emerges at θ ≈ 37° and λ ≈ 902 nm, whereas the TM-BIC appears at θ ≈ 41° and λ ≈ 930 nm. This is somewhat expected since our model only retains the dipolar responses. Still, the agreement is remarkable by all means (other resonant bands and related features are reproduced too), bearing in mind in turn that for high resolution reflectance maps, CEMD calculations take a few hours, about two orders of magnitude less than typical numerical simulations. Our CEMD can thus predict in a fast manner the range of parameters wherein double accidental BICs emerge in real configurations, to be fine tuned later on through full numerical simulations in a limited range of disk/lattice parameters.
Bear in mind that double accidental BICs can be very useful to enhance light-matter interactions within all-dielectric nanodisk metasurfaces: this implies that two resonances extending all over the metasurface with extremely large (in principle, diverging) Q-factors can be simultaneously excited with mixed polarized light, with the advantage that the arrays can be designed to obtain ad-hoc positions in the (ω, k) space (namely, wavelength and angle nearly at will).

Conclusions
To summarize, with the help of our coupled electric/magnetic dipole model for infinite planar arrays, valid for meta-atoms properly described by dipolar resonances, we have determined analytical conditions for the emergence of accidental BICs in all-dielectric metasurfaces. This is demonstrated first for all-dielectric spheres in rectangular arrays through explicit conditions that incorporate their dipolar responses through analytical Mie polarizabilities. We have explored the analytical conditions in normalized dimensions that can be easily extrapolated to any spectral regime as long as the refractive index of the spheres remains constant within the studied spectral domain; in particular, realistic nanosphere metasurfaces supporting accidental BICs in the optical domain can be inferred from them, using either oxides or semiconductors as sphere material. Moreover, such conditions are also exploited to design semiconductor metasurfaces supporting accidental BICs in the optical domain, typically fabricated through standard lithographic means, such as nanodisk arrays, by properly accounting for the polarizabilities of the cylinder-shaped meta-atoms. In this regard, thanks to the additional degree of freedom (aspect ratio) allowed by the cylindrical shape as compared to spheres, we are able to determine conditions for the emergence of, not only single, but also double (both linear polarizations) accidental BICs at the same ω and k || : this in turn implies that two resonances extending all over the metasurface with extremely large (in principle, diverging) Q-factors can be simultaneously excited with mixed polarized light, with the advantage that the arrays can be designed to obtain ad hoc positions in the (ω, k) space (namely, wavelength and angle nearly at will). Such phenomenology opens a new venue for BIC-enhanced light-matter interaction in all-dielectric metasurfaces in the optical domain (and throughout the electromagnetic spectrum), paving the way to a variety of phenomena that may benefit from the emergence of two accidental BICs away from the Γ-point for both polarizations, such as nonlinear processes, lasing, chirality, etc.