The matrix regularization for Riemann surfaces with magnetic fluxes

We consider the matrix regularization of fields on a Riemann surface which couple to gauge fields with a nonvanishing magnetic flux. We show that such fields are described as rectangular matrices in the matrix regularization. We construct the matrix regularization explicitly for the case of the sphere and torus based on the Berezin-Toeplitz quantization, and also discuss a possible generalization to cases with higher genera. We also discuss the matrix version of the Laplacian acting on the rectangular matrices. ∗ e-mail address : adachi@het.ph.tsukuba.ac.jp † e-mail address : ishiki@het.ph.tsukuba.ac.jp ‡ e-mail address : takaki@stp.dias.ie § e-mail address : ksaito@het.ph.tsukuba.ac.jp


Introduction
The matrix regularization plays important roles in the matrix-model formulations of M-theory or superstring theory [1,2]. The first quantized theory of a membrane or a string is mapped by the matrix regularization to the matrix model [3], which is conjectured to give a nonperturbative formulation of M-theory or superstring theory. In the matrix regularization, functions on a closed symplectic manifold (M, ω) are linearly mapped to N × N matrices. In this paper, we consider the case that the manifold M is a closed Riemann surface, which is relevant to the regularizations of closed membranes or strings. In this case, the main property of the matrix regularization is that, for any f, g ∈ C ∞ (M), their images T N (f ), T N (g) ∈ M N (C) of the matrix regularization satisfy [4]  where ||·|| is a matrix norm, { , } is the Poisson bracket on M defined by ω and V = M ω/2π is the symplectic volume. The first two properties show that the matrix regularization approximately preserves two algebraic structures of functions associated with the ordinary pointwise product and the Poisson bracket. For Riemann surfaces, the matrix regularization satisfying (1.1) can be constructed by using the Berezin-Toeplitz quantization [5][6][7][8], which we will review later.
In this paper, we consider a generalization of the above setup to Riemann surfaces with nonzero magnetic flux. Suppose that there exists a U(1) magnetic flux on M as M F/2π = Q with Q a non-zero integer. Note that the gauge field A of the field strength F cannot be globally defined, since any globally defined connection would lead to M dA/2π = ∂M A/2π = 0 for a closed manifold. The gauge field A should be defined on each local patch and, on an overlap of any two patches, they are related to each other by gauge transformations. A typical example is given by the Wu-Yang monopole configuration on S 2 , which we will review in later sections. Complex scalar fields coupling to A through the gauge covariant derivative are also defined locally and receive gauge transformations on the overlaps. In the matrix regularization, only globally defined functions are usually considered. We will consider the matrix regularization of locally defined fields, which couple to the gauge field of the nontrivial magnetic flux. This setup will be relevant for describing D-branes in terms of matrices, on which there can exist nontrivial gauge fluxes.
The locally defined scalar fields are mathematically said to be local sections of the complex line bundle on M with the connection A, where the globally defined fields correspond to the special case of the trivial bundle with A = 0. The local sections form a module of the algebra C ∞ (M).
Here, a left module M L of a unital algebra A is an abelian group such that there exists an operation (f g) · a = f · (g · a), for all f, g ∈ A and a, b ∈ M L , where 1 A is the identity element of A. Similarly, the right module can also be defined with the right multiplication. For the case of the local sections of the line bundle, A = C ∞ (M) and multiplying an element of C ∞ (M) to local sections gives the operation ·. In physical terminology, (1.2) is just the property that U(1) charged fields with the same charge form a vector space and a product of a U(1) charged field and a neutral field gives another charged field with the same charge. In this case, the left and right multiplication gives the same operation, so the local sections give a left and right modules of the algebra C ∞ (M).
The Serre-Swan theorem [9] states that vector bundles on M are dual to modules of the corresponding algebra of functions on M. The fuzzy counterpart of this theorem would suggest a correspondence between the fuzzy version of vector bundles and modules of the matrix algebra M N (C). Any module of M N (C) can be written as a set of rectangular matrices 1 . Thus, it is expected that the matrix regularization should be generalized such that the charged scalar fields are mapped to rectangular matrices.
For the fuzzy sphere, there is indeed such a mapping from local sections to rectangular matrices [10][11][12] (see also [13] for the fuzzy CP n ). In [14,15], it is shown that the map can be formally constructed for Kähler manifolds. The main property of this map is that the relation, holds for any function f and local section a of a complex line bundle, where T N N ′ is the linear map from local sections to N × N ′ rectangular matrices. The difference N ′ − N corresponds to the monopole charge (the Chern number) of the line bundle and this should be kept fixed when one takes the large-N limit. The property (1.3) guarantees that the structure of the module (1.2) is approximated well in terms of the rectangular matrices. Note that when N = N ′ , the charge is vanishing and the local sections are just ordinary functions. In this case, T N N ′ reduces to T N and (1.3) just means the first property of (1.1).
In this paper, after presenting a general construction of the maps T N and T N N ′ , we first show that this construction can be embedded in the Berezin-Toeplitz quantization in a U(2) gauge theory. Then, we explicitly demonstrate the construction for the sphere and the torus. In the case of the fuzzy sphere, this construction gives the well-known fuzzy spherical harmonics [10][11][12][13]. For the fuzzy torus, this provides rectangular matrices written in terms of elliptic functions. We will also construct fuzzy versions of the Laplacians, which act on the rectangular matrices and realize the continuum spectra in the commutative limit.
We also discuss the case of Riemann surfaces with higher genera. In this case, we could not explicitly construct the mappings due to some technical difficulties. In particular, we will discuss that obtaining the orthonormal basis of spinors, which is necessary for defining each matrix element of T N and T N N ′ , is technically difficult to compute, though a non-orthonormal basis can be generally written down. Nevertheless, we derive a general form of the Bergman kernel, which formally defines the map T N N ′ . This paper is organized as follows. In Section 2, we review the Berezin-Toeplitz quantization, which gives systematic constructions of T N and T N N ′ . We also discuss that these constructions are unified in the Berezin-Toeplitz quantization in a U(2) gauge theory. In Section 3 and 4, we explicitly construct this mapping for the case of the sphere and the torus, respectively. In Section 5, we discuss the generalization to surfaces with higher genera. In Section 6, we summarize our results and discuss possible applications. In the appendices, we show some details.

Berezin-Toeplitz quantization
In this section, we review the Berezin-Toeplitz quantization and its generalization to rectangular matrices. We also show that the quantizations with square and rectangular matrices can be reformulated in terms of a U(2) gauge theory. In the following, we denote a closed Riemann surface by M.

Quantization for functions
Let us first briefly outline the Berezin-Toeplitz quantization for C ∞ (M). In this quantization, one first needs to construct zero modes of a certain Dirac operator 2 . Let N be the number of independent zero modes and {ψ I |I = 1, 2, · · · , N} be an orthonormal basis of the zero modes. Then, the Berezin-Toeplitz quantization is given by a map where · means the contraction of spinor indices. This map satisfies (1.1), if appropriate geometric quantities are used in the construction of the zero modes or the Dirac operator, as we will explain below 3 . More detailed setup is as follows. Let (g, ω, J) be a Kähler structure on M, which is a compatible triple of a metric, a symplectic form and a complex structure. The surface M has a spin structure associated with J. Let S be a spinor bundle on M. The fiber of S is C 2 and sections of S are spinors with two components. We define a Dirac operator acting on sections of S by where σ a (a = 1, 2) are Pauli matrices, Ω µab and θ µ a are the spin connection and the inverse of the zweibein for the metric g, respectively, and N is a positive integer corresponding to the charge of the spinor fields. We choose the gauge field A to be the symplectic potential, namely, A is given by ω = V dA on each local patch, where V = M ω/2π. The field strength F = dA then satisfies M F/2π = 1. See Appendix A for a detailed definition of the Dirac operator. With this set up, it follows from the index theorem that the number of zero modes of (2.2) is equal to N. Let {ψ I |I = 1, 2, · · · , N} be an orthonormal basis of the zero modes, with respect to the inner product Then, the Toeplitz operator for f ∈ C ∞ (M) is defined by (2.1). It is shown that with this definition, (2.1) satisfies the main properties (1.1) of the matrix regularization [7,8].

Quantization for local sections
The spinor fields of S transform as ψ → e iN α ψ under a gauge transformation, where α is a local gauge parameter. We can consider complex scalar fields with charge Q, which transform as ϕ → e iQα ϕ. In this subsection, we consider a matrix regularization of such charged scalar fields.
One cannot use (2.1) for the charged fields, since ψ † J ·ϕψ I is not gauge invariant. In order to make a gauge invariant mapping, we introduce two copies of the spinor bundles S (N ) and S (N ′ ) which have the same connection A but different charges, N and N ′ , respectively. Let {ψ (N ) I |I = 1, 2, · · · , N} and {ψ (N ′ ) I |I = 1, 2, · · · , N ′ } be orthonormal bases of the Dirac zero modes, which transform as ψ . Then, for a given charged scalar ϕ with charge Q, we define a rectangular matrix In order for this to be gauge invariant, N and N ′ have to be related by N ′ −N = Q. Note that (2.4) reduces to (2.1) for Q = 0. The map (2.4) can also be formulated in terms of the Bergman kernel, which is a projection operator onto the Dirac zero modes. See Appendix B for this formulation. The map (2.4) has the property (1.3) [15]. See Appendix C for the proof. Thus, it indeed gives a natural matrix regularization of the charged scalar fields.

Quantization in U (2) gauge theory
Here, we show that the map (2.4) is naturally obtained as Toeplitz operators in a U(2) gauge theory.
Note that the definition of (2.1) can be generalized such that the gauge group of A is nonabelian [17,20,21]. Suppose that the gauge group is U(2) and S is the spinor bundle in the fundamental representation. We represent the 2-dimensional vector space of the fundamental representation by using the two-component representation where each upper and lower component is a spinor on M. On these spinors, we can consider actions of adjoint scalar fields of U(2), which can be represented as 2 × 2 matrices: We then consider a quantization of the algebra of the adjoint scalars using the Toeplitz operators. In order to realize the mapping (2.4), let us consider a special case where only the connection of a diagonal U(1) subgroup is nontrivial and the Dirac operator takes the form, where A µ is a U(1) connection satisfying M F/2π = 1. Namely, the upper and lower components of the fundamental representation transform as spinors with charge N and N ′ , respectively, under the diagonal U(1) gauge transformation. The spinor bundle is thus decomposed to a direct sum of S (N ) and S (N ′ ) introduced in the previous subsection. From (2.6), we see that the (1, 2) and (2, 1) elements of an adjoint scalar field behave as fields with charge N − N ′ and N ′ − N, respectively, while the diagonal elements behave as neutral fields. A basis of zero modes of the Dirac operator is given by {Ψ I |I = 1, 2, · · · , N + N ′ } with are bases introduced in the previous subsection. For an adjoint scalar field ϕ, the Toeplitz operator is defined by If we consider a complex adjoint scalar with only (2, 1) element, the Toeplitz operator (2.9) is nonvanishing only for I = 1, 2, · · · , N and J = N + 1, N + 2, · · · , N + N ′ . Thus, we obtain the rectangular map (2.4) as a special case of the Berezin-Toeplitz quantization in the U(2) gauge theory. Note that the quantization of the diagonal elements of scalar fields gives the mapping between neutral fields and square matrices. Hence, this formulation using the U(2) gauge theory gives a unified quantization for charged and non-charged fields.

Fuzzy sphere
In this section, we construct the quantization on a sphere S 2 .

Geometry of S 2 with a magnetic flux
We define two open subsets of S 2 by are the polar coordinates and 0 < φ ≤ 2π. We also define the stereographic coordinates on U 1 and U 2 by z = tan(θ/2)e iφ and w = 1/z, respectively. We define a Kähler metric on S 2 by (3.1) The compatible symplectic form is then given by (A.2) and the volume is V = S 2 ω/2π = 1. With respect to this metric, we can choose a Kähler potential, which is defined by (A.7), as From (A.6) and (A.8), we also choose a spin connection and U(1) gauge field defined on U 1 as respectively. Note that we have Ω ab σ a σ b /4 = iAσ 3 in this gauge. This gauge field is known as the Wu-Yang gauge configuration. On the overlap region U 1 ∩ U 2 , the gauge field A transforms as Let ϕ (Q) be a complex scalar field with charge Q coupling to A. Corresponding to the gauge transformation (3.4), ϕ (Q) transforms as lm . Let ψ = (ψ + , ψ − ) be a spinor field on S 2 , where ± stands for the chirality of each component. In our gauge, ψ ± transform as ψ ± (w) = e ±iφ ψ ± (z), This means that ψ ± are charged scalars with charge ∓1.

Dirac zero modes on S 2
We first construct a Dirac operator (2.2) on S 2 and compute its zero modes. The Dirac operator flips the chirality, so that it has only off-diagonal elements in the chiral representation. From (3.2) and (A.9), the Dirac operator is given by [19,22,23] where D ± are the matrix elements of D acting on the spaces with chirality ±1, respectively. In this section, we also write N = 2J + 1 by using a half integer spin J. Note that spinor fields ψ (N ) = (ψ (N )+ , ψ (N )− ) on which D acts transform as See appendix D for the definition of the monopole harmonics.
In the following, we will derive (3.9). For the decomposition ψ = (ψ + , ψ − ), the Dirac equation Dψ = 0 reduces to two differential equations D ± ψ ± = 0. We can easily solve these equations and find that ψ ± = (1 + |z| 2 ) ∓(N ∓1)/2 h ± , where h + and h − are arbitrary holomorphic and antiholomorphic functions on U 1 , respectively. We focus on the norm of ψ ± given by (3.10) We find that the norm of ψ − does not converge for N ≥ 1 unless h − = 0, whereas the norm of ψ + converges when the degree of h + is less than N. Therefore, we find KerD = KerD + = N, which is consistent with the index theorem. Since h + can be expanded in terms of the basis 1, z, z 2 , . . . , z N −1 , we can choose ψ where we used z = tan(θ/2)e iθ in the second equality. By comparing this with the definition of the monopole harmonics (D.4), we finally obtain (3.9).

Berezin-Toeplitz quantization on S 2
Here, we construct the Berezin-Toeplitz map on S 2 and show that (2.4) relates the monopole harmonics and the so-called fuzzy spherical harmonics. By using these Dirac zero modes derived above, we can construct (2.4) for a charged scalar ϕ (Q) as 12) where N ′ − N = Q. Since any charged scalars with charge Q can be expanded in terms of the monopole harmonics Y lm . We introduce a normalization asỸ lm for convenience. Then, we have Here, we used the formula (D.5) in the first line and the symmetric properties for the Clebsch-Gordan coefficients in the second line. Furthermore, by taking into account the parity 4 for the pairs of J and r, J ′ and r ′ and l and Q/2 as well as the relation 2J ′ − 2J = Q, we find that (3.14) Thus, we have whereŶ lm(JJ ′ ) is the fuzzy spherical harmonics [10][11][12][13]26,31,32]. See Appendix E for the definition ofŶ lm(JJ ′ ) . Here, the Clebsch-Gordan coefficient in (3.15) is given by lm ) coincides with the fuzzy spherical harmonics in the large-N limit except for the trivial overall factor (−1) −l .
Note that (3.12) contains the ordinary matrix regularization for functions as a special case with N ′ = N. For example, the Toeplitz operator for the standard embedding function into R 3 , is given by . This is the well-known configuration for the fuzzy sphere [26].

Laplacian on fuzzy S 2
Here, we construct the matrix Laplacian which acts on the rectangular matrices (3.15).
The Laplacian for functions on S 2 is given by the Casimir operator of SU(2). Since we have the representation of SU(2) on the space of charged scalar fields with charge Q as explained in Appendix D, we can naturally define the Laplacian on the fields by Here, L Similarly, there is another representation of SU(2) on the space of rectangular matrices, which is given in Appendix E. From this structure, it is natural to define the matrix Laplacian bŷ • are the representation of the SU(2) generators defined by (E.1). By the definition (E.2), the fuzzy spherical harmonicsŶ lm(JJ ′ ) are the eigenvectors of ∆ and the spectrum is given by Since we have the relation 2J ′ − 2J = Q, the spectrum of∆ coincides with that of ∆ except for the cutoff J + J ′ , which depends on the matrix size and goes to infinity in the large-N limit. From (3.15), we therefore have∆ for l ≤ J + J ′ .

Fuzzy torus
In this section, we construct the quantization on a torus T 2 .

Geometry of T 2 with a magnetic flux
A complex torus is defined as a quotient space of the complex plane, where ∼ stands for the periodic identifications of a discrete lattice: Without loss of generality, the parameter space of τ can be restricted to the fundamental domain with |τ | > 1, Im τ > 0 and −1/2 < Re τ < 1/2 as usual. We express z in terms of two real variables x, y as where x and y are periodic coordinates, for which we take the fundamental region as x, y ∈ [0, 1).
We introduce the Kähler metric on T 2 as According to (A.2), compatible symplectic form is then given by Thus, the symplectic volume will be V = T 2 ω/2π = Im τ /π. From (A.7), the Kähler potential for the metric (4.4) can be chosen by where ζ := ζ 1 + τ ζ 2 and ζ 1 and ζ 2 are real constants corresponding to the gauge holonomies along the 1-cycles on T 2 . From (A.9), the U(1) gauge field is given by The gauge field A is periodic up to a gauge transformation as where λ 1 and λ 2 are given by the Wilson loop phases for x and y directions, respectively: Im(z + ζ), (4.9) The gauge transformation (4.8) is sometimes called the twisted boundary condition. The spin connection on T 2 is evidently zero from (A.6).
We also introduce spinor fields on T 2 . Let ψ (N ) be a spinor field with charge N ∈ Z. We impose the twisted boundary condition as . (4.10) With this boundary condition, it is easy to see that covariant derivatives of ψ (N ) also satisfy the same boundary condition.

Dirac zero modes on T 2
Let us construct a Dirac operator on T 2 . From (4.6) and (A.9), the Dirac operator is given by We decompose a Dirac spinor as ψ (N ) = (ψ (N )+ , ψ (N )− ) and we introduce χ ± (z,z) by Then, the Dirac equation Dψ (N ) = 0 is equivalent to the following equations: The boundary conditions for χ ± are given from (4.10) and (4.12) as (4.14) Below, we solve (4.13) to determine χ ± [27]. Let us first consider χ +5 . The periodicity of χ + along the x-direction enables us to expand it in a Fourier series: (4.15) By substituting (4.15) into (4.13), we obtain the differential equations for the coefficients c n (y), for ∀n ∈ Z, where the prime denotes the y-derivative. The solution to these equations is given by for ∀n ∈ Z, where k n are complex integration constants. By substituting this into (4.15), we obtain Then, we use the boundary condition for the y-direction. From the second equation in (4.14), we obtain the following recursion relation: The solution to this equation is From (4.12), (4.18) and (4.20), we obtain Because of the condition N n = N n+N , ψ (N )+ can be further decomposed into |N| linearly independent solutions: where r = 0, 1, · · · , N − 1 and ϑ is the Jacobi-theta function defined by The negative chirality mode ψ (N )− can be computed in a similar way and is given by From the definition of Jacobi-theta function and the positivity of Im τ , ψ (N )+ r and ψ (N )− r converge only when N > 0 and N < 0, respectively. In the following, we assume that N > 0, so that ψ (N )− = 0 and the zero modes are finally given as given by (4.22). We can determine the normalization factor N r in such a way that ψ (N ) r become orthonormal. The inner product of the zero modes ψ (N ) r is computed in Appendix F and is given by Thus, if we put r |r = 0, 1, · · · , N − 1} forms an orthonormal basis of the zero modes.

Berezin-Toeplitz quantization on T 2
In this subsection, we construct the Toeplitz operator (2.4) for local sections of a complex line bundle on T 2 .
To construct (2.4), we need two copies of spinor bundles with charges N and N ′ . Let {ψ r ′ |r ′ = 0, 1, · · · , N ′ −1} be the orthonormal bases of the Dirac zero modes, each of which has the form (4.25). We consider another scalar field ϕ (Q) with charge Q = N ′ − N, which satisfies the twisted boundary condition, (4.28) The map (2.4) for ϕ (Q) is then given by Note that the integrand in the inner product on the right-hand side is gauge invariant and hence is a completely periodic function on the torus. The map (4.29) reproduces the well-known configuration of the fuzzy torus given by the clockshift matrices. To see this, we consider the quantization of the functions u(x, y) := exp(i2πx) and v(x, y) := exp(i2πy), which are completely periodic and hence have the vanishing charge Q = 0. The Toeplitz operator (4.29) for these functions are given as where C and S are clock and shift matrices respectively: with q = e i2π 1 N . They satisfy the well-known algebra of the fuzzy torus 6 : C N = S N = id (N ) and CS = qSC.
Similarly, for any periodic function with Q = 0, which can be expanded as the Toeplitz operator takes the form of Note that because of the relations U N , V N ∝ id (N ) as well as the orthogonality under the trace, the matrices U n V m with n, m = 0, 1, · · · N − 1 form a complete basis of M N (C). The coefficients f nm in (4.33) are given as follows. By a direct calculation, we can first obtain where l n andñ are the quotient and remainder of n divided by N, and C ln,lm,ñ,m are given by Then, we find that the coefficientsf nm in (4.33) are given bỹ Now, let us consider the quantization for Q = 0. As in the case of the sphere, it is convenient to consider (4.29) for a basis of the local sections. See Appendix G, where we construct such a basis as eigenstates of the Laplacian. In the following, we focus on for the eigenstates ϕ  By using this relation, we write (4.38) as where the inner product in the last line is used in the sense of (G.11). By using (G.11) 7 , we obtain n,r , we finally obtain (4.42) See also [29,30], in which essentially same computations are done in different contexts.
The rectangular matrices (4.42) are the fuzzy version of the eigenstates of the Laplacian. Though they look very complicated, we will show in the following that they give approximate eigenstates of the matrix version of the Laplacian and the spectrum indeed agrees in the large-N limit with that of the continuum Laplacian.

Laplacian on fuzzy T 2
In this subsection, we construct the matrix Laplacian, which acts on the rectangular matrices (4.42). From now on, we put τ = i and ζ = 0 for simplicity.
We first note that there is a crucial difference between the spectrum of the continuum Laplacian for Q = 0 and that for Q = 0. When Q = 0, the Laplacian is given by The spectrum of this operator is just 4π 2 (n 2 + m 2 ), where n and m are integers. For Q = 0, however, as shown in Appendix H, the spectrum becomes that of the 1-dimensional harmonic oscillator (or equivalently the Landau level), because of the relation [D z , Dz] = const. As we will see below, the matrix Laplacian naturally reproduces both of these spectra in the large-N limit. 7 Note that the range of indices in (G.11) can be extended such that (ϕ We first consider the matrix Laplacian for Q = 0, which is relatively well-known. In this case, the continuum Laplacian can be written in terms of a Poisson bracket and we can construct the matrix Laplacian by replacing the Poisson brackets with the commutators of matrices. Let us introduce the Poisson bracket induced from the symplectic form (4.5): where X f is the Hamiltonian vector field of f , namely, it is defined by ω(X f , v) = df (v). The partial derivatives can be expressed in terms of the Poisson bracket as Thus, we can express the Laplacian as  N 2 ), we obtain the following mapping rule for the Poisson bracket: This suggests that a natural choice of the matrix Laplacian iŝ From the algebra of U and V , we can easily prove that U n V m (n, m ∈ N) are eigenstates of∆ as∆ where [n] q := q n/2 − q −n/2 q 1/2 − q −1/2 = sin(nπ/N) sin(π/N) . (4.50) Since |q 1/2 − q −1/2 | → 4π 2 /N 2 and [n] q → n as N → ∞, the spectrum of the matrix Laplacian reduces to∆ in the large-N limit. This agrees with the continuum spectrum. We next construct the matrix Laplacian for Q = 0. A natural generalization of the Laplacian (4.48) for rectangular matrices iŝ where F is an arbitrary N × N ′ rectangular matrix and the operation • is defined by with Toeplitz operators A (N ) and A (N ′ ) with dimension N and N ′ , respectively. Now, let us investigate the spectrum of (4.52). We compute∆(T N N ′ (ϕ n,s is the eigenfunction of the Laplacian with charge Q obtained in (G.10). We can first show that See Appendix I for the derivation of (4.54). We then make an asymptotic expansion of (4.54) in the large-N limit aŝ This shows that the spectrum of (4.52) agrees with that in the continuum limit: In Appendix J, we show that the exact eigenvalue problem of (4.52) can be mapped to the so-called Hofstadter problem. Numerical studies of this problem also show that the eigenvalues are given as (4.56). Note that if we write U = X 1 + iX 2 and V = X 3 + iX 4 with the four Hermitian matrices X i corresponding to an embedding into R 4 , (4.48) is proportional to [X i , [X i , F ]], which is the natural Laplacian appearing in the matrix models. The matrix Laplacian (4.52) for rectangular matrices also naturally appears in the matrix models. For example, let us consider a block diagonal matrix configuration in the matrix models, are the configurations of the fuzzy torus with size N and N ′ , respectively. Then, the Laplacian [X i , [X i , F ]] in the matrix models reduces to (4.52) for the off-diagonal blocks of F , while it reduces to (4.48) for the diagonal blocks. Thus, (4.52) can be seen as the Laplacian of the open string modes connecting the two fuzzy tori. In fact, the same structure can also be found for the case of the sphere [31,32]. It is interesting that the natural matrix Laplacian [X i , [X i , F ]] reproduces in a unified way the both spectra of the charged and non-charged fields.

Discussion on fuzzy Riemann surfaces with higher genera
In this section, we discuss cases with higher genera.

Construction of Riemann surfaces with higher genera
Any Riemann surface with the genus greater than 1 can be constructed as the Poincaré disk with some identifications imposed. We first review this construction. See e.g. [33] for more details.
Let us consider a unit disk D 2 = {z ∈ C | |z| < 1} on the complex plane. We adopt the Poincaré metric on D 2 , which is the Kähler metric compatible with the standard complex structure on D 2 . The space (D 2 , g) is called the Poincaré disk. We consider a group SU(1, 1), which acts on D 2 as where a, b are complex numbers satisfying |a| 2 − |b| 2 = 1 and γ represents an element of SU(1, 1), For any γ ∈ SU(1, 1), the map z → γ(z) is an automorphism on D 2 preserving the Kähler structure. Note that γ and −γ give the same transformation on D 2 , so the automorphism group is isomorphic to P SU(1, 1) = SU(1, 1)/Z 2 . Let Γ be a Fuchsian group, which means a discrete subgroup of P SU(1, 1). Compact Riemann surfaces with genera greater than 1 are known to be constructed as a coset space M = D 2 /Γ. In this construction, all the information about the genus or the moduli of M is contained in Γ, and M is given by a set of all orbits on D 2 with respect to actions of Γ. By analogy with the torus, it is also useful to regard M as the disk D 2 with a nontrivial boundary condition imposed by actions of Γ. Note that this construction also gives a natural metric on M. Since Γ preserves g, the metric (5.1) also gives a local Kähler metric on M.

Geometric structures on M
We next consider charged scalars and spinor fields on M. For simplicity, we assume that the symplectic volume is V = M ω/2π = 1, where ω is defined by (A.2). From (5.1), a Kähler potential on M is given by Then, from (A.6) and (A.8), a spin connection and U(1) gauge field are In this gauge, we have Ω ab σ a σ b /4 = −iAσ 3 . For γ ∈ Γ given by (5.3), the gauge field A transforms as This is analogous to (4.8) on T 2 . Let ϕ (Q) be a complex scalar field coupling to A with charge Q and ψ = (ψ + , ψ − ) a spinor field on M. For γ ∈ Γ, they transform as ψ ± (z,z).

Dirac zero modes and automorphic forms
We next construct Dirac zero modes on M.
In our case, from (5.4) and (A.9), the Dirac operator is locally given by We find that the zero modes ψ (N ) = (ψ (N )+ , ψ (N )− ) take the following form: where h + and h − are holomorphic and anti-holomorphic functions. As we discussed in the previous subsection, their gauge transformations are given by ψ (N )± (γ(z),γ(z)) = bz + ā bz +ā  Thus, h ± are given by automorphic forms 8 . The automorphic forms h ± can be represented in terms of the Poincaré series [6,36] as Here, the summations are taken over all elements of Γ and f + and f − are arbitrary holomorphic and anti-holomorphic functions on D 2 , respectively. Note that for any γ ′ ∈ Γ, we have In order to construct (2.1) and (2.4), we need to construct an orthonormal basis of the zero modes by choosing h + (or f + ) in (5.14) appropriately. This should be done case-by-case, since it highly depends on the structure of Γ. For example, for the Bolza surface [34,35], which is the simplest example of surfaces with genus 2, the structure of Γ is relatively well-studied and it might be possible to obtain an orthonormal basis in this case. However, this is in general very difficult and is beyond the scope of this paper (see [36] for a formal construction) 9 .

Quantization on M
As explained in Appendix B, the quantization map can be defined as (B.2) in terms of the projection operators onto zero modes of the Dirac operator. This definition does not explicitly depends on the orthonormal basis of zero modes, which is only needed to write down the matrix elements of (B.2). Here, we show that (B.2) can be constructed for the Riemann surfaces considered above.
The projection is given by the Bergman kernel (B.3). From (5.14), one finds that the Bergman Kernel is given by [6]

Conclusion and discussion
In this paper, we considered the Berezin-Toeplitz quantization for local sections of nontrivial complex line bundles on Riemann surfaces. This corresponds to the matrix regularization of fields with nonvanishing U(1) charges in a nontrivial gauge flux. We argued that such fields are naturally mapped to rectangular matrices, while fields with vanishing charge are mapped to square matrices. We also showed that these mappings are embedded in the Berezin-Toeplitz quantization in a U(2) gauge theory in a unified way. We then explicitly constructed those mappings for the sphere and the torus and also discussed a possible extension to Riemann surfaces with higher genera. For the case of the sphere, we showed that this mapping reproduces the well-known fuzzy spherical harmonics. For the case of the torus, we found that the mapping produces rectangular matrices written in terms of elliptic functions and we also proposed a matrix Laplacian, which reproduces the spectrum on the commutative torus in the large-N limit.
In our examples of the sphere and the torus, the operation •, which is defined in (4.53) or (E.1), was used in constructing the matrix Laplacians. This operation gives another module structure of the matrix algebra satisfying where f and a are arbitrary smooth function and local section with charge Q, respectively, and W µν is a Poison tensor, given by the inverse of the symplectic form. { , } 1 only satisfies the first and the second properties of (6.1), while { , } 2 satisfies the first and the third, provided that the quantization condition [D µ , D ν ]a = −iQω µν a is satisfied. These operations cannot be the classical counterpart of •, because the violations of the classical counterpart of (6.1) contradict with the properties (1.1) and (1.3). Probably, we would need another structure to construct the classical operation in general. Nevertheless, we found that the continuum Laplacian on the sphere is proportional to {x i , {x i , } 2 } 2 , while that on the torus is proportional to {ū, {u, } 1 } 1 + {v, {v, } 1 } 1 , and we consider that these brackets may still have meanings when some special functions are put in the first slots. It should be important to understand this correspondence to find the description of the Laplacian in a general matrix geometry. The mapping for the local sections we considered in this paper should also be relevant for describing D-branes with gauge fluxes. In particular, it is known that the Berezin-Toeplitz quantization naturally appears in non-BPS D-brane systems with tachyon condensations [18,19] (see also [40,41]). In this context, introducing non-Abelian gauge groups seems to be quite natural, so that our formulation using the U(2) gauge theory will be naturally understood. This will be studied elsewhere. It is important to understand properties of the quantization with a general gauge group in order to understand its implications in the D-brane systems.
Our results will also shed light on the problem of describing curved spaces in the matrix models. For example, by using the mapping for rectangular matrices, it will be possible to generalize the work [31,32] and construct the large-N reduction for nontrivial U(1) bundles on Riemann surfaces, called Seifert manifolds.

A Dirac operator on Riemann surfaces with magnetic fluxes
In this appendix, we construct a Dirac operator on a general Riemann surface M with a magnetic flux.
Let z be a local complex coordinate on an open subset U ⊂ M. We define the standard complex structure J on M by J(∂ z ) = i∂ z and J(∂z) = −i∂z. Note that this definition does not depend on the choice of the local coordinate. Let g be a Kähler metric on M compatible with J.
On U, we have g zz = gzz = 0 and we can write g as g = 2g zz dzdz, (A.1) where g zz = g(∂ z , ∂z). We define a symplectic form on M by ω( · , · ) = g( J· , · ). In terms of the local coordinate, we can write ω as We choose a U(1) gauge field A which satisfies ω = V dA, as explained in section 2.1. Let e a be the zweibein for the Kähler metric (A.1). They are explicitly given by Note that from the positivity of the metric, g zz is always positive. The inverse θ a of e a is given by The spin connection is determined by and Ω ab = −Ω ba . By solving this equation, we find Any Kähler metric is locally written as in terms of the Kähler potential ρ, which is a real function defined locally on U. The geometric structures we introduced above can also be expressed in terms of ρ. For example, from (A.2) and ω = V dA, the gauge field A is given by up to the gauge transformation. For M, the Dirac operator D defined by (2.2) flips the chirality, so that it has only off-diagonal elements. Using the above data, we can express D as where D ± are the matrix elements of D acting on the spaces with chirality ±1, respectively.

B Bergman kernel
In this section, we give a formulation of the maps (2.1) and (2.4) in terms of the Bergman kernel. Let us consider a product ϕ (Q) ψ (N ) I of a charged scalar fields ϕ (Q) with charge Q and a Dirac zero mode ψ (N ) I with charge N, where I = 1, 2, · · · , N. This product has the total charge N ′ = N + Q, and can be expanded in terms of the Dirac eigen modes with charge N ′ as where · · · stands for the terms of non-zero modes. If ψ  [15]: The projections are given by the so-called Bergman kernel, Here, the spinor indices are not contracted, so that K (N ) (z, w) is a 2 × 2 matrix with those indices. The projection Π (N ) is then defined by Note that the original expression (2.4) is just the matrix representation of (B.2) 10 : In this appendix, we give a proof of (1.3) following [15]. We first show that Dirac eigen modes with charge N have a large energy gap in the large-N limit. Let D be a Dirac operator for the charged spinors. Let us consider the action of D on a two-component spinor χ. Since D generally flips the chirality, D can be represented as where χ + and χ − are the positive and negative chirality modes of χ, respectively. Below, we assume that D is normalized such that it is Hermitian and (D ± ) † = D ∓ . If χ is a normalized eigen mode of D with eigenvalue E, we have Then, we can estimate the energy E as Thus,  N 0 ). Also, the norm a is finite in the large-N limit. The only nontrivial factor is [Π (N ) , f ] and we will estimate this in the following. Let us consider an operator (1 + αD 2 ) −1 , where α is an Nindependent positive number. Note that 1 + αE 2 > 0 for any eigenvalue E, so that the operator (1 + αD 2 ) −1 is well-defined. Since E ≥ O(N 1/2 ) except for the case E = 0 as shown in (C.4), the operator (1 + αD 2 ) −1 behaves as a projection onto KerD for sufficiently large values of N. Hence, we obtain Thus, the problem of estimating [Π (N ) , f ] reduces to that of [(1 + αD 2 ) −1 , f ]. In order to evaluate the latter, we first rewrite this as (C.10) Thus, we obtain This gives the following estimation: where E 1 is the smallest nonzero eigenvalue of D. The last inequality is obtained as follows. For any eigenvalue E, we have the relation, This implies that which, together with the obvious relation (1 + αD 2 ) −1 ≤ 1, leads to the second inequality in (C.12). By applying (C.12) to (C.6), we finally obtain Since |E 1 | ≥ O(N 1/2 ), the right-hand side vanishes in the large-N limit and we find that (1.3) is indeed satisfied.

D Monopole harmonics
In this appendix, we review the definition of the monopole harmonics. See [24,25] for more details. We first introduce linear operators which is locally defined on U 1 and U 2 as where the upper and lower signs represent the expressions on U 1 and U 2 , respectively. These operators are the angular momentum operator in the presence of a monopole with magnetic charge Q/2 at the origin and reduces to the ordinary angular momentum operators when Q = 0. In fact, L and the orthonormal condition for a fixed l, where ω is the volume form for the metric (3.1). The concrete expression of Y A are the (2l + 1)-dimensional irreducible representation of the generators of SU(2) and |lm are the standard basis of the representation space. Again, the upper and lower signs represent the expressions defined on U 1 and U 2 , respectively.
The following formula is very useful.
where C l 1 m 1 l 2 m 2 l 3 m 3 is the Clebsch-Gordan coefficient. For the gauge invariance of the left hand side, Q 1 = Q 2 + Q 3 must hold.
We first define linear operators on M N ×N ′ (C) by and the orthonormal condition for a fixed l. Here, the trace is defined over N × N ′ matrices. In terms of the basis { |Jr J ′ r ′ |}, they are expressed asŶ F Detailed calculation of the normalization factor N r In this appendix, we derive (4.26). From of the explicit form of the Jacobi-theta function, we first write Then, by shifting the integration variable as z → z − ζ and substituting z = x + τ y, we obtain dy × e −2N π Im(τ )y 2 e i2π{(r+N l)τ −(r ′ +N l ′ )τ }y e i2π{r−r ′ +N (l−l ′ )}x .

(F.2)
The integration over x just produces the Kronecker delta factor δ r,r ′ δ l,l ′ . Thus, by taking the summation over l ′ and we obtain This can also be written in a compact form as By again shifting the integration variable as y → y − l − r/N, we obtain Since the l-dependence appears only in the integration range, summing up all l ∈ Z is equivalent to extending the integration range to (−∞, ∞). Thus, we finally arrive at a simple Gaussian integral. The final result is

G Orthonormal basis of local sections on the torus
In this appendix, we construct an orthonormal basis of local sections of the nontrivial line bundle with charge Q.
As the orthonormal basis, we consider a set of eigenfunctions of the Laplacian. Let us consider the Laplacian for charge Q ∈ N given by In the following, we will solve the eigenvalue problem, to find the orthonormal eigen modes ϕ (Q) n as well as the eigenvalues E n . Here, the eigen modes shall be ordered as E n < E n+1 (∀n ∈ Z ≥0 ).
We first introduce the creation-annihilation operators aŝ which satisfy the commutation relation Then, (G.3) can be expressed in terms of the number operatorN :=â †â as

H Useful relations for the eigenstates of the Laplacian
In this appendix, we show some useful identities for the eigenstates of the Laplacian. Let us consider the product of the eigenstates (G.10), whereâ † y andâ † z stand for the creation operators (G.4) acting on the complex variables y and z, respectively 11 . By using the identity of the Jacobi-theta function, where we used N ′ = N + Q and X and Y are defined by (H.4) If we regard (H.4) as a change of variables from (y, z) to (X, Y ), we can also convertâ † y andâ † z to those in the X-and Y -coordinates:â † (H.5) 11 In this appendix, we use y as a complex variable exceptionally, while it is used as a real variable in the other sections.

J Laplacian for rectangular matrices and Hofstadter problem
In this appendix, we consider the exact eigenvalue problem of the matrix Laplacian (4.52) for rectangular matrices. We show that the problem is equivalent to a special case of the Hofstadter problem, which we will review below. The eigenvalue equation for N × N ′ matrices is written aŝ In terms of the matrix elements of F , this is equivalent to F r+1,r ′ +1 + 2 cos 2π whereẼ is given byẼ The periodic structure of (J.2) enables us to extend the range of indices as F r+N,r ′ +N ′ = F r,r ′ . With this notation, assuming that N and N ′ are coprime, we relabel the matrix elements as F r := F r,r (J.4) for r = 0, 1, · · · , NN ′ − 1. In this notation, (J.2) reduces to F r+1 + 2 cos 2Qπr NN ′ F r + F r−1 =ẼF r (J.5) for r = 0, 1, · · · , NN ′ − 1, where F −1 := F N N ′ −1 and F N N ′ := F 0 . This is also equivalent to the following eigenvalue problem: The eigenvalue problem of H is what is known as the Hofstadter problem [42]. Finding an exact solution to this problem is still an open problem, though some numerical analyses have been done [42] and revealed a fractal structure of the spectrum, known as a Hofstadter butterfly. It is interesting that the same Hofstadter problem also arises in a system of tight-binding Bloch electrons under a constant magnetic flux in a periodic two-dimenisonal surface, which has the following Hamiltonian: where t is the hopping parameter, q is the number of lattice sites, d i,σ ( k) is a creation and annihilation operator for the wave number k and spin σ, satisfying anti-commutation relations {d i,σ ( k) , d † j,σ ′ ( k ′ )} = δ ij δ σσ ′ δ k k ′ . The eigenvalue ω i ( k) is obtained by solving the eigenvalue problem of where φ = p q is the the strength of U(1) flux per unit plaquette and p is the Chern number. Here, we assumed that p and q are coprime for simplicity. Readers may refer to [43] for the derivations of the Hamiltonian (J.9) and the matrix (J.10). If we put φ = Q N N ′ , q = NN ′ and k = 0, the matrix (J.10) reduces to the Hamiltonian (J.7) for the matrix Laplacian.
The spectrum of (J.7) or (J.10) has been studied numerically. In the large-q limit, it is shown that the spectrum of (J.7) indeed coincides with the Landau level [44], and this is consistent with our result.

K Evaluating the norm of ψ −
In this appendix, we show that the norm of ψ − , which is given by (5.9) and (5.12), does not converge.
First, for the inner product (2.3), the norm of ψ − can be rewritten as − (w)f − (η(w)). (K.1) To obtain the second equality we changed the dummy variable from γ ′ to η = γ ′ γ −1 and to obtain the third equality, we changed the integral variable by w = γ(z). Note that ω is invariant under actions of Γ. To obtain the last equality, we used the fact that for any γ ∈ Γ, the relation This shows that for N ≥ 1, the integration of |w| 2 does not converge. Thus, ψ − 2 is not convergent for N ≥ 1.

L Bergman kernel on disk
In this appendix, we construct a Bergman kernel on the Poincaré disk D 2 . See [5] for more details.
On the Poincaré disk, an orthonormal basis of the Dirac zero mode is given by ψ n (z,z) = (1 − |z| 2 ) (N +1)/2 N 2π 1/2 N + n N 1/2 z n . (L.1) Here, n = 1, 2, . . . , ∞, so the dimension of KerD is infinity. This comes from the noncompactness of D 2 . We can check the orthonormality as follows. x n , we find that the Bergman kernel is given by