Lohe Systems of matrix Riccati equations , linear fractional transformations , partial integrability and synchronization

 After publication please use: “This article may be downloaded for personal use only. Any other use requires prior permission of the author and AIP Publishing. This article appeared in (citation of published article) and may be found at (URL/link for published article abstract).  Prior to publication please use: “The following article has been submitted to/accepted by [Name of Journal]. After it is published, it will be found at Link.”  For Creative Commons licensed material, please use: “Copyright (year) Author(s). This article is distributed under a Creative Commons Attribution (CC BY) License.”

 You may deposit the accepted manuscript immediately after acceptance, using the credit line formatting below  You may deposit the VOR 12 months after publication, with the credit line and a link to the VOR on AIP Publishing's site Format for credit lines  After publication please use: "This article may be downloaded for personal use only. Any other use requires prior permission of the author and AIP Publishing. This article appeared in (citation of published article) and may be found at (URL/link for published article abstract).  Prior to publication please use: "The following article has been submitted to/accepted by [Name of Journal]. After it is published, it will be found at Link."

I. INTRODUCTION
Synchronization phenomena in physical, biological and other complex systems have been extensively investigated in a wide variety of contexts [1][2][3][4][5][6] , including opinion dynamics 7 with related concepts of consensus 8 . In typical models the complex system consists of a large number of elements interacting nonlinearly over a network of connections, with each node comprising an independent dynamical system. The nodes become correlated in time as synchronization occurs, leading to the evolution of the system as a collective entity. In large systems with N nodes there are a correspondingly large number of differential equations to be analyzed and solved, with behaviours that are sensitive to underlying parameters such as coupling strengths, phase lag angles, time delays of interactions between the nodes, as well as the effects of dynamic network topologies. The investigation and analysis of these systems and their synchronization properties for large N remains an ongoing challenge, particularly if the underlying parameters take distributed values that vary with time.
The widely-studied Kuramoto-type models of synchronization 3,4,6,9 , with its many extensions and generalizations 5,10 , are useful because they are amenable to both numerical and analytic investigations, while at the same time displaying a diverse range of behaviours.
In particular, a special case of the Kuramoto model with identical frequencies ω, but of arbitrary size N , can be partially integrated by means of the Watanabe-Strogatz (WS) transform [11][12][13][14][15] , which reduces the number of independent equations from N to only 3. This allows for a much more detailed investigation of synchronization properties as a function of the underlying parameters 16,17 .
The purpose of this paper is to generalize such models and their property of partial integrability to systems where the variable Z i at each node i of the network is a rectangular matrix of arbitrary size, by firstly formulating the defining differential equations as a system of cross-coupled matrix Riccati equations with time-dependent coefficients that are independent of the node. Whereas Kuramoto models and their higher-dimensional extensions generally have cubic nonlinearities, in special cases these can be regarded as quadratic nonlinearities with time-dependent coefficients, and therefore take a Riccati form. These equations can be partially integrated by means of linear fractional transformations, leaving a reduced system of a size that depends only on the matrix dimensions, independent of N .
We find that for square matrices and also for unitary rectangular systems, the constants of integration are related to the eigenvalues of matrix cross ratios.
Special cases of interest are firstly square matrix variables Z i which are elements of a classical Lie group, such as GL(d, C) or the compact unitary and orthogonal subgroups, which have previously been investigated in detail [18][19][20][21][22] . It is known that complete synchronization can occur in these models for identical frequency matrices 18,19 , i.e. Z i Z −1 j approaches the identity matrix exponentially quickly for all i, j = 1, . . . N as the system evolves, for restricted initial values. The second case of special interest, also extensively investigated, is that in which the variables are real vectors that lie on the unit sphere in any dimension, for which complete synchronization again occurs 20,[23][24][25][26][27] . Partial integration has previously been performed for these models 28,29 , although the reduced equations differ from those derived here, see Sec. VII C; evidently the method of partial integration is not unique. Partial integation extends to general systems of rectangular matrices under the restricted circumstances mentioned such as identical frequency matrices, and we may then in principal use the reduced equations to derive detailed properties of the trajectories as has been done for the Kuramoto system 17 . Here, we focus on the methods of partial integration and the subsequent reduced equations, along with properties of linear fractional transformations, matrix cross ratios and the constants of motion, and it remains to show how to use the reduced equations to analyze the specific behaviour of the system for various underlying parameters.

A. Matrix Riccati equations
There is an extensive body of literature on Riccati equations for rectangular matrices of arbitrary size p × q, with applications to random processes, optimal control, and diffusion problems, as well as to other engineering science applications such as robust stabilization, and network synthesis [30][31][32][33][34] , see also several reviews [35][36][37][38] . We refer in particular to Ref. 31 for a comprehensive development with further references, although the applications are in these cases to single matrix Riccati equations, sometimes with constant coefficients.
By contrast, we consider cross-coupled Riccati matrix systems of unlimited size N with coefficients of arbitrary time-dependence and these, it appears, have not previously been investigated, either in general or for the specific application to synchronized systems, except for the Kuramoto model 13,39 and the n-sphere generalizations 28,29 . Although the partial integration that we outline in Sec. III has not been previously described in its generality, many known results for a single matrix Riccati equation using for example linear fractional   transformations and matrix cross ratios, together with corresponding group properties and classical domains, apply also to our matrix Riccati systems. This is to be expected in view of the fact that the properties of a single Riccati matrix equation are a special case of those There are four families of classical domains 48 of which type I, which is diffeomorphic to U (p, q)/U (p) × U (q), is relevant to our models of synchronization. The group U (p, q) acts by means of linear fractional transformations, as we show explicitly in Sec. V D.
It remains an open question as to whether the other classical domains also lead to matrix Riccati systems with synchronization properties. We allow trajectories to lie either inside  [11][12][13][14][15] and also preserves the unit disk, and is an example of a linear fractional transformation 13,39 applied to the complex variables z i . The classical domain in this case is the unit disk on which SU (1, 1) acts transitively by means of the Möbius transformation, and trajectories lie either on the unit circle or within the unit disk, depending on the initial values. The cross ratios C ijkl are constants of motion for all distinct nodes i, j, k, l = 1, . . . N and are invariant under linear fractional transformations.

B. Summary
In Sec. II we review the Kuramoto model and the WS transform as the simplest example of the method of partial integration, together with properties of linear fractional transformations (the Möbius map), cross ratios, and the identification of the classical domain as the unit disk. We show by numerical example that the system can synchronize whether trajectories lie either inside the unit disk, or on the boundary, or both. In Sec. III we consider a general Riccati system of N rectangular matrix equations which we partially integrate and although the derivation is elementary (Theorem 2), this leads to a significant reduction to a set of equations which in number is independent of N . In Sec. IV we restrict our attention to square matrix systems of size d × d for which synchronization is known to occur 19 , and develop properties of matrix cross ratios, the eigenvalues of which are constants of motion.
The existence of these nontrivial constants of motion explains why the system is partially integrable. In Sec. V the system is furthered restricted to unitary matrices which we relate to the well-studied hermitean Riccati matrix equation. The matrix cross ratios lead in special cases to constants of motion for synchronization models on the unit sphere (Corollary 8). We show directly in Sec. V D that trajectories lie in the classical domain U (d, d)/U (d) × U (d).
In Sec. VI we further restrict the system to orthogonal matrices, firstly to point out that matrix cross ratios are defined only for even d, and secondly to show that properties of partial integration extend to noncompact orthogonal groups. As an example we consider the Lorentz group SO(1, 1) and the associated hyperbolic Kuramoto equations, which are known to describe synchronization of relativistic systems in Minkowski space 49 . We find a hyperbolic form of the WS transform and derive the relevant classical domain which in this case is unbounded. In Sec. VII we return to the rectangular matrix system, restricted to satisfy unitarity conditions, which are of interest because they include the case of complex vector unknowns which can be regarded as quantum wavefunctions. The matrix equations describe quantum synchronization, which has been well-studied [50][51][52][53][54][55][56] , and we show in particular in Sec. VII B that partial integration formally extends to infinite-dimensional equations which constitute a nonlinear system of Schrödinger equations. For real vectors the system has trajectories which lie on the unit sphere or within the unit ball, and we show by numerical example that synchronization can occur even for the mixed case (Sec. VII C). A point of particular interest is that the method of partial integration in Sec. VII C, derived using linear fractional transformations, differs from that derived specifically for the n-sphere models 28 ; the relationship between these two methods evidently deserves further investigation. We conclude with a summary and final remarks in Sec. VIII.

MODELS
The WS transform was introduced by Watanabe and Strogatz 11,12 in 1993 as a trigonometric substitution in order to solve systems of globally coupled oscillators described by the Kuramoto model, but was later seen to be equivalent to the Möbius group of transformations which map the unit disk D = {z ∈ C : |z| < 1} onto itself 13,15 . We refer to Ref.
where the frequency ω of oscillation is fixed across all nodes, and where we have included multiplicative parameters λ j . One could also include arbitrary node-dependent phase lag angles α j as well as a time delays τ j > 0 for the interacting variables θ j . Although the WS transform applies to such general models, it is nevertheless restricted in that, as well as requiring identical frequencies, the coupling constant κ must be uniform across all nodes, at least in magnitude, and the network topology must be complete, i.e. we allow only global connectivity. While these restrictions are severe, they result in a system which is expressible in Riccati form and can be partially integrated and analyzed in detail 16,17 .
We perform the WS transform on (1) by introducing N angles ξ i (t) and a complex variable z(t) according to then the N equations for ξ i which follow by substitution into (1) are linear and can be integrated directly, leaving only 3 real equations which govern the time evolution of the system.
The mapping (2) is an automorphism of the unit circle which extends to automorphisms of the unit disk, generating the 3-dimensional Möbius group 13,14 . These transformations generalize in higher dimensions to linear fractional transformations which are automorphisms of the group manifolds of the classical Lie groups.

A. Linear fractional transformations and Riccati equations
We rewrite (1) as a system of Riccati equations for z i = e iθ i , but more generally we consider a system in which z i ∈ C evolves from the initial value z i (0) = z 0 i according tȯ where ω, Γ 1 , Γ 2 are complex functions of t, independent of i. The Kuramoto model (1) is obtained by setting z i = e iθ i together with Γ 1 = Γ 2 = ∑ j λ j z j /(2N ), as previously observed 13,39 . The system (3) in this case allows solutions that can lie inside the unit disk as well as on the unit circle, depending on the initial values z 0 i . Synchronization is measured by the order parameter r = | ∑ j z j /N |, with r → 1 for complete synchronization, for which we have |z i − z j | → 0 as t → ∞. Numerically we find that this indeed occurs for random initial values z 0 i , provided that ∑ j λ j > 0 and more generally, complete synchronization occurs provided that ∑ j λ j cos β j > 0, where β j are phase lag angles 16,17 . Rigorous results have been obtained for the standard Kuramoto model in which all trajectories lie on the unit circle 22,57,58 .
If we allow some of the parameters λ i to be negative with ∑ j λ j < 0 then the system can settle into asymptotic configurations consisting either of asynchronous or bipolar states. This is already known to occur for trajectories on the unit circle 17 , and here we provide numerical examples to show that this occurs also for trajectories in the unit disk. In Fig. 1 we choose N = 15 with ω = 0, and plot 12 trajectories in the unit disk, with the remaining 3 confined to the unit circle. The initial values θ 0 i and the parameters λ i ∈ [−1, 1] are generated randomly. The trajectories are shown in red, the initial points in blue, and the final locations for which asymptotic configurations have been attained are marked in black. In (a) we have ∑ j λ j > 0 and complete synchronization of all nodes occurs, i.e. z i → e iϕ for some fixed angle ϕ for all nodes i. For (b) we have ∑ j λ j < 0, and in this case synchronization does not occur, but the nodes are phase-locked into an asynchronous configuration, with some nodes remaining well inside the unit circle.
In (c) we solve the Kuramoto system (1) with distributed frequencies ω i and for ∑ j λ j > 0, and observe that phase-locked synchronization occurs, with the final configuration similar to those well-known for the simple Kuramoto model. The only difference is that here some of the trajectories start inside the unit circle. Although this system is not expressible in the Riccati form (3) and so is not partially integrable, we present this example in order to demonstrate that synchronization can still occur for general extended Kuramoto models whether trajectories start inside or on the unit circle, possibly in mixed form. Now let us briefly review how the system (3) can be partially integrated 13 . Firstly, we observe that the cross ratios are conserved functions of t for any distinct indices i, j, k, l = 1, . . . N and for arbitrary complex functions ω(t), Γ 1 (t), Γ 2 (t), as follows by direct computation using (3). Secondly, these cross ratios are invariant under linear fractional transformations defined by where ad − bc ̸ = 0, and where the complex coefficients a, b, c, d, as well as z i , are timedependent. We can use these linear fractional transformations therefore to transform the system (3) whilst leaving the constants of motion (4) invariant. Define (assuming that ad ̸ = 0, βγ ̸ = 1), then we rewrite (5) in terms of ζ i , β, γ and substitute into (3). The resulting equations are solved by: as is verified by direct substitution. The first two of these are themselves Riccati equations and the last, which is linear, is partially integrated by writing ζ i = ζ 0 i exp iα, where α is a complex function satisfyingα = ω − i γΓ 1 + i βΓ 2 with α(0) = 0. The initial values β 0 , γ 0 for β, γ are arbitrary, and ζ 0 i is obtained from (7) by imposing z i (0) = z 0 i . The system (3) of N complex equations is therefore reduced to just three equations for the complex variables α, β, γ.
A special case of (3) is when ω is real and Γ 1 = Γ 2 , which we refer to as a unitary system (since the variables z i are unitary for all t > 0 provided initially they are unitary), and is related to hermitean Riccati systems as discussed in Sec. V A. In this case we can choose β = γ = z, and so (7) agrees with (2) upon identifying ζ i = e iξ i and z i = e iθ i . Such a system is the Kuramoto model for which , and in this case α is real.
unitary Riccati systems, g preserves the unit circle provided that |g(z i )| = 1 whenever |z i | = 1, which holds provided that which in turn implies that g ∈ U (1, 1). It follows from these relations that for any g ∈ U (1, 1) and so the transformation z i → g(z i ) preserves both the unit disk and the unit circle 48 , and (7) is equivalent to the Möbius transformation 13 . Trajectories z i (t) of the system either lie entirely inside the unit disk or on the unit circle, depending on the initial value z 0 i for each node i.

III. RECTANGULAR MATRIX RICCATI SYSTEMS
We consider a globally connected network of N nodes in which the variable Z i , a p × q complex matrix, is located at the ith node and evolves according to the equation: where the matrix dimensions are given by The coefficient matrices Γ 1 , Γ 2 , Ω 1 , Ω 2 , depend on t, not necessarily as fixed functions of t, but through the time-dependent variables Z j for j = 1, . . . N . The system (9) comprises therefore a set of N cross-coupled nonlinear equations. We specify initial values Z i (0) = Z 0 i , and assume that solutions exist at least locally. Equations (9) reduce to (3) for p = 1 = q, As examples of coefficient matrices which are used in models of synchronization, Ω 1 , Ω 2 are typically constant hermitean matrices, which therefore have real eigenvalues which can be regarded as frequencies of oscillation, and Γ 1 = ∑ N j=1 a j Z j b j where a j , b j are any fixed set of p × p and q × q complex matrices, respectively, together with Γ 2 = ∑ N j=1 c j Z † j d j where c j , d j are also any set of q × q and p × p complex matrices, respectively, where Z † j denotes the hermitean conjugate of Z j . For the case p = q = d of square matrices an alternative we choose H i = Ω 1 to be independent of i, then we obtain (9) with Ω 2 = 0 and Γ 1 = . Solutions for these particular equations exist locally (Proposition 2.1 19 ), and also globally provided that the initial values are suitably restricted.
Under these conditions the solutions completely synchronize, meaning that

A. Linear fractional transformations
In order to partially integrate (9) we define the linear fractional transformations where the dimensions of the complex matrices A, B, C, D are given by A : p×p, B : p×q, C : then successive linear fractional transformations (10), in which are equivalent to a single linear fractional transformation Z i → g 1 g 2 (Z i ), where g 1 g 2 denotes the matrix product for the matrices (11). Hence successive linear fractional transformations correspond to group composition, and so we view (10) as a group transformation with g ∈ GL(p + q, C), provided that g is invertible. We note the following decomposition of g as a product of lower, diagonal, and upper triangular matrices: is nonzero. The inverse transformation g −1 can be deduced from this decomposition, or directly from (10).
Next, we define (generalizing (6)): where the complex matrices ζ i , β, γ are of size p × q, p × q, q × p respectively. The transformation (10) now reads: We substitute for Z i into (9), and therefore replace the N equations for Z i by N +2 equations for the matrices ζ i , β, γ, which leaves the degrees of freedom residing in the two matrix variables β, γ yet to be fixed; only then is ζ i determined uniquely for all i. These extra degrees of freedom arise from the group properties of the linear fractional transformations (10) with respect to GL(p + q, C). If we choose β = 0 = γ, we regain Z i , however we now define β, γ as the solutions of specific equations, from which the equations for ζ i follow: Proof. The proof is by direct substitution. We choose the equation for β by requiring that (9) be satisfied to lowest order in ζ i which, in effect, means we set ζ i → 0 in (14) and hence we simply replace Z i → β in (9), which gives (15). Next, we substitute (14) into (9), and postmultiply both sides by (γζ i + I q ). Inverse matrices are differentiated using the identity for any invertible square matrix A. Then we substitute forβ from (15). The resulting equation can be manipulated into the form: and we then set the quantity in brackets to zero. The remaining terms in this equation lead directly to (17). On back-substituting forζ i into the term in brackets in (18), and setting it to zero, we obtain (16).
Evidently (15,16) are themselves Riccati equations if we regard the coefficients as having a fixed time-dependence. Similarly, (17) is linear in ζ i and is sometimes referred to as a system of Sylvester differential equations (in homogeneous form), see Ref. 31, Chapter 1, . The fact that solutions of the matrix Riccati equation (9), for any fixed i, can be associated with a linear matrix system is well-known, and is sometimes known as Radon's Lemma (see Theorem 3.1.1 31 ). We emphasize, however, that in applications to models of synchronization the coefficient matrices Γ 1 γ + i Ω 1 and Γ 2 β − i Ω 2 are not actually fixed functions of time, but depend on t through the unknowns Z j by means of expressions such as Γ 1 = ∑ N j=1 a j Z j b j , as described above and hence, upon substituting for Z j using (14), are in fact nonlinear functions of β, γ, ζ j . Nevertheless, we may perform N integrations to solve for ζ i by regarding (17) as a linear system. Firstly, we define the p × p complex matrix S and the q × q complex matrix T by means of: which are square matrix equations of size p×p and q ×q respectively. Hence, if we regard the coefficients as fixed functions of t, then S, T are principal matrix solutions which exist over any interval for which the coefficient matrices exist. Despite the fact that (17) is nonlinear in the variables ζ j , j = 1, . . . N , we can now solve for ζ i for each i in terms of S, T as follows: where S, T are determined by (19).

The solution of the N matrix equations (9) with initial values
where β, γ are determined by (15,16) with arbitrary initial values β 0 = β(0), γ 0 = γ(0), and where the initial values ζ 0 i are given in terms of Z 0 i by: Proof.
Because we have introduced two extra variables β, γ, the system is not fully determined until we specify initial values β(0) = β 0 , γ(0) = γ 0 , which then fixes β, γ as the unique solutions of (15,16). The choice of β 0 , γ 0 is arbitrary, reflecting the group transformation properties of the linear fractional transformations (10), although we do not demonstrate this explicitly. A convenient choice is β 0 = 0 = γ 0 , for then we have from (22) ζ 0 i = Z 0 i . In summary, by partial integration we have reduced the N matrix equations (9) to the four matrix equations (15,16,19) for β, γ, S, T , with associated initial values. For unitary systems, which we define in Secs. V A and VII, we may also choose β = γ † , which then leaves only three complex matrix equations to be solved. The partial integration is therefore useful for N ⩾ 4, but is particularly significant for large N ; numerical simulations of synchronization models have been performed, for example, with N > 25, 000 59 . In terms of matrix elements, (9) comprises pqN complex equations, whereas (15,16,19) comprise (p + q) 2 complex equations, independent of N .
In typical models of synchronization Ω 1 , Ω 2 are constant matrices, and Γ 1 is a complex linear combination of all elements Z j , for j = 1, . . . N , i.e.
where a j , b j are a fixed set of p×p and q×q complex matrices, respectively. This expresses Γ 1 in terms of the four unknowns β, γ, S, T , and similarly for Then (15,16,19) together with their corresponding initial values comprises a set of four coupled nonlinear equations for the unknowns β, γ, S, T . Having solved for β, γ, S, T we construct Z i from (21).

SYNCHRONIZATION
The general results of Section III for rectangular matrices apply in particular to square matrices, and so the system can be partially integrated as before, however we now obtain explicit constants as the eigenvalues of cross ratio matrices, which restrict the possible dynamics of the system. In this section, therefore, we investigate in detail the properties of matrix cross ratios together with linear fractional transformations, which have been wellknown since the work of Siegel 45 and Hua 47 . We set p = q = d in (9) with Z i (0) = Z 0 i , and we assume that the unknowns Z i are invertible at all times. A specific model of synchronization which is defined for a classical Lie group G, which could be noncompact 19 , takes the form: where The matrices H i are elements of the Lie algebra of G, and therefore are also elements of GL(d, C), and K is a coupling constant. Solutions of (25) exist locally (Proposition 2.1 19 ), and global existence also follows provided that some a priori conditions are imposed (Proposition 2.2 19 ). If G is compact the equations (25) admit a unique smooth global solution for any initial values Z 0 is constant for all i, j = 1, . . . , N , then it is known that asymptotic phase-locking occurs for sufficiently large K > 0, and for restricted initial conditions (Theorem 5.1 19 ). Further results have recently been obtained for systems with nontrivial topologies 60 .
Of particular interest here is the case of identical oscillators in which H i = H is independent of i, since then the system (25) takes the matrix Riccati form (24) with In this case exponential complete entrainment occurs, i.e. Z i (t)Z −1 j (t) approaches the identity matrix for all i, j = 1, . . . , N exponentially fast as the system evolves (Theorem 4.2 19 ).
As outlined in Section III, equations (25) can be partially integrated and so can be reduced to (15,16,17). A specific example is the well-studied Kuramoto model (1) with identical frequencies, which corresponds to (24) with d = 1 and Z i = e iθ i , together with the identifications (26) with Ω 1 = ω, and can be partially integrated by means of the WS transform as discussed in Sec. II.

A. Matrix cross ratios
For square matrices the partial integration described in Sec. III can be viewed as a consequence of the existence of conserved quantities for the system (24). We define the d × d where Proof. See Appendix A for the proof by direct calculation.
Γ 1 does not appear in the commutator on the RHS of (28), as is evident from the defining equations (24) and the fact that the cross ratios depend only on the differences Z i − Z j .
Less obvious is the fact thatĊ ijkl is also independent of Ω 2 . If we fix the indices i, j, k, l and denote A = C ijkl , H = iZ k Γ 2 + Ω 1 , then (28) then the solution of (28) is Proof. We can regard S k as the principal matrix solution of (29) which exists over any interval for which the coefficient matrix i Ω 1 − Z k Γ 2 exists. We havė Hence C ijkl (t) ∼ C 0 ijkl for all t > 0 and so the eigenvalues of C ijkl are constant in time and are therefore determined by the initial values Z 0 i for the system (24). In particular the trace and determinant of C ijkl are constant. For the Kuramoto model the cross ratios reduce to (4), which are constant for all trajectories whether on or inside the unit circle.
The fact that the eigenvalues of C ijkl are constants of motion constrains the possible trajectories of the system, in particular any two nodes which are initially co-located, i.e.

Theorem 5. Under linear fractional transformations
Proof. See Appendix B. This result is similar to that derived by Zelikin 35 (Section III). Hence

B. Symmetries and equivalences of the matrix cross ratios
There are 4! = 24 permutations of the fixed indices i, j, k, l, however the corresponding cross ratios are either similar to C ijkl , or are explicit functions of C ijkl , and hence do not lead to independent constants of motion. The symmetries and similarities that we derive here generalize known properties for the case d = 1 of complex cross ratios 13 .

Theorem 6. For the matrix cross ratios C ijkl (27), and for any indices
3. Set j = i in 2., use 1., then relabel m → j.
Next we show that the 24 matrix cross ratios comprising C ijkl and its permutations split into 6 equivalence classes with respect to similarity, each containing four elements.

Corollary 7.
We have the following similarities: Proof. The equalities 3.-7. in Theorem 6 express each of C jikl , C ljki , C likj , C jlki , C ilkj as explicit functions of C ijkl , and so these particular cross ratios are not similar to C ijkl . It is sufficient to show that the elements listed in item 1 are all similar to C ijkl , and then the similarities for the other equivalence classes follow from the functional relations of 3.-7. in Theorem 6. From (27): From Z jl = Z il + Z jk − Z ik we obtain: where implies that P + Q − QP ∼ P + Q − P Q, and hence from (31,32) we obtain C ijkl ∼ C klij .
Next, we use the first equivalence C ijkl ∼ C jilk applied to C klij to obtain C klij ∼ C lkji , which completes the proof of item 1. Similarly, we generate the four elements for each of the remaining five equivalence classes.
It is perhaps useful to observe from items 3.-7. in Theorem 6 that C ijkl commutes with each of C jikl , C ljki , C likj , C jlki , C ilkj , but not with any other elements of the corresponding equivalence class.

V. UNITARY GROUP MODELS OF SYNCHRONIZATION
The variables Z i of the matrix models considered in Section IV are elements of the with initial values U i (0) = U 0 i . These models have been extensively investigated 18-21 , both numerically and analytically, and have synchronization properties similar to those of the Kuramoto model. We have chosen the coupling constant κ to be uniform across the network and the d×d hermitean matrix H to be independent of i, in order to allow partial integration of the system. It is usually supposed that H is constant in time, however partial integration can still be performed for any time-dependent H. Provided that the initial matrices U 0 i are unitary, it follows from (33) that U i remains unitary as the system evolves, i.e. all trajectories remain on the group manifold 20 . In the application to quantum mechanics U i can be regarded as the unitary time evolution operator at the ith node, and H is the Hamiltonian at each node of the quantum system, the real eigenvalues of which comprise the energy levels.
We write the system (33) in the forṁ In the form (34) unitarity is still preserved, since it follows from (34) and its hermitean conjugate that d dt Eq. (34) is a special case of the Riccati system (24) with Γ 1 = Γ = Γ † 2 , Ω 1 = −H = Ω † 1 , Ω 2 = 0, and variables Z i = U i . The partial integration outlined in Section III remains valid, but instead of introducing two independent functions β, γ as in (14), we can choose β = γ † and still satisfy the reduced equations (15)(16)(17) in Lemma 1. The equations for γ and ζ i may be integrated as in Theorem 2. The condition Ω 2 = 0 is not necessary in order to partially integrate models of synchronization, see for example Eq. (26) in Ref. 20 and Eq.

A. Hermitean and unitary Riccati matrix systems
We generalize the properties of the system (34) to define a unitary Riccati matrix system as a square matrix system (24) in which These conditions are typical of synchronization models such as the Kuramoto model (1) and more generally (34). A consequence of the defining properties (36) is that U i evolves as a unitary matrix for all t > 0, provided that U i (0) = U 0 i is unitary. We discuss this more generally in Sec. VII, where we extend (36) to rectangular matrices.

B. Cross ratio matrices and constants of motion
Since the system of Riccati equations (34) is a special case of that for square complex matrices Z i in Section IV, we can define matrix cross ratios C ijkl as in (27), then the eigenvalues are constants of motion as stated in Corollary 4. Specifically, we have where From Corollary 7, item 1. we also have C ijkl ∼ C lkji , and therefore C † ijkl ∼ C ijkl , hence the conserved eigenvalues of C ijkl are either real, or appear as complex conjugate pairs. Both the conserved trace and determinant of C ijkl are therefore real.

C. Trajectories on U(2)
The case d = 2 is of particular interest as the simplest example of a non-Abelian model with synchronization properties 18,20,21 that is also partially integrable. We parametrize U i ∈ U (2) according to: where σ k denotes the Pauli matrices, and where , and in order that U i ∈ U (2) we impose the constraint x i x i = 1, hence x i ∈ S 3 . Since H is hermitean we have the expansion H = ν I 2 + ∑ 3 k=1 ω k σ k for frequencies ν and parameters (ω 1 , ω 2 , ω 3 ), and we define the antisymmetric matrix Ω by: The equations of motion (33) or (34) reduce to 20 : with initial values θ i (0) = θ 0 i , x i (0) = x 0 i ∈ S 3 . If we choose θ 0 i = 0 for all i, together with ν = 0, then the unique solution of (42) is θ i = 0 for all t > 0, and so from (41) this implies that U i ∈ SU (2). This reduction from U (2) to SU (2) is possible because the RHS of (33) has zero trace for d = 2 provided that tr H = 0, which ensures that det U i = 1 for all i and all t > 0, provided that det U 0 i = 1. That trajectories on SU (2) can be equivalently formulated as trajectories on S 3 is of course due to the fact that the group manifold of SU (2) is diffeomorphic to S 3 . Complete synchronization for the system (33) of identical oscillators occurs from any initial configuration [18][19][20][21][22] .
In order to explicitly evaluate the constants of motion, namely the determinant and trace of C ijkl for any i, j, k, l = 1, . . . N , we calculate from (41): where we have used x i , x j ∈ S 3 . From (39) we obtain: and we may also evaluate the conserved quantities tr C ijkl . These are not independent constants, but are related to det C ijkl as follows from the identity trA = 1+det A−det(I 2 −A) which is valid for all 2 × 2 matrices A. On setting A = C ijkl and using Theorem 6, item 4., we obtain tr It is evident that the system (42,43) can be extended to any dimension n with trajectories on S n−1 × S 1 , by allowing x i ∈ S n−1 ⊂ R n , and so we deduce: Corollary 8. For the system (42,43) in which x i ∈ S n−1 ⊂ R n , where Ω is any n × n antisymmetric matrix, the parameters λ ijkl = det C ijkl defined by (44) are constants of motion for any i ̸ = l, j ̸ = k.
For θ i = 0 these constants λ ijkl are precisely those found for the unit sphere models, which are partially integrable by means of transformations which preserve the unit sphere 28 .

D. Linear fractional transformations, unitarity and classical domains
We partially integrate the system (34) in any dimension d as in Theorem 1, and hence we substitute (setting β = γ † ) into (34) which is satisfied provided thaṫ The equations for ζ i can be integrated as shown in (20), giving ζ i = S ζ 0 i T −1 , and so we obtain the three matrix equations, namely (46) for γ and (19) for S and T , together with the initial determine the conditions which ensure that U i is unitary, and hence the conditions which are satisfied by the trajectories of the variable γ which appears in (45). We firstly return to the definition (10), U i → g(U i ) = (AU i + B)(CU i + D) −1 and consider conditions on the d × d coefficient matrices A, B, C, D which ensure that g(U i ) is unitary. As is well-known 63 , unitarity is satisfied provided that in which case the 2d × 2d matrix g defined by (11) Hence g ∈ U (d, d) and the mapping U i → g(U i ) = (AU i + B)(CU i + D) −1 preserves unitarity for any square matrices A, B, C, D satisfying (47). It also follows from (47) that is also positive definite, i.e. g preserves the bounded classical domain defined by the condition U † i U i < I d . Next, we define as in (13), ζ i = AU i D −1 , γ = CA −1 , then it follows from (47) that Hence, the d×d matrices ζ i are not themselves unitary unless d = 1. We suppose that I d −γ † γ and I d − γγ † , which are each hermitean, are positive definite matrices; if for example we choose γ(0) = 0, then both I d − γ † γ > 0 and I d − γγ † > 0 hold for all t > 0, and so the square roots of these matrices are well-defined. Then unitary as a consequence of (49). Define also the d × d matrices V, W according to then from (47) both V and W are unitary, and these expressions represent polar decompositions of A, D. Also from (50) we have: From the definition of g we find which defines matrices g γ and g V W . Evidently g V W ∈ U (d) × U (d). By factorizing g γ according to: As a specific example, for d = 1 denote γ = z, and set V = e iψ , W = e iϕ , then from (50,51) we have: Since I d − γ † γ = I d − γγ † = 1 − |z| 2 is positive definite, the trajectories z(t) lie within the unit disk, and the group element is which satisfies g † Jg = J, and so g ∈ U (1, 1). The matrix is an element of SU (1, 1) and acts transitively on the unit disk, which is precisely the classical domain for the Kuramoto model. By means of the matrix g z therefore we can transform any initial value z 0 = γ(0) to lie at the origin. The substitution (45), on setting U i = e iθ i , ζ i = e iξ i , reduces to the Möbius transformation (2).

VI. ORTHOGONAL GROUP MODELS OF SYNCHRONIZATION
The square matrix models considered in Section IV are also of particular interest for the group SO(d), since it is known from numerical and analytic investigations 19,20 that synchronization occurs in such models. The variables R i in this case are d×d real orthogonal matrices with unit determinant which evolve according to (33) . The case d = 2 reduces to the Kuramoto model. From the general considerations of Sect. III, this system can be partially integrated by the substitution where γ, ζ i are real d × d matrices, having set β = γ T .
Matrix cross ratios C ijkl are defined as in Section IV A for complex square matrices, by replacing Z i → R i , except that now the dimension d is restricted to be even, since otherwise R i − R j is singular for any two elements R i , R j ∈ SO(d). This follows from: using det A = det A T for any square matrix A, and also For odd d this implies det(R i −R j ) = 0, and so the matrix cross ratios are singular. For even d, however, the matrices C ijkl evolve as in Theorem 3, and the eigenvalues are constants of motion.

A. Noncompact orthogonal group models
The unitary groups considered in Sect. V are compact and so the trajectories U i (t) are bounded, and synchronization occurs under conditions which generalize those for the Kuramoto model. Synchronization can also occur for noncompact groups in which case trajectories are unbounded, although in some cases solutions exist only locally 19 . We demonstrate here that partial integrability and properties of linear fractional transformations and matrix cross ratios extend to the noncompact case by considering the noncompact group SO(m, n), A special case is n = 0 for which K is the identity matrix, when we regain the SO(d) model. (25) in the Riccati forṁ where the d×d real matrix Ω is an element of the Lie algebra of SO(m, n), and hence satisfies (25). We choose Γ = κ ∑ j λ j M j /(2N ) with multiplicative real constants λ j . We obtain the Riccati system (24)  The system (55) can be partially integrated as described in Section III A, where now we set β = Kγ T K and substitute into (55) to obtain: The equations for ζ i may be partially integrated as in Theorem 2. Let us consider this now in detail for d = 2, for which there are applications to synchronization in relativistic dynamical systems.

B. The Lorentz group SO(1,1) and partial integrability
For d = 2 with n = 1 = m the elements M i of the system are SO(1, 1) matrices which evolve according to (55), and complete synchronization occurs 19 in the sense that M i M −1 j → I 2 as t → ∞. This system describes the synchronization of a relativistic cluster of particles in Minkowski space-time 49 , for distributed matrices Ω i . We consider here in particular identical matrices Ω for which the Riccati system (55) can be partially integrated. We parametrize elements of SO(1, 1) according to for hyperbolic angles α i , hence M i ∈ SO + (1, 1) (the connected component of the identity) and then M T with JK + KJ = 0. Define Ω = ωJ for some real parameter ω, then Ω is an element of the Lie algebra of SO(1, 1) as required. The evolution equations (55), in which we choose which evidently is a hyperbolic form of the Kuramoto model (1). This model has been discussed in the context of relativistic mechanics in 1 + 1 dimensions 49 , in which the particle where α i is the rapidity of the particle of unit mass at the ith node. For identical parameters ω i = ω and with κ i = κ/λ i , complete synchronization occurs exponentially quickly, with α i − α j → 0 as t → ∞.
The system (59) may be partially integrated by substituting for M i as in (56). Because SO(1, 1) is abelian we have ζ i ∈ SO(1, 1) and hence we parametrize ζ i = e β i J for hyperbolic angles β i , similar to (58). Since γ commutes with ζ i for all i we can parametrize γ according to γ = v 1 I 2 − v 2 J for time-dependent real coordinates v 1 , v 2 , and hence γ is a symmetric matrix. The formula (56) reduces to which is a hyperbolic form of the WS transform (2).
The trajectory of the vector (v 1 , v 2 ) ∈ R 2 is restricted to a region in R 2 as determined by the positivity of the LHS of (60). By taking the cases in which e β i is either arbitrarily small or arbitrarily large, positivity requires both v 1 + v 2 ⩾ 0 and v 1 − v 2 ⩾ 0. The inverse mapping to (60) is given by and is of the same form as (60), but exists only if 1 − v 2 1 + v 2 2 ̸ = 0. Hence the trajectories (v 1 , v 2 ) are restricted to the region R in the plane defined by which is the region bounded by the lines v 1 = ±v 2 for v 1 ⩾ 0, and the right branch of the hyperbola v 2 1 − v 2 2 = 1, provided that the initial value (v 1 (0), v 2 (0)) lies in R. The mapping (60) does not give all possible values for e α i as (v 1 , v 2 ) evolves, rather trajectories α i are such that, from (61) Partial integration of the system (59) is achieved by substituting for α i using (60) to obtain equations for v 1 , v 2 , β i , as may also be derived directly from (57). The N equations which are the given initial values for (59). Hence we obtain a closed system of 3 equations The matrix cross ratio C ijkl given by (27) .
These cross ratios are conserved with respect to the evolution equations (59) g γ in the factorization g = g γ g V W is given explicitly by and satisfies g T γ K ′ g γ = K ′ , as well as det g γ = 1, where K ′ = diag[1, −1, −1, 1], and so we deduce that g γ ∈ SO (2,2). Trajectories (v 1 , v 2 ) are therefore constructed on the homogeneous (1,1) which is diffeomorphic to the unbounded classical domain R defined in (62) via the mapping γ → g γ .

VII. RECTANGULAR UNITARY MATRIX RICCATI SYSTEMS
We return now to the general Riccati system of equations (9) for rectangular p×q matrices Z i with dimensions as in Sec. III, but with coefficients restricted according to which extend the conditions for square matrices that define a unitary Riccati system in (36) Writing Γ 1 = Γ, we consider therefore the systeṁ where a j , b j are any set of p × p and q × q complex matrices, respectively, together with constant hermitean matrices Ω 1 = Ω † 1 , Ω 2 = Ω † 2 . Eqs. (65) can be partially integrated as described in Section III with β = γ † , and hence we substitute into (65), which is satisfied provided thaṫ The latter N equations can be integrated as in (20) to give ζ i = S ζ 0 i T −1 , and so we obtain The three matrix equations, (67) for γ and (19) for S with T , together with the initial conditions ζ 0 i = Z 0 i , γ 0 = 0, S(0) = I p , T (0) = I q comprise a system of nonlinear matrix equations, independent of N , which determine γ, S, T , and hence the unknowns Z i for all i = 1, . . . N and any t > 0.
It follows from (65), as shown in (C2,C3) by setting i = j, that d dt If we set Z † i Z i = I q initially, then we must have Z † i Z i = I q for all t > 0, and similarly for Z i Z † i = I p , by uniqueness of solutions. For square matrices this means that if Z i is initially chosen to be unitary for all i, then this unitarity is preserved as the system evolves, as discussed in Sec. V. For q = 1 with Z i a real column p-vector, if we specify Z † i Z i = 1 at t = 0 then Z i remains a unit vector for all t > 0, i.e. the trajectory of the ith node is restricted to the unit sphere.
On the other hand, if we choose I q − Z † i Z i to be positive definite at t = 0, then it also follows from (68)

A. Matrix cross ratios
The partial integrability of the system of square matrix Riccati equations discussed in Sec. IV can be viewed as a consequence of the existence of conserved quantities, namely the eigenvalues of the cross ratio matrices C ijkl defined in (27). Although the expression for C ijkl requires that the matrices Z i be square, we can manipulate C ijkl into a form which, upon replacing the inverse Z −1 i by the hermitean conjugate Z † i , is also well-defined for rectangular matrices, and satisfies conservation properties with respect to the system (65), similar to those of C ijkl . Define therefore the p × p matrix D p ijkl for any i, j, k, l = 1, . . . N : and the q × q matrix (39).

Lemma 9.
The matrix cross ratios D p ijkl , D q ijkl evolve according to: for any i, j, k, l = 1, . . . N .
Proof. See Appendix C.
It follows as in Corollary 4 that the eigenvalues of D p ijkl , D q ijkl are constants of motion. Of the relations and symmetries for C ijkl listed in Theorem 6, those for items 1-3 remain valid for D p ijkl , D q ijkl , and numerical evaluations indicate that the similarities in Corollary 7 also hold.
We consider next the special case q = 1 with p = d for which Z i is a real or complex d-vector. General properties such as partial integrability and the existence of constants of motion remain valid for q = 1, but we consider this particular case here for several reasons; firstly, complex vectors Z i are relevant to quantum synchronization where they can be viewed as wavefunctions which evolve according to finite-dimensional nonlinear Schrödinger equations 64 . The system of equations can be extended to infinite dimensions in which the wavefunctions are elements of a Hilbert space in which the Hamiltonian at each node acts as a self-adjoint operator. It is known 50-56 for models with identical potentials that solutions exist globally, and that complete synchronization occurs under specified conditions 56 . We show that partial integration extends to this infinite-dimensional case, with a fixed number of nodes N . Secondly, we consider in Sec. VII C the formulas for partial integration for models with real d-vectors lying on the unit sphere which have been extensively investigated, mostly without using the properties of partial integration. As we will see, the method of partial integration using linear fractional transformations is not unique, and we briefly compare two alternative schemes.

B. Quantum synchronization and partial integration
We first consider models of quantum synchronization 64 in which a quantum oscillator such as a quantum particle of fixed spin is located at each node i of the network with a corresponding wavefunction |ψ i ⟩. Each d-dimensional wavefunction evolves according to where we have set ℏ = 1 and a ij = 1 (for all-to-all coupling), and where the Hamiltonian is therefore a complex function ψ i (x, t) of both x ∈ R n and t. The Hamiltonian H, which is a sum of kinetic and potential terms, acts in H as a self-adjoint operator with an identical potential for each node, and the time-dependent scalar product ⟨ψ j |ψ i ⟩ is defined in the usual way as a Hilbert space inner product: Equations (72) therefore comprise a system of N coupled partial differential-integral equations. It is known 50-56 that complete synchronization can occur in both finite and infinite dimensions.
We write (72) in the form where the averaged wavefunction |Φ⟩ is defined by |Φ⟩ = κ ∑ j |ψ j ⟩/(2N ). In finite dimensions we identify Z i = |ψ i ⟩, Ω 1 = −H, Ω 2 = 0, Γ = |Φ⟩, and these equations take the form of a unitary Riccati system (65). Partial integration proceeds by means of the linear fractional transformations (66), which in turn leads to the system (67) which may be partially integrated. This procedure formally extends to infinite dimensions in which (73) is regarded as an infinite-dimensional Riccati system. Hence, denoting ζ i = |ζ i ⟩ and γ † = |β⟩, which are now elements of a Hilbert space H and are therefore functions of x ∈ R n and t, we substitute into (65), which is satisfied provided that The equations for |ζ i ⟩ are satisfied by the wavefunctions M t=0 equal to the identity operator. The system is reduced therefore to partial differentialintegral equations for the operator M and for the wavefunction |β⟩ ∈ H. The cross ratios D q=1 ijkl in (71) are constant in time, as follows from Lemma 9, with the explicit form: That these cross ratios are conserved has been independently observed in recent work 56 , where various synchronization and stability properties of the system are also derived but without using the partially integrated (reduced) system (75). Again, our aim in this section is to demonstrate that partial integration extends formally to infinite-dimensional systems.

C. Vector models on the unit sphere
Models of synchronization on the unit sphere S d−1 have been extensively developed, including proofs of synchronization properties together with various applications 20,23-27,65-67 , including consensus properties in opinion dynamics [68][69][70][71] . It has been previously observed that the system with identical frequency matrices can be partially integrated in any dimension d using the vector transform 28 , and also for S 3 by using quaternionic variables 29 . We show here that the general formalism applied to these vector models leads to known properties such as conserved cross ratios, but in particular we observe that the details of partial integration differ from those in Ref. 28. Evidently partial integration can be performed in several distinct ways. Here, we use linear fractional transformations which differ from the vector transform, which algebraically maps the unit sphere to itself in any dimension and also preserves the unit ball.
We choose q = 1, p = d and denote the real column d-vector Z i = x i , then x i ∈ R d satisfies (65) which reads:ẋ where we have set Ω 2 = 0 and have defined the real antisymmetric d × d matrix Ω = i Ω 1 , and have denoted Γ = X, a real column d-vector with Γ † = X T . A typical choice could be where λ j are multiplicative real parameters, but more generally could be any set of constant rotation matrices in SO(d) which implement rotational phase lag.
We choose initial values x i (0) = x 0 i ∈ S d−1 which ensures that x i ∈ S d−1 for all t > 0, but we also allow x 0 i ∈ B d for some or all nodes i, in which case the corresponding trajectories are strictly confined to the unit ball B d . In order to partially integrate (77) by means of the linear fractional transformation (66), we substitute for x i according to: where w = γ T and ζ i = ζ i are real column d-vectors. The equations for w, ζ i which solve (77) are given by (67), which reads: where the d × d matrix M is defined by M = Ω + Xw T − (w X)I d , and so is independent of i. The equations for ζ i can be integrated as before, and we obtain ζ i = Sζ 0 i where the initial values are ζ i (0) = ζ 0 i , and S is the principal matrix solution ofṠ = M S with S(0) = I d . Since S has d 2 independent elements, there are d(d +1) coupled nonlinear equations for w, S which remain to be solved.
By contrast, the vector transform proceeds 28 by means of the mapping u i → x i defined by: which maps the unit sphere to itself and also preserves the unit ball for any v ∈ R d . The defining equations (77) are partially integrated by specifying equations to be satisfied by v and R, where u i = Ru 0 i ; in this case the reduced equations are d(d +1)/2 in number because R ∈ SO(d). The mapping (80) differs from (78), which does not algebraically preserve the unit sphere since ζ i is not a unit vector for arbitrary w ∈ R d , whereas (80) preserves the unit sphere for any v ∈ R d . Nevertheless, because the equation is satisfied, x i remains a unit vector as the system evolves, whether we use (78) or (80), provided that x 0 i is a unit vector. Similarly, if we choose x 0 i ∈ B d , then x i remains inside B d for all t > 0 for that particular node i.
The cross ratios defined in (71) take the explicit form: and are constants of motion as follows from Lemma 9, in agreement with Ref. 28 (equation (27)). These are a special case of the cross ratios discussed in Sec. V C, see in particular the cross ratios (44) which are constants of motion for the system (42,43) as stated in Corollary 8, which holds for any time-dependent coefficient X(t).
Let us briefly discuss synchronization properties of the system (77) as determined numer- where λ j are real parameters. In particular we demonstrate that synchronization can occur even if some trajectories are confined to the unit ball, and others lie on the sphere. For d = 2 it is known 28 that complete synchronization occurs whenever ∑ j λ j > 0, but for ∑ j λ j < 0 either asynchronous states appear in which all asymptotic positions x i are distinct, or else (N − 1, 1) states, sometimes known as bipolar states 22 , occur in which all nodes except one approach the same position. Bipolar states also occur in models of quantum synchronization 56 and are unstable if λ i = 1 for all i. Our numerical studies indicate that similar behaviour occurs for d > 2, and that these properties are maintained whether trajectories lie inside the unit ball or on the unit sphere.
Such behaviour occurs for d = 2 as shown in Fig. 1 (a,b).
In Fig. 2 with d = 3 we choose N = 15, Ω = 0, with randomly generated values for each λ j ∈ [−1, 1], together with randomly generated initial values x 0 i that lie either inside or on S 2 . We plot trajectories (shown in red) in the unit ball B 3 , with just one node confined to the surface S 2 . The initial points are shown in blue, except for the single node on S 2 in green, and the asymptotic (final) locations are marked in black. In (a) we have ∑ j λ j > 0 and complete synchronization occurs with all trajectories converging to a single point which, in the limit, lies on S 2 . For (b) we have ∑ j λ j < 0, and in this case asynchronous states appear in which the asymptotic positions x i are distinct, with all nodes except one on S 2 (in green) remaining well inside the unit ball. The plots Fig. 2 (a,b) are generalizations to d = 3 of those in Fig. 1 (a,b) for d = 2. In (c), for which again ∑ j λ j < 0, we obtain a bipolar state in which all nodes except one completely synchronize to a final position which is arbitrarily close to S 2 . The single node which lies exactly on S 2 , for this example, is diametrically opposite the remaining asymptotic nodes. We have obtained these configurations numerically by solving both the original equations (77) and, as a consistency check, the reduced equations

VIII. CONCLUSION
We have partially integrated a system of rectangular matrix Riccati equations of arbitrary size which in special cases is known to describe the synchronization behaviour of a nonlinear leading to: where we have cancelled multiplicative factors and additive terms wherever possible. In particular, the terms involving Ω 2 cancel out completely. Collecting all terms, as well as subtracting and adding Z ik Z −1 il Z l Γ 2 Z jl Z −1 jk , we obtain: We also have (using which, after substitution into (A1), leads to the required resultĊ ijkl = [C ijkl , Z k Γ 2 − i Ω 1 ].

Appendix B: Transformation of matrix cross ratios
We prove here Theorem 5. We decompose g as shown in (12) and so C ijkl again remains invariant. For B = 0 = C and D = I d (10) reads Z i → AZ i , in which case Z ij → AZ ij and Finally, for A = I d = D and B = 0 (10) reads and so which gives (B4) From the decomposition (12), by applying successive factors from the right, we obtain: which, on substituting for Z ′′ k , correctly reproduces (10), while for C ijkl we find: We can rearrange this expression by means of the Sherman-Morrison-Woodbury formula 74 , which reads: We set U V = (AZ k + B)(CZ k + D) −1 C in this formula to obtain (30).