I'm trying to figure out the regularity for solutions to the following integral equation: \begin{align} \begin{cases} \displaystyle{\frac{d}{dt}}u(t,x)&=\displaystyle{\frac{1}{\sigma_d\epsilon^{d+2}}}\int_{\mathbb{T}^d}\sin\left(u(t,y)-u(t,x)\right)1_{\{|y-x|\le\epsilon\}}dy, \label{eq:int}\\[10pt] u|_{t=0}&=u_0. \end{cases} \end{align} Here $\epsilon>0$ is a parameter and $\mathbb{T}^d$ is the $d-$dimensional torus. One can construct solutions that live in the space $C([0,\infty),L^1(\mathbb{T}^d))$ via a fixed point argument. My question is: are solutions to this equation as good as the initial condition $u_0$ in the spatial domain (i.e. $u_0\in C^\infty(\mathbb{T}^d)\Rightarrow u\in C^\infty(\mathbb{T}^d)$)?
My first guess is yes, becuase event if the integral term does not regularize solutions (since there is no derivative on the right hand side as, say, in the heat equation) it doesn't make them worst. I thought of something like an Euler's method to prove it: $$ u_n(t,x)=u_{n-1}(t,x)+h\displaystyle{\frac{1}{\sigma_d\epsilon^{d+2}}}\int_{\mathbb{T}^d}\sin\left(u_{n-1}(t,y)-u_{n-1}(t,x)\right)1_{\{|y-x|\le\epsilon\}}dy. $$ Since $u_0$ is smooth in $x$, so is $u_1$ and so on.