I ran the first configuration for 1550 nm earlier, which is why I was getting different number of modes. I reran the simulations and my results are very close to yours. I don’t think the third mode is due to numerical errors. I believe it is a higher order TE mode.
Also, one thing that confused me in the paper is that at some places they have mentioned 1310 nm while at others 1330 nm. I did find a more recent textbook, that had results for these configurations and they are very close to yours. But the caveat is that they have also used Lumerical for simulating these structures.
And sorry I wasn't able to look at your files as we are not allowed to exchange files here.