Since the RHS is $2 \pi$-periodic in all variables, one can consider it on the three-dimensional torus $(\mathbb{R}/2 \pi \mathbb{Z})^3$.

Assume $a_{21} = a_{31}$. Then the two-dimensional torus $$ T := \{\, (x_1, x_2, x_2): x_1, x_2 \in \mathbb{R}/2 \pi \mathbb{Z} \,\} $$ is an invariant submanifold. To investigate its stability, I propose to use the function $u := x_2 - x_3$.

By subtracting the third equation from the second and performing some transformations we obtain (if my computations are O.K.) $$ \tag{1} u'(t) = - 2 a_{21} \cos(x_1(t) - \tfrac{1}{2}(x_2(t) + x_3(t))) \sin(\tfrac{1}{2}u(t)) - (a_{23} + a_{32}) \sin(u(t)). $$ The zero stationary solution of (1) corresponds to $T$. Formally linearizing (1) along the zero solution we obtain $$ \tag{2} v'(t) = - a_{21} \cos(x_1(t) - \tfrac{1}{2}(x_2(t) + x_3(t))) v(t) - (a_{23} + a_{32}) v(t). $$ In the general (multidimensional) case a sufficient condition for the global asymptotic stability of the zero solution of (1) is the existence of $c \ge 1$ and $\lambda > 0$ such that $$ \tag{3} \lVert \Phi(t, s) \rVert \le c \exp(- \lambda (t -s)), \quad s \le t, $$ where $\Phi(t,s)$ denotes the solution operator (transition or Cauchy matrix). (Caveat: in the multidimensional case, for $y'(t) = A(t) y(t)$, the property that the real parts of the eigenvalues of $A(t)$ are, for all $t \in \mathbb{R}$, less than some negative number is not a sufficient condition for (3).) However, in the one-dimensional case we can make use of the comparison property of solutions to an ODE: $-1 \le \cos(x_1(t) - \tfrac{1}{2}(x_2(t) + x_3(t))) \le 1$, so we can compare equation (2) with the autonomous equation $$ w'(t) = (\lvert a_{21} \rvert - (a_{23} + a_{32})) w'(t). $$ From the properties of the cosine function it follows that all the estimates appearing above are uniform both in $t$ and w.r.t. any solution $(x_1(\cdot), x_2(\cdot), x_3(\cdot))$.

So, if $$ \lvert a_{21} \rvert < a_{23} + a_{32} $$ then the invariant torus $T$ is uniformly asymptotically (exponentially) stable.

Since the RHS is $2 \pi$-periodic in all variables, one can consider it on the three-dimensional torus $(\mathbb{R}/2 \pi \mathbb{Z})^3$.

Assume $a_{21} = a_{31}$. Then the two-dimensional torus $$ T := \{\, (x_1, x_2, x_2): x_1, x_2 \in \mathbb{R}/2 \pi \mathbb{Z} \,\} $$ is an invariant submanifold. To investigate its stability, I propose to use the function $u := x_2 - x_3$.

By subtracting the third equation from the second and performing some transformations we obtain (if my computations are O.K.) $$ \tag{1} u'(t) = - 2 a_{21} \cos(x_1(t) - \tfrac{1}{2}(x_2(t) + x_3(t))) \sin(\tfrac{1}{2}u(t)) - (a_{23} + a_{32}) \sin(u(t)). $$ The zero stationary solution of (1) corresponds to $T$. Formally linearizing (1) along the zero solution we obtain $$ \tag{2} y'(t) = - a_{21} \cos(x_1(t) - \tfrac{1}{2}(x_2(t) + x_3(t))) y(t) - (a_{23} + a_{32}) y(t). $$ In the general (multidimensional) case a sufficient condition for the global asymptotic stability of the zero solution of (1) is the existence of $c \ge 1$ and $\lambda > 0$ such that $$ \tag{3} \lVert \Phi(t, s) \rVert \le c \exp(- \lambda (t -s)), \quad s \le t, $$ where $\Phi(t,s)$ denotes the solution operator (transition or Cauchy matrix). (Caveat: in the multidimensional case, for $\xi'(t) = A(t) \xi(t)$, the property that the real parts of the eigenvalues of $A(t)$ are, for all $t \in \mathbb{R}$, less than some negative number is not a sufficient condition for (3).) However, in the one-dimensional case we can make use of the comparison property of solutions to an ODE: $-1 \le \cos(x_1(t) - \tfrac{1}{2}(x_2(t) + x_3(t))) \le 1$, so we can compare equation (2) with the autonomous equation $$ z'(t) = (\lvert a_{21} \rvert - (a_{23} + a_{32})) z'(t). $$ From the properties of the cosine function it follows that all the estimates appearing above are uniform both in $t$ and w.r.t. any solution $(x_1(\cdot), x_2(\cdot), x_3(\cdot))$.

So, if $$ \lvert a_{21} \rvert < a_{23} + a_{32} $$ then the invariant torus $T$ is uniformly asymptotically (exponentially) stable.