You get the boundedness from the Hardy space $\mathscr H^1$ into $L^1$ with the atomic decomposition and by duality follows the boundedness from $L^\infty$ to BMO (which is the dual of $\mathscr H^1$).

In fact I refer to the fact that the maximal function sends $L^1$ in $L^1_w$ and then I use the $L^2$ bounedness and the Marcinkiewicz Interpolation Theorem to get the $L^p$ boundedness for $p\in (1,2]$ and duality for finite exponents.

Take $a=1/2$, $f=0$, so that the equation reads $2x z'=z$. Multiplying both sides by $z$, calling $y=z^2$, you get $xy'=y$ and thus $y=\alpha x$. If $z$ is analytic, $\alpha\not=0$ and $\zeta= \alpha^{-1/2}z$, you found an analytic function $\zeta$ such that $\zeta^2=x$, which is not possible. Thus $z=0$.

Take $a=1/2$, $f=0$, so that the equation reads $2x z'=z$. Multiplying both sides by $z$, calling $y=z^2$, you get $xy'=y$ and thus $y=\alpha x$. If $z$ is analytic, $:alpha\not=0$,and $\zeta= \alpha^{-1/2}z$, you found an analytic function $\zeta$ such that $\zeta^2=x$, which is not possible. Thus $z=0$.

@Jack Yes, but it is not completely obvious, you have to write$$u(x,y)=u(a,y)+\int_a^x(\partial_1u)(s,y) ds$$ and the similar formula for $u(x,y)=u(x,b)+\dots$. Then you have an expression for $\vert u(x,y)\vert^2$ that you can apply to $u=v^2$.

No, you just have to think about the continuous canonical injection mentioned above.

I have added a comment on periodic functions in my answer.

Thanks a lot, I am indeed very grateful for your help. I do not believe that Weyl's law is enough for the stability of gaps since it probably could be the same with eigenvalues spread on the interval $(k,k+1)$ instead of being clustered at the integers. Also, it is surprising to me that the references that you kindly indicate are so recent; I should say that I thought that, for a smooth perturbation it was either wrong or classical.