Let $N\geq1$ be an integer and let $H:[0,1]^2\to\mathbb C^{N\times N}$ be a pointwise hermitean matrix valued function. For $y\in[0,1]$ and $0\leq a\leq b\leq 1$, let $U_y(b,a)$ be the time evolution operator (solution operator) to the Schrödinger equation $i\psi'(x)=H(x,y)\psi(x)$. More precisely, the operator is defined so that for any $\phi\in\mathbb C^N$ and $a\in[0,1]$ the solution of the initial value problem $$ \begin{cases} i\psi'(x)=H(x,y)\psi(x)\\ \psi(a)=\phi \end{cases} $$ is $\psi(x)=U_y(x,a)\phi$. (Alternatively, we could define $U_y(b,a)$ to be the solution of $U_y(a,a)=I$, $\partial_aU_y(a,b)=-iH(a,y)U_y(a,b)$ and $\partial_bU_y(a,b)=iU_y(a,b)H(b,y)$.)
I'm interested in differentiating this time evolution operator with respect to the parameter $y$. A sketchy calculation and physical intuition suggest that $$ \partial_yU_y(1,0) = -i\int_0^1U_y(1,x)\partial_yH(x,y)U_y(x,0)dx. $$ What regularity of $H$ do I need to assume for this to be true? Is there a reference for this result? I feel that giving a proof instead of a reference would be a distraction from what I want (readers of what I'm writing) to focus on, unless there is a particularly short and neat proof.