$\mathbb Q$-linear dependence of certain Bessel moments

Let $I_0$ and $K_0$ be modified Bessel functions of the zeroth order. We use Vanhove's differential operators for Feynman integrals to derive upper bounds for dimensions of the $\mathbb Q$-vector space spanned by certain sequences of Bessel moments $$\left\{\left.\int_0^\infty [I_0(t)]^a[K_0(t)]^b t^{2k+1}\mathrm{d}\, t\right|k\in\mathbb Z_{\geq0}\right\},$$where $a$ and $b$ are fixed non-negative integers. For $a\in\mathbb Z\cap[1,b)$, our upper bound for the $\mathbb Q$-linear dimension is $\lfloor (a+b-1)/2\rfloor$, which improves the Borwein-Salvy bound $\lfloor (a+b+1)/2\rfloor$. Our new upper bound $\lfloor (a+b-1)/2\rfloor$ is not sharp for $a=2,b=6$, due to an exceptional $\mathbb Q$-linear relation $\int_0^\infty [I_0(t)]^2[K_0(t)]^6 t\mathrm{d}\, t=72\int_0^\infty [I_0(t)]^2[K_0(t)]^6 t^{3}\mathrm{d}\, t$, which is provable by integrating modular forms.