General formula for mixed models

I'm trying to wrap my head around the general formula of mixed models and how it relates to the system of equations I'm used to.

The general formula read like this:

$$\mathbf{Y_{j}}=\mathbf{X_{j} \beta} + \mathbf{Z_{j}b_{j}} + \epsilon_{j},$$

where:

• $\mathbf{J}$ is a cluster ${j = 1,...,J}$ and for each cluster $j$ there are ${i = 1,...,n_{i}}$ observations
• $\mathbf{Y}$ is the response vector of length $n_{i} \times 1$, where $n_{i}$ is the number of observations.
• $\mathbf{X}$ is the fixed-effects design matrix of dimensions $n_{i} \times p$, where $p$ is the number of fixed-effect parameters.
• $\boldsymbol{\beta}$ is the vector of fixed-effects coefficients of length $p \times 1$.
• $\mathbf{Z}$ is the random-effects design matrix of dimensions $n_{i} \times q$, where $q$ is the number of random effects.
• $\mathbf{b}$ is the vector of random effects of length $q \times 1$.
• $\epsilon$ is the vector of residual errors of length $n_{i} \times 1$, often assumed to follow $\epsilon \sim N(0, \sigma^2 \mathbf{I})$, where $\sigma^2$ is the variance of the residual errors.

However, I'm unsure, how the general formula relates to e.g. a model like this. A simple random slope model with a Level 1 predictor $X$:

L1:

$y_{ij} = \beta_{0j} + \beta_{10} \cdot X_{ij} + \epsilon_{ij} $

L2:

$\beta_{0j} = \gamma_{00} + U_{0j}$

$\beta_{10} = \gamma_{10} + U_{1j} $

For one cluster I would write the same model in the general formula like this. However I'm unsure about the random effects design matrix $\mathbf{Z}$.

\begin{equation} \begin{bmatrix} y_{1j}\\ \vdots \\ y_{n_{i}j} \end{bmatrix} = \begin{bmatrix} 1 & X_{1j} \\ \vdots \\ 1 & X_{n_{i}j} \\ \end{bmatrix} \cdot \begin{bmatrix} \gamma_{00} \\ \gamma_{10} \\ \end{bmatrix} + \begin{bmatrix} 1 & X_{1j}\\ \vdots \\ 1 & X_{n_{i}j} \end{bmatrix} \cdot \begin{bmatrix} U_{0j} \\ U_{1j} \end{bmatrix} + \begin{bmatrix} \epsilon_{1j}\\ \vdots \\ \epsilon_{n_{i}j} \end{bmatrix} \end{equation}

For a more complex model like this - a random slope model with a cross level interaction and a level 2 predictor $M$.

L1:

$y_{ij} = \beta_{0j} + \beta_{10} \cdot X_{ij} + \epsilon_{ij} $

L2:
$\beta_{0j} = \gamma_{00} + \gamma_{01} \cdot M_{j} + U_{0j}$

$\beta_{10} = \gamma_{10} + \gamma_{11} \cdot M_{j} + U_{1j} $

I would write this model for one cluster like so:

\begin{equation} \begin{bmatrix} y_{1j}\\ \vdots \\ y_{n_{i}j} \end{bmatrix} = \begin{bmatrix} 1 & X_{1j} & M_{j} & M_{j} \cdot X_{1j} \\ \vdots \\ 1 & X_{n_{i}j} & M_{j} & M_{j} \cdot X_{n_{i}j} \end{bmatrix} \cdot \begin{bmatrix} \gamma_{00} \\ \gamma_{10} \\ \gamma_{01} \\ \gamma_{11} \end{bmatrix} + \begin{bmatrix} 1 & X_{1j} \\ \vdots \\ 1 & X_{n_{i}j} \end{bmatrix} \cdot \begin{bmatrix} U_{0j} \\ U_{1j} \end{bmatrix} + \begin{bmatrix} \epsilon_{1j}\\ \vdots \\ \epsilon_{n_{i}j} \end{bmatrix} \end{equation}

Is my understanding/usage of $\mathbf{Z}$ correct? Why do I have to specify the $X_{ij}$ again for the $\mathbf{Z}$ matrix?