Failed to calculate a BIC mode in 1D grating

[Vimim3]

Recently, I read a PRB paper with the title "Bound states in the continuum in symmetric and asymmetric photonic crystal slabs". I try to repeat the result in Fig. 2(b) and (c). I found the results given by legume do not include the second BIC mode. The paper and my program are attached. Is there any problem in my program?
PhysRevB.101.155303.pdf
BIC in 1D grating.zip

[momchil-flex]

I think you may have hit a case which is strongly affected by the main source of inaccuracy in GME: the incompleteness of the basis. Specifically, we only use fully guided modes in the expansion, which means we do not use slab modes which are above the light line. Let's number the bands in the paper 1 through 4 with increasing frequency at Gamma (so first two black bands, then green, then blue). I think we can sort of identify bands 1, 3, and 4 in the legume figure. However, band 2 is completely missing: I suspect that this one may be coming from the leaky modes of a higher-order slab band (e.g. gmode_inds = 2). These leaky modes are note included in the expansion however; only the fully guided ones past the light cone are. This may also be why nearby bands are not right, and specifically the BIC of band 4 (the blue band) appears at a higher k in legume. I checked the field profiles and that band does look like what's plotted on the right in panel (d).

We've had some in-depth discussions on remedying this GME deficiency, but the problem is that the basis then becomes quasi-normal, and it's... not easy to define things then. It's a bit of an open research project how to make this work.

In short the good news is that I don't think anything is wrong in your script, but the bad news is that I don't think the results can be further improved...

[Vimim3]

I see. I hope this problem can be solved in the future. Thank you very much!