## 3 Deriving Input-Output Mappings for the Two Classes of Systems
### The Variational Equation Approach
Denote with \(\mathbf{\Sigma}_{\rm B}^{(k)}\) the \(k\)th subsystem whose dynamics are determined by \(\mathbf{x}_{k}\) (for \(k\geqslant 2\)), i.e.,  
\[\mathbf{\Sigma}_{\rm B}^{(k)}:\ \begin{cases}\mathbf{E}\hat{\mathbf{x}}_{k}(t)=\mathbf{A} \mathbf{x}_{k}(t)+\mathbf{N}\mathbf{x}_{k-1}(t)\hat{\mathbf{u}}(t),\\ \mathbf{y}_{k}=\mathbf{C}\mathbf{x}_{k}.\end{cases}  \tag{20}\]Note that \(\mathbf{\Sigma}_{\rm B}^{(1)}\) is indeed a LTI system, while \(\mathbf{\Sigma}_{\rm B}^{(k)}\) is a LTV (linear time-varying) system for all \(k\geq 2\). Next, we continue with splitting the dynamics of the QBTI model in (13), similar to the BTI models above. The following equality holds true:  
\[\mathbf{Q}\big{(}\mathbf{x}(t)\otimes\mathbf{x}(t)\big{)}=\mathbf{Q}\Big{(} \sum_{k=1}^{\infty}\alpha^{k}\mathbf{x}_{k}(t)\Big{)}\otimes\Big{(}\sum_{k=1} ^{\infty}\alpha^{k}\mathbf{x}_{k}(t)\Big{)}=\sum_{k=2}^{\infty}\alpha^{k} \Big{[}\sum_{j=1}^{k-1}\mathbf{Q}\big{(}\mathbf{x}_{j}(t)\otimes\mathbf{x}_{k- j}(t)\big{)}\Big{]}.\]  
By substituting (19) and \(\mathbf{u}(t)=\alpha\hat{\mathbf{u}}(t)\) into the differential equation in (13), we write  
\[ \begin{split}\mathbf{E}\sum_{k=1}^{\infty}\alpha^{k}\dot{ \mathbf{x}}_{k}(t)&=\mathbf{A}\sum_{k=1}^{\infty}\alpha^{k} \mathbf{x}_{k}(t)+\mathbf{Q}\Big{(}\sum_{k=1}^{\infty}\alpha^{k}\mathbf{x}_{ k}(t)\Big{)}\otimes\Big{(}\sum_{k=1}^{\infty}\alpha^{k}\mathbf{x}_{k}(t)\Big{)}\\ &+\mathbf{N}\Big{(}\sum_{k=1}^{\infty}\alpha^{k}\mathbf{x}_{k}(t )\Big{)}\alpha\hat{\mathbf{u}}(t)+\alpha\bar{\mathbf{B}}\hat{\mathbf{u}}(t) \\ &=\sum_{k=1}^{\infty}\alpha^{k}\mathbf{A}\mathbf{x}_{k}(t)+\sum_{ k=2}^{\infty}\alpha^{k}\Big{[}\sum_{j=1}^{k-1}\mathbf{Q}\big{(}\mathbf{x}_{j}(t) \otimes\mathbf{x}_{k-j}(t)\big{)}\Big{]}\\ &+\sum_{k=2}^{\infty}\alpha^{k}\mathbf{N}\mathbf{x}_{k-1}(t) \mathbf{u}(t)+\alpha\bar{\mathbf{B}}\hat{\mathbf{u}}(t).\end{split} \]  
Now, similarly as for the bilinear case, since this differential equation is valid for all \(\alpha\in\mathbb{R}\), coefficients corresponding to the same powers of \(\alpha\) can be equated and doing so for the coefficients of \(\alpha^{k},\;k>1\), we can write  
\[\mathbf{E}\dot{\mathbf{x}}_{k}(t)=\mathbf{A}\mathbf{x}_{k}(t)+\sum_{j=1}^{k-1 }\mathbf{Q}\big{(}\mathbf{x}_{j}(t)\otimes\mathbf{x}_{k-j}(t)\big{)}+\mathbf{ N}\mathbf{x}_{k-1}(t)\hat{\mathbf{u}}(t). \tag{21}\]  
Decouple the differential equation above into infinitely many differential equations as  
\[ \begin{split}\mathbf{E}\dot{x}_{1}(t)&=\mathbf{A} \mathbf{x}_{1}(t)+\mathbf{B}\hat{\mathbf{u}}(t),\\ \mathbf{E}\dot{x}_{2}(t)&=\mathbf{A}\mathbf{x}_{2}(t )+\mathbf{Q}\big{(}\mathbf{x}_{1}(t)\otimes\mathbf{x}_{1}(t)\big{)}+\mathbf{ N}\mathbf{x}_{1}(t)\hat{\mathbf{u}}(t),\\ &\quad\quad\vdots\\ \mathbf{E}\dot{\mathbf{x}}_{k}(t)&=\mathbf{A}\mathbf{x }_{k}(t)+\mathbf{Q}\big{(}\mathbf{x}_{1}(t)\otimes\mathbf{x}_{k-1}(t)+\ldots+ \mathbf{x}_{k-1}(t)\otimes\mathbf{x}_{1}(t)\big{)}+\mathbf{N}\mathbf{x}_{k-1 }(t)\hat{\mathbf{u}}(t),\\ &\quad\quad\vdots\end{split} \]  
so that (19) holds true. Hence, the initial QBTI system is written equivalently as infinitely many coupled LTV subsystems. Denote with \(\mathbf{\Sigma}_{\rm QB}^{(k)}\) the \(k\)th subsystem in variable \(\mathbf{x}_{k}\):  
\[\mathbf{\Sigma}_{\rm QB}^{(k)}:\; \begin{cases}\mathbf{E}\dot{\mathbf{x}}_{k}(t)= \mathbf{A}\mathbf{x}_{k}(t)+\sum_{j=1}^{k-1}\mathbf{Q}\big{(}\mathbf{x}_{j}(t) \otimes\mathbf{x}_{k-j}(t)\big{)}+\mathbf{N}\mathbf{x}_{k-1}(t)\hat{\mathbf{u} }(t),\\ \mathbf{y}_{k}=\mathbf{C}\mathbf{x}_{k}.\end{cases}  \tag{22}\]  
The advantage of the _variational equation approach_ is that the solution \(\mathbf{x}(t)\) can be derived by iteratively solving a series of LTV systems, which are non-linearly coupled. For instance, starting with the first LTI subsystem, we can directly write the solution as \(\mathbf{x}_{1}(t)=\int_{0}^{t}e^{\mathbf{A}(t-\tau)}\mathbf{B}u(\tau)d\tau\) (assuming \(\mathbf{x}_{1}(0)=\mathbf{0}\)).