Quantum Complexity as Hydrodynamics

As a new step towards defining complexity for quantum field theories, we map Nielsen operator complexity for $SU(N)$ gates to two-dimensional hydrodynamics. We develop a tractable large $N$ limit that leads to regular geometries on the manifold of unitaries as $N$ is taken to infinity. To achieve this, we introduce a basis of non-commutative plane waves for the $\mathfrak{su}(N)$ algebra and define a metric with polynomial penalty factors. Through the Euler-Arnold approach we identify incompressible inviscid hydrodynamics on the two-torus as a novel effective theory of large-qudit operator complexity. For large $N$, our cost function captures two essential properties of holographic complexity measures: ergodicity and conjugate points.

As a new step towards defining complexity for quantum field theories, we map Nielsen operator complexity for SU (N ) gates to two-dimensional hydrodynamics. We develop a tractable large N limit that leads to regular geometries on the manifold of unitaries as N is taken to infinity. To achieve this, we introduce a basis of non-commutative plane waves for the su(N ) algebra and define a metric with polynomial penalty factors. Through the Euler-Arnold approach we identify incompressible inviscid hydrodynamics on the two-torus as a novel effective theory of large-qudit operator complexity. For large N , our cost function captures two essential properties of holographic complexity measures: ergodicity and conjugate points.
1. Introduction.-Quantum computational complexity [1], referred to as complexity hereafter, quantifies the number of simple gates required to synthesize a given unitary operation in quantum computing. In recent years, the geometric approach to complexity of Nielsen et al. [2][3][4] has proven immensely useful for investigating the complexity of n-qubit systems, independently of the particular state of the system. In its original manifestation, this geometric framework relied on the SU (2 n ) manifold of unitaries acting on n-qubits. In this Letter, we consider a generalization of Nielsen's approach that incorporates quantum circuits acting on N -level systems, i.e. qudits of dimension N . For qudit systems, the unitaries of interest belong to SU (N ). The key ingredient in Nielsen's approach is the choice of metric on this manifold, assigning penalty factors to unitaries departing from the identity operator I. The complexity of U ∈ SU (N ) is then identified with the length of the minimal geodesic connecting I and U . These unitaries are generated by a control Hamiltonian H tangent to SU (N ). The corresponding group algebra characterizes fully these Hamiltonians and their geodesics, through the Euler-Arnold method [5][6][7][8].
Recent progress on complexity measures via the AdS/CFT correspondence (or holography) [9][10][11] further motivates the present work. In general, holographic complexity measures should be ergodic and exhibit conjugate points: Ergodicity ensures all points on the group manifold can be reached in finite time, thus implying a linear [12] growth of complexity with time [13,14]. In contrast, conjugate points, i.e. the meeting points of equal-length geodesics, provide bounds on this growth [4]. This conjecture is supported by work on the curvatures of complexity measures in holographic CFTs [7,8].
In this Letter, we show how a judicious choice of basis and metric for the su(N ) algebra with polynomial penalty factors leads to well-defined non-singular geometry for SU (N → ∞) [26]. We show that, on the SU (N ) manifold at infinite N and at low energies, Nielsen complexity can be equivalently evaluated on the manifold of volume-preserving diffeomorphisms SDiff(T 2 ) of the torus. We show that the Euler-Arnold equation on the resulting manifold coincides with the Euler equation of a two-dimensional ideal fluid [27]. This permits identifying control Hamiltonians with diffeomorphism generators in two-dimensional hydrodynamics. Moreover, it suggests a natural cost function, based on the two-dimensional Laplacian, with a smooth dependence on N . This smooth dependence suggests the prevalence of particular characteristics of the hydrodynamic theory even at finite N . In particular, two-dimensional ideal hydrodynamics is classically chaotic due to the hyperbolic geometry of its phase space [28]. We quantify this instability at finite, large N by numerically computing the sectional and Ricci curvatures of SU (N ).
Finally, we find that our formulation of complexity exhibits both ergodicity and conjugate points, whose presence is necessary for a proper holographic complexity measure (although see [26]). In this way, our results arXiv:2109.01152v2 [hep-th] 22 Feb 2022 constitute a new step towards understanding quantum complexity in QFT's with a holographic dual.
2. The algebra of N -level qudits.-At the heart of our construction lies a new and non-trivial choice of anti-Hermitian generators for the su(N ) Lie algebra. We employ known results for su(N ) to investigate the large N limit of our basis and explain the subtleties it carries. The corresponding structure constants determine the Riemannian curvature of SU (N ), which we compute in section 5. We outline our construction below [29].
We first introduce the N × N "shift" matrix h ij = δ i+1,j and "clock" matrix g ij = ω j δ i,j , with ω = exp 2πi N a primitive N th root of unity, and i, j = 0, . . . , N − 1 mod N . These matrices commute up to a phase, i.e. hg = ωgh. Then, following [30][31][32][33], we define a basis of unitary, but not necessarily anti-Hermitian, matrices J m = ω m 1 m 2 2 g m1 h m2 , indexed by a two-vector m = (m 1 , m 2 ) on the Z 2 -lattice. These can be thought of as a non-commutative version of plane waves, with the vector index m playing the role of the wave vector, and h and g the momentum and position modes, respectively [32]. Their commutator is given by where m × n ≡ m 1 n 2 − m 2 n 1 . In [32] it was shown that the algebra (1) is, in the N → ∞ limit, isomorphic to the algebra SVect(T 2 ) of the group SDiff(T 2 ) of volume-preserving diffeomorphisms on the standard twotorus T 2 . To see this, note that SVect(T 2 ) admits a symplectic structure in terms of divergence-free vector fields which, in two-dimensions, are Hamiltonian vector fields X f . The X f are uniquely determined by their associated stream function f [34] , which can be expanded in plane waves on T 2 as f m ∝ exp(i(m 1 x + m 2 p)) [28]. The isomorphism then between SVect(T 2 ) and su(N ) in the large N limit is obtained by expanding the sine in (1) to first order in 1/N and identifying The isomorphism (2) thus relates Hamiltonian vector fields (i.e. elements of SVect(T 2 )) with the basis elements of su(N ) given by (1). A detail not addressed in [30][31][32][33], but already mentioned in [35], is that the Taylor expansion truncation is invalid for several classes of vectors m, n; There are vector pairs defined for all N , e.g. m = ( N −1 2 , 0) and n = (0, N −1 2 ), for which the cross product is of order O(N 2 ) or O(N ). We must restrict the isomorphism to only the pairs of O(1). These turn out to be precisely the low-momentum modes relevant for hydrodynamics, see Sec. 4 and [29] for more details.
We now proceed with the definition of anti-Hermitian basis elements, capable of constructing SU (N ) operators through exponentiation, as required by Nielsen's ap-proach. To this end, we introduce for each m, The generators (3) obey commutation relations inherited from (1) [29]. The structure constants thus obtained are more involved than those in (1) due to the overcompleteness of (3). However, we can make this basis complete via modularity symmetries and linear dependences enjoyed by the generators [29]. Most importantly, due to linearity, the isomorphism (2) carries over to this basis, which hence exhibits a well-defined large N limit. We exploit this in our curvature computations in section 5.
3. Euler-Arnold framework.-Nielsen's approach identifies complexity with the length of the minimal geodesic U (s), with s parametrizing the position along the trajectory [36], connecting the identity element with the desired U on the manifold of unitaries. Each geodesic U (s) is generated by an (anti-Hermitian) control Hamiltonian H(s) via the Schrödinger equation dU ds = H(s)U (s). The Euler-Arnold formalism [28,37] exploits the group structure by identifying H with a Lie algebra element H(s) = U −1 (s)U (s) ∈ su(N ), i.e. with the pullback of the vectorU (s) onto the tangent space at the identity. The time-evolution of H(s) within this approach is then given by the Euler-Arnold (EA) equatioṅ where κ is a quadratic bilinear two-form defined via [X, Y ], Z = κ(Z, X), Y , for X, Y, Z ∈ su(N ) and ·, · the Lie algebra inner product [37]. Distinct κ-forms are induced by different inner products on the algebra. This is equivalent to choosing a metric on SU (N ), and hence a cost function in Nielsen's setup. Using the EA equation is advantageous, since it can be easier to solve (4) than to compute the nested commutators appearing in a solution of the Schrödinger equation [19]. In fact, we explain in the next section how the EA equation drastically simplifies in the large N limit, which allows for the direct calculation of the control Hamiltonian.
4. Inner product and penalty factors.-We now show how at large N , su(N ) leads within the EA framework to the ideal-fluid Euler equation. Solutions to this equation are control Hamiltonians in the sense of Nielsen, which allows for a hydrodynamic interpretation of the standard computation of Nielsen complexity [29]. We also discuss in detail how this suggests a natural extrapolation of the hydrodynamic cost function to finite N .
In the large N limit and for low-energy O(1) generators, SU (N ) is identified via (2) with the manifold of volume-preserving diffeomorphisms SDiff(T 2 ) [38]. We consider the standard inner product on its algebra SVect(T 2 ), given by the L 2 -inner product between Hamiltonian vector fields. This can be rewritten in terms of the Laplacian acting on stream functions as with ω the symplectic form on T 2 [29]. For f = g, this is (twice) the kinetic energy of a flow. The innerproduct (5) induces a metric on SDiff(T 2 ) which defines the length, and thus the cost, of geodesics. Consequently, it defines a κ form given by κ(f, f ) = −∆ −1 {f, ∆f }, with {·, ·} the usual Poisson bracket [29]. This κ form leads to the following EA equation for the control Hamiltonian H with h the "control stream function" associated to H via the symplectic form [29]. Equation (6) constitutes a main result of this work. Considering the large N limit of the su(N ) algebra, we obtain the stream function form of the Euler equation for a (2 + 1)-dimensional ideal fluid [28].
We find that the Nielsen complexity of a large-qudit unitary can be evaluated via this effective hydrodynamic theory, for which the control Hamiltonian can be straightforwardly computed (6). Moreover, the computation of Nielsen complexity can be recast in the hydrodynamic setting: The Schrödinger equation defining the target unitary U in terms of the control Hamiltonian, is now the equation relating the Eulerian and Lagrangian frames of reference of the fluid df ds = H(s, f (s)) [39]. H(s) is the control Hamiltonian obtained from the solution h of (6) and f the corresponding diffeomorphism. Imposing the boundary conditions f (0) = id, the identity map, and f (1) = f target , the target diffeomorphism [40], yields initial velocities v m ≡ v m (s, f target ) [41] for the geodesic as functions of f target . These are inserted into the length functional = 1 0 ds G m n v m v n . Here,G m n is the metric induced by the inner product (5). Minimizing this functional over all solutions (6) yields the complexity of f target , C(f target ). In summary, the Nielsen complexity of large-qudit unitaries is given by the length of the minimal geodesic, generated by a solution to the EA equation (6), that connects the Lie algebra SVect(T 2 ) with the desired target element of SDiff(T 2 ). Our construction, hence, provides a smooth geometry at large N , the manifold SDiff(T 2 ), on which complexity can be calculated. Thus, we avoid the singular geometries encountered in previous manifestations of large N complexity models. Additionally, SDiff(T 2 ) has a clear physical interpretation as the phase space of a well-known theory, namely two-dimensional hydrodynamics.
Motivated by our construction for complexity at N → ∞, and exploiting the generator isomorphism (2), we now formulate new results for the complexity geometry at finite N by ensuring a smooth transition between the finite-and infinite-dimensional setups. We adapt the inner product (5) at N → ∞ to finite N by defining the action of the Laplacian on non-commutative waves as ∆J m = −m 2 J m . This choice transitions smoothly to infinite N , where the action of the standard Laplacian on plane waves f k is ∆ f k = (∂ 2 x + ∂ 2 p ) f k = −k 2 f k [42]. We define an inner product on su(N ) for finite N as with T ∈ {C, S}. By means of group translation, this inner product induces a right-invariant metric G m n on SU (N ). Its components are the penalty factors for different directions on the tangent space, given by the eigenvalues of the Laplacian acting on the generators, i.e. G m n = m 2 δ m n , with m = | m|. These penalties render the metric homogeneous but not isotropic, since not all directions get penalized equally. Equal penalty factors are assigned only to those vectors related by parity or conjugation e.g. m = (1, 2), n = (2, 1) and l = (−2, −1). This reflects the Hamiltonian structure of the problem, with the Laplacian being invariant under symplectic transformations, and is visually manifest in our curvature results shown in Fig. 1. Due to this symmetry, every direction on the Lie algebra gets assigned a different penalty with at most eight-fold degeneracy, yielding a maximally anisotropic metric for the manifold of unitaries. This is an essential property, since anisotropy leads to negative curvature on the manifold and negative curvature is a strong indicator of ergodic geodesic flow [15]. In terms of N -level qudits, our choice of metric ensures high-energy excitations with large wave vector receive larger penalty factors. These high-energy sectors effectively decouple in the strict large N limit.
Our choice of penalty factors is fundamentally different from the majority of previous work on the subject, e.g. [4,14,16,17]. In particular, the penalty factors in these setups grow exponentially p ∼ α k ∼ e k ln α instead of polynomially [43]. Here, k is the Pauli weight of the many-qubit gate in the Pauli basis [44]. Although exponential penalty factors are well motivated from the point of view of local quantum operations, they typically lead to singular geometries in the N = 2 n → ∞ limit [16][17][18]. Instead, our penalty factors remain finite in the large N limit, by transitioning to the position and momentum modes of Hamiltonian vector fields on T 2 via (2). This identification with the hydrodynamical phasespace is naturally restricted to the low-energy sector of O(1) vectors, the so-called admissible directions. Thus, there are no infinitely penalized directions on SDiff(T 2 ), resulting in a regular geometry. The remaining directions of the large SU (N → ∞) manifold containing the high-energy sectors with O(N ) vectors are inadmissible and effectively decouple from the geometry since they are assigned penalty factors that are at least infinite. This situation is captured in the framework of sub-Riemannian geometry [45], also recently mentioned in the context of complexity in [46]. A fundamental theorem due to Chow and Rashevskii [47,48] asserts that geodesics can still reach every point by only accessing admissible directions [49] [50]. That is, trajectories from the hydrodynamic phase space can still reach every large qudit unitary. Since it is infinitely expensive to move in inadmissible directions, the hydrodynamic trajectories have an overall smaller cost, i.e. smaller complexity, cf. [29].
5. Average Ricci curvature-Curvature computations in recent literature regarding Nielsen complexity geometries [14][15][16][17][18]25] typically focus on the sectional curvatures of the manifold of unitaries. The sign of the sectional curvatures is an indicator for convergence (positive sign) or divergence (negative sign) of nearby geodesics [51]. However, the stability of a geodesic and, hence, the emergence of ergodic behaviour does not only depend on the sign of the sectional curvature in the direction parallel to its velocity, but rather on the sign of the sectional curvatures of all two-planes containing its velocity vector [28]. For this reason, we believe that a more telling quantity to describe the stability of a geodesic with velocity vector v is given by the normalized Ricci curvature [52], with K the sectional curvature tensor and m running over the algebra directions. Eq. (8) can be thought of as an average sectional curvature across an orthonormal basis for the tangent space and is well-defined as N → ∞. Importantly, Ric(v) ≤ 0 for SDiff(T 2 ) [52], which is the quantitative reason for the chaotic behavior of two-dimensional hydrodynamics. The smooth limit of our su(N ) basis for large N , given by SVect(T 2 ), indicates we can compute Ricci curvatures of SU (N ) at large N and compare to the hydrodynamic result. We evaluate the Ricci curvature for every direction m in su(N ) for odd values of N ∈ [3,39] by first calculating the corresponding sectional curvatures [29]. Note that, since the dimensionality of the tangent space grows with N , a given velocity m can be defined only after it appears within the distribution of directions at N = N 0 ( m). Its corresponding Ricci curvature Ric( m) is thus defined only after N = N 0 ( m) and will continue to change with N as more and more directions contribute to the average in (8). Our numerical data shows that Ricci curvatures of newly introduced directions at a given N = N 0 are always positive, but all eventually turn negative at some critical value N = N c ( m). The resulting data for N c as a function of the direction m ∈ su(N ) is shown in Fig. 1 and constitutes a second main result of our work. We interpret this figure and its extrapolation at large N as a visual definition of the O(1) sub-sector (the blue region) from which the hydrodynamic theory emerges in the strict large N limit [29].
Our results have the following implications for the complexity geometry of SU (N ) at large N . The large N geometry of the low-energy sector has negative Ricci curvature, thus numerically confirming previous mathematical results [52], as well as indicating emergent chaotic behaviour [53]. This implies geodesics can reach every point of the manifold in finite time, resulting in an ergodic geodesic flow. Therefore, our canonical cost function has a property characteristic of any proper holographic complexity measure as conjectured by [13,14].
Furthermore, the numerically-evaluated distribution of sectional curvatures of our model always contains positive sectional curvatures at large N [29]. The presence of positively sectional curvatures is also a necessary property of our cost function from two points of view: First, while a strictly negative geometry is indeed ergodic [17,54], complexity metrics on Lie groups without positive curvatures are necessarily flat [18,55]. Our choice of metric thus combines the negative average Ricci curvature beneficial for ergodicity with the necessity of having positively curved directions. Second, strictly negative geometries lack an important feature of geometric complexity, namely conjugate points [56]. These are points on a manifold where a geodesic ceases to be globally minimizing, e.g. antipodal points on the sphere. Conjugate points seem to be a necessary feature of complexity geometries, since they forbid complexity from exhibiting an unbounded linear growth with time [4,[16][17][18]. It is wellknown that SDiff(T 2 ), our effective geometry at N → ∞, indeed exhibits conjugate points [28,57]. All in all, our results indicate that even though the large N limit con-sidered here is more similar to the vector large N limit than to the matrix large N limit relevant in holography [58], our setup still exhibits desirable properties of a holographic complexity measure for large enough values of N .
Finally, our results for the sectional curvatures [29] indicate the existence of a universality class of SU (N ) metrics, as defined in [46], indexed by N 2 . These SU (N ) metrics are conjectured to be equivalent, i.e. leading to the same complexity, at late geodesic times. We infer from this conjecture that the complexity of SU (N ) scales at large N and at large geodesic distances as C(SU (N )) C(SDiff(T 2 )) + O(1/N ). Moreover, this implies a finite critical N = N C such that C(SU (N C )) at short distances equals C(SDiff(T 2 )) at long distances, even if the manifolds are not isomorphic.
6. Discussion and Outlook.-For the first time, we provide a definition for Nielsen operator complexity of SU (N ) with a well-defined large N limit. This is realized by using the ideal hydrodynamics equation as the geodesic equation on the low-momentum sector of SU (N → ∞) ∼ = SDiff(T 2 ). The natural choice of cost function is the kinetic energy of the fluid, which we derive within the Euler-Arnold approach and adapt for every value of N . Our construction provides a simple way of computing the control Hamiltonian, thus simplifying one of the main obstacles in the computation of Nielsen complexity. In particular, our setup allows to reach every point of the large SU (N → ∞) manifold via admissible geodesics within the hydrodynamical phase space, thus drastically simplifying the complexity geometry.
From the perspective of quantum information, we find a basis for qudit unitaries that scales nicely with the qudit size. This scalability allows for the synthesis of large qudit gates as long as they only implement O(1) transitions, with respect to N . This corresponds to a locality property of our basis, implying its usefulness for constructing qudit lattices. This is particularly interesting in view of computing complexity of fault-tolerant quantum error-correction qudit architectures [59][60][61].
It is possible to include 1/N corrections for large, but finite, qudit unitaries in order to confirm the conjectured scaling of C(SU (N )) with N in terms of C(SDiff(T 2 )). This will also allow to derive the critical N C at which the complexities for the two manifolds coincide. This approach is closely related to integrable systems in noncommutative geometry [35,[62][63][64], see [29] for a first step in this direction.
Finally, our cost function captures two essential properties in view of holographic complexity, ergodicity and conjugate points. Both are consequences of the phase space geometry of hydrodynamics. This suggests that cost functions based on the Laplacian acting on infinitesimal gates are a promising new avenue for describing operator complexity also in holographic CFTs.
Acknowledgements. We thank Stefan Waldmann and Knut Hüper for useful discussions. We thank R. Auzzi  In this section we explain in detail the construction of a basis for the su(N ) algebra out of the shift and clock matrices h and g, introduced in the beginning of section 2 of the main text. These have the matrix form where we recall that ω = exp 2π N . It is readily confirmed that g and h satisfy a modularity condition g N = h N = 1, and commute up to a phase factor, hg = ωgh. The nomenclature for these matrices comes from their action on momentum and position eigenstates in a finite-dimensional Hilbert space. The "shift" matrix h shifts us from one momentum eigenvector to the next one, while the "clock" matrix g acts on position eigenstates by multiplying with a phase. This phase is a power of a primitive root of unity, all of which can be visualized as points on the complex unit circle, hence the allusion to a clock. In the continuum limit N → ∞, these matrices become precisely the position and momentum modes of a plane wave. As introduced in [30][31][32][33], we can use h and g as building blocks to define a basis of unitary, but not necessarily Hermitian, matrices indexed by a two-vector m = (m 1 , m 2 ) on the Z 2 -lattice, We may regard these matrices as non-commutative plane waves, with the vector index m playing the role of the wave vector. One can readily verify that the J's are traceless and satisfy J † m = J − m . The N -modularity of h and g is inherited by the J matrices and allows the identification of generators on the Z 2 lattice via suitable N -translations with t = (t 1 , t 2 ) ∈ Z 2 being an integer vector. This divides the Z 2 -lattice into N × N cells and allows restriction to generators in what we will denote as the fundamental cell, defined by m 1 , m 2 = − N 2 , · · · , N 2 + 1 for even N and for N odd. Note that the odd N case is technically easier to implement due to the symmetry of the fundamental cell around the origin. All numerical computations were thus performed for odd N , but the setup is just as valid for even N . The generators J (0 mod N,0 mod N ) are proportional to the identity and thus decouple from the algebra. The matrices (S.2) then close under multiplication to the algebra provided in eq. (1) of the main text where m × n ≡ m 1 n 2 − m 2 n 1 .
As they were first introduced in [30][31][32][33], the matrices (S.2) are a unitary basis for the su(N ) algebra thought of as a vector space. This construction, however, neglects the relation between Lie algebra elements and elements of the Lie group. In order for a basis of su(N ) to also act as generators for unitary elements of the Lie group SU (N ), it needs N=5 ℤ 2 to be a set of N 2 − 1 traceless anti-Hermitian generators [65]. To achieve this, we define new generators via the two possible anti-Hermitian linear combinations of J's, one for each vector-index m, The notation should be suggestive of the usual Euler decomposition for the Cosine and Sine functions. However, we have now doubled the number of elements in the basis to 2(N 2 − 1), since every vector in the fundamental cell is associated to both a C and an S matrix. Moreover, some of these C's and S's are linearly dependent.
To construct the reduced set of N 2 −1 linearly independent generators, which we denote as the final distribution, we must associate each point in the (N − 1) × (N − 1) fundamental cell on the Z 2 -lattice to an anti-Hermitian generator. In general, we have two possible choices for this generator for each vector m in the Z 2 lattice, given by C m and S m . Starting from the unique distribution for N = 2, one can systematically find the distribution for the (N − 1) × (N − 1) fundamental cell for N > 2 by filling the new points of the cell with linearly independent generators, as will be described in the following. Let us first consider the case N = 2. By construction, )} vanish identically for all even N . Thus for N = 2, we have only C's as generators and there is no ambiguity in the distribution within the fundamental cell. For N = 3, the fundamental cell extends from −1 to 1 in both directions and there are non-vanishing S generators. To insure compatibility [66] with the N = 2 case, we appoint a C generator to the vectors (0, 1), (1, 0) and (1, 1) of the final distribution. We then fill the remaining points with S generators. However, not all of the remaining points in the 2 × 2 lattice around the origin host an independent S m generator. This is due to the symmetry conditions We denote (S.6) as the linear dependence condition (LDC). According to the LDC, S (1,−1) ∝ S (−1,1) . Thus, the point (1, −1) must be associated to a C generator, which is guaranteed to be linearly independent from the S by construction. Once the first two distributions are set, this procedure can be performed iteratively to find the final distribution for (N + 1). This construction ensures compatibility with the previous distribution and it completes the remaining points with linearly independent generators according to the LDC symmetry (S.6). In the general odd N case, then, we assign C generators to all points on the right half-plane together with the positive m 2 -axis and S generators to all points on the left half-plane and the negative m 2 -axis. A precise list of the points corresponding to C's and those corresponding to S's can be generated automatically for any given N > 3 (even or odd). As an example, the final distribution for the case N = 5 is shown in Fig. 2.

II. Computation of Structure Constants
Once the final distribution has been set, as described in the previous section, we can derive the algebra fulfilled by the C's and S's directly from the commutation relation for the J's in (S.4). These are given by The mixed commutators in (S.7) seem to not be antisymmetric, due to the relative minus sign inside the bracket on the right-hand side of (S.7). However, the subtle substitution procedure giving rise to the correct structure constants explained below also ensures that [T m , T n ] = −[T n , T m ], with T either C or S, depending on m. Therefore, let us now explain how to compute this (N 2 − 1) × (N 2 − 1) × (N 2 − 1) anti-symmetric array of structure constants α m n l of the algebra in (S.7). There are four possible commutation relations in (S.7); which one of these is to be used is uniquely determined by the vector indices m and n. Recall that each m and n is uniquely associated to either a C or an S generator. For example, consider a fixed Z 2 vector m. Then, either C m or S m will be in the final distribution, but not both. It is practical to define the overall prefactor which is always a part of the structure constants. The subscript (±, ±) denotes which of the four possible sign combinations (overall sign and sign in front of the generator on the RHS of (S.7), respectively) is to be used. The correct one depends on the commutator and the structure constant we wish to compute. For example, for the m − n structure constant of the [C, C] commutator, the correct sign combination would be β (+,−) . Further, the commutator uniquely defines the class of generators on the right-hand side of (S.7). The correct structure constants are found when the right-hand side is a linear combination of the right class of generators and these are in the final distribution. For example, the right hand side of a [C, C] commutator must be a function of only C generators in the final distribution. However, the m ± n indices on the right-hand side of (S.7) might lie outside of the fundamental cell, in which case one has to substitute them on the right-hand side of the commutation relations with their equivalent generators from the final distribution. To relate the generators outside the unit cell to those within it, we use the N -modularity of the J's which is inherited by the anti-Hermitian generators as for T ∈ {C, S, J}. We denote these equations as the modularity conditions (MC). Thus, the generators corresponding to integral points on the lattice that can be connected via N -translations can also be related. The precise translation required is controlled by the t ∈ Z 2 vector and depends on the specific quadrant of the Z 2 -lattice that one considers. A case-by-case distinction is tedious but can be automatized. Therefore, the overall factor introduced by eq. (S.9) when replacing generators on the right-hand side of the commutation relations is part of the corresponding structure constant α m n( m± n) . Note that the point reached after applying the MC condition might not correspond to a generator in the final distribution. In that case, one needs to further apply the LDC condition given in (S.6) in order to replace it by the correct generator. This substitution will also carry a prefactor which is to be included into the structure constant. Both MC and LDC conditions are needed to find the final distribution. They are also needed for the computation of the correct structure constants. At first glance it might appear as if (S.6) is just a special case of (S.9), but these are indeed independent conditions originating from two distinct symmetries: Eq. (S.9) imposes the modular symmetry of the vectors on the Z 2 -lattice, while (S.6) is a result of having two copies of the original set of generators.
The overall procedure for calculating the structure constants can be summarized as follows: 1. For given input vectors m and n, find the associated generators according to the final distribution and select the corresponding commutator (S.7). This also defines the class of generators that should appear on the right-hand side of (S.7). Compute the corresponding generator for m ± n.
2. If m± n lies inside the fundamental cell and it corresponds to the right class of generators, the structure constant is simply given by β (±,±) (S.8).
3. If m ± n lies inside the fundamental cell but (its corresponding generator is) not in the final distribution, apply the LDC condition (S.6) and the structure constant is then given by β (±,±) times the prefactor from the LDC condition.
4. Finally, if m ± n does not lie in the fundamental cell, apply the MC condition (S.9) to find the corresponding point k in the fundamental cell. Repeat steps 2 and 3 with the resulting k. The structure constant is then given by β (±,±) times the total prefactor from the substitutions.

III. Curvature Calculations
In this section, we elaborate on the computations of sectional and normalized Ricci curvature based on our choice of basis and inner product for the su(N ). We explain how the sectional curvatures were evaluated in our setup and why a further quantity, the normalized Ricci curvature, is a better quantitative description of geodesic stability.

Sectional Curvature
Sectional curvature K m n is the scalar curvature of a manifold along the direction of the two-plane spanned by two tangent vectors m, n and is defined as [28]  Based on our finite-dimensional setup, we evaluated (S.11) for odd values of N ∈ [2,39]. The value N = 2 is simple enough to be carried out analytically and was thus taken into account. The resulting percentages of non-trivially flat [67], negatively and positively curved directions between the generators (S.5) are shown in Fig. 3. We find that the percentages seemingly stabilize for large N to values of 4 %, 44 %, and 52 % respectively. Although it seems as if the flat directions will vanish, this is just an artifact. There will always be non-trivial vanishing structure constants, and thus flat sectional curvatures, since for all N > 2 there will be at least N − 1 vectors in the unit cell which are linearly dependent (but whose corresponding generators are not) that yield identically vanishing structure constants due to m × n. Note that the result for the single-qubit case N = 2, where there are only three non-trivial directions, we find a ratio of 1 3 : 2 3 of negative to positive sectional curvatures, which coincides with previous analysis of complexity for other cost functions [17].
The interest of previous works in sectional curvatures of complexity metrics lies in their role as indicator for geodesic instability. In particular, if the majority of sectional curvatures in the tangent space is negative, nearby geodesics 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37  have a high probability of exponentially diverging and thus lead to classically chaotic trajectories. Unfortunately, our numerical data on the distribution of the sectional curvatures on SU (N ) with respect to our complexity metric, though robust, is inconclusive as to whether typical trajectories on the manifold of unitaries will behave chaotically. This is not in contradiction with the results of Arnold [28], who showed that (S.19) evaluated on SDiff(T 2 ) hosts a large class of solutions with negative sectional curvature. To see this, we note that while Arnold concluded SDiff(T 2 ) is, in his words, "mostly negative" [28], he does not provide a quantitative description of the sectional curvatures or a rigorous, measure-theoretic description of his results.

Normalized Ricci curvature
The present subsection is concerned with evaluating the normalized or average Ricci curvature introduced in section 5 of the main text, which we find to be a better measure of stability than the sectional curvature alone.
Using our results for the sectional curvatures, we use Eq. (8) of the main text to compute the normalized Ricci curvature for every direction in the su(N ) algebra and track its evolution with N up to N = 39. Through this calculation, we are able to calculate the critical N c of each direction, at which the normalized Ricci curvature becomes negative. The results are shown in Fig. 1 of the main text.
This figure can be used to divide the generators into O(1) and O(N ) subsectors. We observe from Fig. 1 of the main text, that the value of N c seems to grow very slowly with the distance of the vector to the origin. However, it remains finite. This means that every possible direction on the su(N ) algebra for finite N will initially have a positive normalized Ricci curvature, but will eventually turn negative at a higher, finite value of N . For any finite N , there will always be a next higher N at which new directions will enter the final distribution and will have a positive Ricci curvature. These belong to the O(N ) subsector. Simultaneously, there will be directions whose curvature will have already reached the critical value and have turned negative. These will belong to the O(1) sector. However, in the strict large N limit, a transition to an infinite-dimensional manifold takes place. This SU (∞) manifold will have an O(1) subsector, which contains all the directions that can be identified isomorphically with Hamiltonian vector fields. Indeed, this is the entirety of the hydrodynamic phase-space, corresponding to all vectors on Z 2 .

IV. The Large N Limit and Sub-Riemannian geometry
As mentioned in the main text, the isomorphism identifying generators of su(N ) and Hamiltonian vector fields at large N is subtle. Some families of vectors have cross products of order O(N ) or O(N 2 ) large enough to impede a truncation of the sine to first order. In these cases, an identification with Hamiltonian vector fields is not proper. However, for all vectors whose cross product is of order O(1), the argument of the sine is indeed small and a truncation of the Taylor expansion can be made. This restricts the implementation of the isomorphism to only those vectors in this O(1) subsector.
In terms of qudit unitary gates, our basis allows the implementation of any gate describing the transition between two energy levels inside the qudit, cf. fig. 4. However, the restriction to low-momentum O(1) vectors means that the large N limit, and thus, the hydrodynamic description, impose a certain notion of locality. Though it is not a hard cut-off locality prescription, such as the notion of k-locality often imposed when considering complexity of many-qubit gates in the Pauli basis, this constraint indicates that gates implementing O(N ) are highly penalized while O(1) gates are preferred in terms of cost. Of course, the precise meaning of O(1) and O(N ) changes for every finite value of N . Some vectors which are in the O(N ) sector for a given N will be contained in the O(1) sector for a higher N , which gives rise to the evolution of Ricci curvatures explained in section III of the supplemental and section 5 of the main text.
The situation explained above can be described within the framework of sub-Riemannian geometry [45]. The relevance of sub-Riemannian geometry for complexity was originally mentioned by Nielsen [2] and has been recently invoked again by Susskind and collaborators [46]. Our setup can be thought of as defining a family of metrics {(g µν ) N } for every finite N . Each of these metrics will impose a weak locality condition by applying larger penalties to infinitesimal gates in O(N ) directions. The transition to the strict large N limit carries a transition to infinitedimensional manifold with it. In this limit, the family of smooth, approximately sub-Riemannian metrics defined above converge to a proper sub-Riemannian metric on the SU (∞) manifold, In particular, this indicates that hydrodynamics emerges in the strict large N limit as the low-energy, O(1) sub-Riemannian manifold of the much larger SU (N → ∞) group manifold. The metric G µν on the latter consists of two objects: First, the Laplacian, which applies a penalty m 2 in all directions m corresponding to the hydrodynamical phase-space. Note that this encompasses all finite penalties available on Z 2 . Second, a hard-wall metric, imposing infinite penalties in all remaining directions which arose from the strict large N limit of the O(N ) sector. These directions are denoted as inadmissible.
The notion of controllability is fundamental to Nielsen's geometric approach to complexity. It implies that every point of the manifold can be reached via a geodesic following the admissible directions defined by the cost function. In this regard, a main result of sub-Riemannian geometry is the Chow-Rashevskii theorem [47,48], which guarantees controllability on finite-dimensional manifolds provided that the basis for the tangent space is a bracket-generating distribution. This means that every direction in the tangent bundle can be equivalently described by a finite chain of commutators involving only generators associated to admissible directions. For any finite N , our basis (S.5) fulfills the bracket-generating condition, even when restricting to the O(1) subsector. This is because su(N ) is a simple Lie algebra. However, the proper sub-Riemannian metric G µν is defined only in the strict large N limit. Fortunately, an analogue of the Chow-Rashevskii theorem on infinite-dimensional manifolds exists [49]. This suggests that controllablity survives the transition to infinite N in our construction. In other words, we can reach any point of the large SU (∞) manifold via geodesics that flow only through admissible directions given by the hydrodynamic phase-space. Thus, the computation of Nielsen complexity as explained in the main text and in section VI proves to be well-posed.
The growth of the penalty factors with N , implies that it becomes increasingly more "expensive" to move directly in an O(N ) direction, as we veer further and further away from the identity operator, in terms of geodesic distance on SU (N ). Instead, the system will prefer to move indirectly towards an O(N ) direction through the "cheaper" O(1) directions [68]. This heuristic understanding is made more quantitative via the recent conjecture of [46]. This conjecture involves a family of metrics, of the same dimension, {(g µν ) I }, indexed by the penalty factor of the hard directions, I [69]. Namely, the conjecture states that for large values of I and long geodesic distances The index I can also be implicitly defined from the leading order behaviour of the sectional curvatures for large I, schematically K ∼ I. In our case, we can confirm both through the metric and our numerical calculation of the sectional curvatures that, I is identified with N 2 . With this identification, Eq. (S.13) leads to the relation between C(SU (N )) and C(SDiff(T 2 )) mentioned at the end of section 5 in the main text.
Note that, we use the conjecture of Susskind et al. [46] for a family of metrics of different dimension. This is not a violation of the conjecture, but rather a confirmation of the expectations of the authors as seen in the discussion section of [46].
A further aspect of this conjecture is the existence of a critical I, for which the long geodesic distance behaviour matches the short one. At this critical value, it is in general easier to calculate complexity and confirm explicitly the leading scaling behaviour (S.13). On that note, let us expand on the comment made in the conclusions section of the main text regarding an explicit derivation of (S.13). A rigorous derivation of such scalings requires a more detailed analysis of the isomorphism relating su(N → ∞) and SVect(T 2 ). This is achieved by considering a non-commutative version of plane waves. The mathematical framework formalizing this problem is that of non-commutative geometry (NCG) [62]. The standard procedure is to introduce the so-called Weyl-Moyal star product [64], This defines a way of multiplying classical smooth functions in a non-commutative manner. The deformation parameter is θ = 1/N in our case. This leads to a deformed Poisson bracket, the Moyal bracket {f, g} , which in turn defines non-commutative algebras for vector fields. This also yields a non-commutative version of the Euler-Arnold equation In order to properly investigate the scaling of Nielsen complexity on SU (N ) in terms of that on SDiff(T 2 ), we have to determine higher orders in the algebra isomorphism given in eq. (2) of the main text. These are encoded in the expansion of the Weyl-Moyal product, which provides the higher-order in 1/N corrections to the commutative Euler-Arnold equation given in eq. (6) of the main text. A proper identification of finite-and infinite-dimensional algebras at higher orders in 1/N would provide a scaling relation of the type discussed here. There have also been attempts of discretizing the Laplacian in non-canonical ways as a way to define finite-dimensional approximations to two-dimensional hydrodynamics [35], which would imply the definition of a new cost function in Nielsen's approach. We leave such investigations for future work.

V. Derivation of the Euler-Arnold equation for hydrodynamics
The derivation of the Euler-Arnold equation for the hydrodynamics system given by eq. (6) in the main text relies on the following steps. We use Arnold's conventions [28] and define Hamiltonian vector fields as ι X f ω = −df , with ω the symplectic form on T 2 and f the Hamiltonian function dual to X f . d is the exterior derivative and ι denotes the interior product. The symplectic form further defines the Poisson bracket of Hamiltonian functions as {f, g} = ω(X f , X g ). This implies that the commutator of Hamiltonian vector fields is itself a Hamiltonian vector field. This new vector field has Hamiltonian function given by the Poisson bracket of the Hamiltonian functions f and g, Here, L Y h = Y (h) denotes the Lie derivative of a smooth function h in the direction of the vector field Y and we made use of Cartan's magic formula L X = d ι X + ι X d. We now use relation (S. 16), together with the definition of the L 2 -inner product in eq. (5) of the main text, to compute the κ form (S.17) The first equality in (S.17) implements (S. 16). The insertion of the definition of the L 2 -inner product then yields the second line. To arrive to the third line, we make use of the Leibniz rule for the Poisson bracket and the fact that {f, g}ω can be rewritten as the exterior derivative of a 1-form d(g ι X f ω) via Cartan's formula, for any Hamiltonian functions f and g. Then, the generalized Stokes' theorem says that M dα = ∂M α for any 1-form α. For the torus, ∂M = 0 and, hence the integral over α vanishes identically. In the second integral of the third line of (S.17), we insert an identity operator in the form of ∆∆ −1 , to recover the expression for the L 2 -inner product. From this we can read off the resulting κ form which directly determines the EA equation given by eq. (6) in the main text.

VI. Control Hamiltonians, Physical Hamiltonians and Euler-Arnold equations
We now detail the precise role of the control Hamiltonian in the computation of complexity and the technical difficulties that arise in its calculation. First, we clarify some notation: even though Nielsen's approach is valid for general unitaries, a usual choice for the target unitary U in the context of holography is the time evolution operator U = e −iHt of the system under some physical Hamiltonian H up to time t. Control Hamiltonians H(s) and physical Hamiltonians H are distinct objects; the former generates a trajectory, while the latter specifies the boundary condition to this trajectory, meaning that we impose U (1) = e −iHt on the geodesics generated by some H. We note that the geodesic parameter s is not necessarily identical with the physical time t. The computation of complexity can be generally divided into two main steps: the computation of possible control Hamiltonians and the implementation of boundary conditions to find those that actually generate the desired unitary, i.e. solving the Schrödinger equation. We now review these steps.
Within the Euler-Arnold approach [28,37] introduced in the main text, the defining equation for the control Hamiltonian is the so-called Euler-Arnold equation [3], with the quadratic bilinear form κ defined on the Lie algebra as 20) and ·, · the inner product on the algebra. Solutions to (S. 19) in principle yield direct expressions for the control Hamiltonian. However, such a direct computation of the Hamiltonian is in general rather involved, since the κ form with respect to metrics other than the canonical Killing metric yield complicated differential equations. Moreover, the subsequent evaluation of the Schrödinger equation contains expressions with nested commutators [19], which are generally difficult to evaluate. For this reason, the usual procedure is to change perspective and expand the control Hamiltonian in a basis, in order to obtain more tractable equations for the expansion coefficients. On the manifold of unitaries, a chart is provided by a basis of the tangent space, given by anti-Hermitian generators. Due to the group structure of the manifold, every tangent space is isomorphic to the Lie algebra, and thus we can expand the control Hamiltonian in terms of our generators of su(N ), These velocities obey a geodesic equation with respect to the metric associated to the cost function on SU (N ), i.e. our inner product defined in eq. (7) of the main text. Solutions will yield the initial velocities as a set of integration constants. These have to be inserted into the Schrödinger equation and the boundary conditions need to be imposed. This way we can relate the initial velocities to the actual target unitary, i.e. we can find Y m = Y m (s, U ). In the case of unitary time-evolution, this step will relate the velocities to the couplings in the physical Hamiltonian of the system, which can also be expanded in terms of the generators of the SU (N ) group as H = m A m T m . In the end, one wishes to have a set of functions Y m (s, A k ) which define the correct initial velocities of the geodesic. In summary, the overall procedure for computing quantum complexityá la Nielsen consists of the following steps.
1. Choose a cost function F which will define a metric on the manifold of unitaries, thus also defining geodesic paths with respect to it. In our case this is given by the inner product on su(N ) given by eq. (7)  3. Insert these control Hamiltonians into the Schrödinger equation dU ds = HU and enforce the correct boundary conditions, i.e. U (0) = 1 and U (1) = U target . This will fix the velocities of the geodesic, through which the length can be calculated via the standard length functional of Riemannian geometry with respect to the chosen metric.
4. If multiple control Hamiltonians are found, one has to choose the one which minimizes the length, and thus the cost. This yields the complexity C.
The benefit arising from our analysis is that, in the large N limit, we identify the su(N ) algebra with the algebra of volume-preserving diffeomorphisms SVect(T 2 ) via the isomorphism between generators in eq. (2) of the main text. This naturally carries over to any other tangent space, such that we can, locally and at low energies, identify the SU (N ) manifold with the manifold of volume-preserving diffeomorphisms SDiff(T 2 ). Our results show how to define the general setup for Nielsen complexity explained above in such a way that it survives the large N limit. The whole construction can be interpreted in terms of objects in two-dimensional hydrodynamics. Most importantly, our setup simplifies the first step of the computation considerably, providing a tractable and well-known Euler-Arnold equation for the control Hamiltonian, namely the Euler equation for ideal hydrodynamics in (2+1) dimensions given in eq. (6) of the main text. The control stream function h and the corresponding control Hamiltonian vector field H are related in the usual way, ι H ω = −dh. It is important to note that both the control stream function and vector field have an implicit dependence (x, p) on the coordinates of the fluid on the torus. Solutions of this equation can then uniquely be decomposed in order to find the velocities with which the length functional can be evaluated. Naturally, the correct boundary conditions need to be imposed. The Schrödinger equation in the hydrodynamic setup corresponds to the mapping between Lagrangian and Eulerian frames of reference for the fluid [39],

Physical Hamiltonian
We now discuss some details on the role of physical Hamiltonians in Nielsen complexity. Even though the approach explained above and our results are valid for any target operation U , the literature on holographic complexity generally focuses on the target unitary being the time evolution operator. In recent years it has been found [70] that entanglement is not enough to describe all the late time properties of a holographic QFT and the inner black hole geometry dual to it. Thus, complexity was proposed [13] as a suitable additional measure of the information content of the QFT at late times. To focus on the late time behavior, one considers the unitary operation to be the time evolution operator U = e −iHt of the system at time t under some physical Hamiltonian H and calculates its complexity. This physical Hamiltonian is usually assumed to be k-local or chaotic. Thus, we detail some properties that such a Hamiltonian should have if it ought to be compatible with our setup. The physical system we consider is given by a single N -level qudit, see Fig. 4. Notions of locality such as k-locality may be interpreted in this system as the type of internal transitions between the levels that are allowed by the Hamiltonian, e.g. only transitions with energy smaller than some threshold k. However, we take a different approach. We allow the physical Hamiltonian to contain, at first, all possible couplings A m between levels, i.e. H qudit implements all-to-all transitions on the qudit. Written in terms of generators of the algebra in the final distribution, we have T ∈ {C, S} .

(S.24)
A natural assumption would be that there exists no hierarchy between the couplings, i.e. all couplings are assumed to be of order O(1). It is worth noting that although we do not impose any locality restrictions on the Hamiltonian itself, our choice of cost function has an approximate locality property built into the penalties, which is responsible for the restriction to the O(1) sector.
Similarly to the previous discussion for the qudit Hamiltonian, we can interpret these assumptions for volumepreserving diffeomorphisms. More precisely, we can choose the target diffeomorphism f target introduced in the previous section to be the time evolution of the fluid under some physical Hamiltonian H fluid after some time t. We assume this Hamiltonian to obey similar properties as H qudit , albeit adapted to the infinite-dimensional Hilbert space. In particular, this means there are infinite many levels, each one excited by a particular plane wave with wave vector m. The components of this wave vector remain integer, as they have to be compatible with the periodicity of the torus. This Hamiltonian is still assumed to implement all-to-all interactions. Expanding again in terms of the generators, this time of SDiff(T 2 ), we have the schematic structure 25) with the couplings B f indexed by the corresponding plane wave f m with wave vector m. For our analysis, we do not assume a direct connection between A m and B f , instead leaving them free, subject only to the general assumptions explained above.