Related to the question : https://mathoverflow.net/questions/242136/why-m-1-subset-m-2-not-rightarrow-n-m-1-lambda-leq-n-m-2-lambda
The Neumann eigenvalues of the rectangle with sides $a$ and $b$ are $$\nu_{k,l}=\frac{(\pi k)^2}{a^2}+\frac{(\pi l)^2}{b^2},$$ with $k,l \in \mathbb{N}_0$. So assuming that $a>b$, the first $3$ eigenvalues are $\nu_1=0$, $\nu_2=\frac{\pi^2}{a^2}$, and $\nu_3 = \frac{\pi^2}{b^2}.$ We pick $1 < a < \sqrt{2}$, and choose $b>0$ sufficiently small, so that the rectanglecan be place inside the unit squre. For the unit squre, the first $3$ Neumann eigenvalues are $\nu_1 ' = 0$, $\nu_2 ' = \pi^2$, and $\nu_3 ' = \pi^2$. Since $a>1$, we have $\nu_2 < \nu_2 '$, which could not happen if domain monotonicity were true.
Does this example work? If so, since the spectrum of the rectangle is the same as Dirichlet condition, why it is a counterexample for NBC but not for DBC (explanation in details please)
As mickep observes in the comments, the spectra are not the same! The Dirichlet spectrum is $$ \lambda_{kl} = \bigg(\frac{\pi k}{a}\bigg)^2 + \bigg(\frac{\pi l}{b}\bigg)^2,\ k,l = 1,2,3,\ldots $$ while the Neumann spectrum is $$ \nu_{kl} = \bigg(\frac{\pi k}{a}\bigg)^2 + \bigg(\frac{\pi l}{b}\bigg)^2,\ k,l = 0,1,2,3,\ldots $$ The intuition is that the Neumann spectrum explicitly "hears" the individual side lengths $a,b$, while in the Dirichlet spectrum the two are muddled together. This is the idea behind the counterexample given in Christian Remling's answer to your MathOverflow question, where given a rectangle you can create a smaller rectangle with a longer side that embeds.