Nonlocal Hormander's hypoellipticity theorem

Consider the following nonlocal integro-differential operator: for $\alpha\in(0,2)$, $$\cal L^{(\alpha)}_{\sigma,b} f(x):=\mbox{p.v.} \int_{\mathbb{R}^d-\{0\}}\frac{f(x+\sigma(x)z)-f(x)}{|z|^{d+\alpha}}d z+b(x)\cdot\nabla f(x),$$ where $\sigma:\mathbb{R}^d\to\mathbb{R}^d\times\mathbb{R}^d$ and $b:\mathbb{R}^d\to\mathbb{R}^d$ are two $C^\infty_b$-functions, and p.v. stands for the Cauchy principal value. Let $B_1(x):=\sigma(x)$ and $B_{j+1}(x):=b(x)\cdot\nabla B_j(x)-\nabla b(x)\cdot B_j(x)$ for $j\in\mathbb{N}$. Under the following H\"ormander's type condition: for any $x\in\mathbb{R}^d$ and some $n=n(x)\in\mathbb{N}$, $$\mathrm{Rank}[B_1(x), B_2(x),\cdots, B_n(x)]=d,$$ by using the Malliavin calculus, we prove the existence of the heat kernel $\rho_t(x,y)$ to the operator $\cal L^{(\alpha)}_{\sigma,b}$ as well as the continuity of $x\mapsto \rho_t(x,\cdot)$ in $L^1(\mathbb{R}^d)$ for each $t>0$. Moreover, when $\sigma(x)=\sigma$ is constant, under the following uniform H\"ormander's type condition: for some $j_0\in\mathbb{N}$, $$\inf_{x\in\mathbb{R}^d}\inf_{|u|=1}\sum_{j=1}^{j_0}|u B_j(x)|^2>0,$$ we also prove the smoothness of $(t,x,y)\mapsto\rho_t(x,y)$ with $\rho_t(\cdot,\cdot)\in C^\infty_b(\mathbb{R}^d\times\mathbb{R}^d)$ for each $t>0$.