text:
# 1 Introduction

## 1 Introduction

A recent excitement in string theory is that we finally arrive at concrete proposals for nonperturbative definitions of superstring theory . In particular, Matrix Theory  and the IIB matrix model , which are candidates for a nonperturbative definition of M-theory and type IIB superstring theory, respectively, have attracted considerable interest (for reviews, see Refs. ). The proposed formulations take the form of large $`N`$ reduced models , which can be obtained by dimensional reduction of large $`N`$ gauge theories; a reduction to one dimension for M-theory, and to zero dimension (one point) for type IIB superstring theory. These proposals are supported by some evidences such as the similarity of the Hamiltonian (or the action) to that of membranes or strings, the appearance of soliton-type objects known as D-branes with consistent interactions, and the consistency with string dualities upon compactification. For the IIB matrix model, even an attempt to establish a direct connection to perturbative string theory has been made by deriving the light-cone string field Hamiltonian from loop equations of the model . This attempt was indeed successful, albeit with the aid of symmetry and power-counting arguments. Further variants of that model have been proposed in Refs. .

As another approach to analyze these proposals we can investigate the dynamical properties of large $`N`$ reduced models of this kind, and verify if they really have the potential to describe nonperturbative string theory. In Ref.  a two-dimensional reduced model with unitary matrices has been studied in this context. There, a large $`N`$ limit, which differs from the planar limit (or ’t Hooft limit), has been discovered numerically 111This non-planar large $`N`$ limit has recently been re-interpreted as a continuum limit of non-commutative gauge theory .. A Hermitian matrix model obtained by simply omitting the fermions in the IIB matrix model, and its generalizations to arbitrary dimensions larger than two, have been studied in Refs. . In Ref. , Monte Carlo simulations up to $`N=256`$ have been reported and analytical methods such as perturbation theory, Schwinger-Dyson equations and $`1/D`$ expansions have been applied, providing a comprehensive understanding of the large $`N`$ dynamics of that model. A new type of Monte Carlo technique was used to extract the value of the partition function . This technique has been further applied to extract the asymptotic behavior of the eigenvalue distribution for large eigenvalues .

In the present paper, we make a first attempt to extract the large $`N`$ dynamics of a supersymmetric large $`N`$ reduced model obtained by dimensional reduction of 4D supersymmetric Yang-Mills theory. The model can be regarded as a 4D counterpart of the IIB matrix model. The bosonic model is well understood , but the inclusion of fermions makes the system far more complicated. An attempt to study the model analytically maps it onto a soluble system , but the relevance to the original model is unclear due to a nontrivial change of variables in an analytic continuation. Here we take a direct approach, based on Monte Carlo simulations. Fermions are completely included by the use of the so-called Hybrid-R algorithm , which is one of the standard methods in QCD simulations with dynamical quarks.

One of the features that makes the IIB matrix model most attractive as a nonperturbative definition of string theory is that space-time is dynamically generated as the eigenvalue distribution of the bosonic matrices . In Ref.  a low energy effective theory of the model is constructed, where the authors discuss some possible mechanisms that may induce a collapse of the eigenvalue distribution to a four-dimensional manifold. We extract the large $`N`$ behavior of the space-time extent in our model and compare the result with the prediction obtained by the low energy effective theory. Another dynamical issue to be addressed in this context is the space-time uncertainty relation, which was proposed as a principle for constructing nonperturbative string theory . We extract the large $`N`$ behavior of the space-time uncertainty of our model and confirm that the model indeed satisfies the proposed principle.

Another attractive feature of the IIB matrix model as a nonperturbative definition of string theory is that its only parameter $`g`$ is a simple scale parameter 222This means in particular that the string coupling constant, which is related to the vacuum expectation value of the dilaton field, is not a tunable parameter. We come back to this point in Section 4.3.. One has to tune $`g`$ suitably as one sends $`N`$ to infinity, so that the correlation functions have finite large $`N`$ limits. According to Ref. , Wilson loop operators can be interpreted as the string creation and annihilation operators, and it was found that $`g^2N`$ should be fixed in order to obtain the light-cone string field Hamiltonian in the large $`N`$ limit. It is a non-trivial test of the model to verify if the correlation functions of Wilson loops really have a universal large $`N`$ scaling. We address this issue in the present model and show that there is indeed a universal large $`N`$ scaling at fixed $`g^2N`$.

We also address yet another important dynamical issue in this model, namely the question of equivalence to ordinary super Yang-Mills theory in the sense of Eguchi and Kawai , which is exactly the way large $`N`$ reduced models first appeared in history. The crucial observation is that large $`N`$ gauge theory does not depend on the volume (under some assumptions), which inspired Eguchi and Kawai to propose the zero-volume limit of large $`N`$ gauge theory as a model equivalent to the gauge theories in an infinite volume . One of the assumptions is that the $`(\text{}_N)^D`$ symmetry of the model is not spontaneously broken, where $`D`$ is the space-time dimension. However, in the purely bosonic case in $`D>2`$, the symmetry is spontaneously broken at weak coupling , thus preventing one from taking a continuum limit. This led to modifications of the model  so that the $`(\text{}_N)^D`$ symmetry is not spontaneously broken while keeping the equivalence valid. In the supersymmetric case, the effective action which induces the spontaneous symmetry breaking of the (N)D symmetry is naively cancelled by the contributions of fermions. Indeed, in the scalar field case, it has been shown that the reduced model is equivalent to the field theory without such modifications . We observe in the present supersymmetric model that the Eguchi-Kawai equivalence indeed holds at least in a finite range of scale. What is rather remarkable is that actually this is true also for the bosonic case, which is contrary to what has been expected.

In Section 2 we describe the model we are going to investigate. In Section 3 we study the space-time structure of the model. In Section 4 we present our results for correlation functions of Wilson loop and Polyakov line operators, and we discuss the Eguchi-Kawai equivalence as well as the universal scaling behavior. Section 5 is devoted to a summary and discussion. In Appendix A we comment on the algorithm we used for the simulation. In Appendix B we present the corresponding results for the bosonic case for comparison.

## 2 The model

The model we investigate is a supersymmetric matrix model obtained by dimensional reduction of 4D SU($`N`$) super Yang-Mills theory. The partition function is given by

$`Z`$ $`=`$ $`{\displaystyle \text{d}A\text{e}^{S_b}\text{d}\psi \text{d}\overline{\psi }\text{e}^{S_f}},`$
$`S_b`$ $`=`$ $`{\displaystyle \frac{1}{4g^2}}\text{tr}[A_\mu ,A_\nu ]^2,`$
$`S_f`$ $`=`$ $`{\displaystyle \frac{1}{g^2}}\text{tr}\left(\overline{\psi }_\alpha (\mathrm{\Gamma }^\mu )_{\alpha \beta }[A_\mu ,\psi _\beta ]\right),`$ (2.1)
where $`A_\mu `$ ($`\mu =1,\mathrm{},4`$) are traceless $`N\times N`$ Hermitian matrices, and $`\psi _\alpha `$, $`\overline{\psi }_\alpha `$ ($`\alpha =1,2`$) are traceless $`N\times N`$ complex matrices. The measure is defined as

$`\text{d}\psi \text{d}\overline{\psi }`$ $`=`$ $`{\displaystyle \underset{\alpha =1}{\overset{2}{}}}\left[{\displaystyle \underset{i,j=1}{\overset{N}{}}}\left[\text{d}(\psi _\alpha )_{ij}\text{d}(\overline{\psi }_\alpha )_{ij}\right]\delta \left({\displaystyle \underset{i=1}{\overset{N}{}}}(\psi _\alpha )_{ii}\right)\delta \left({\displaystyle \underset{i=1}{\overset{N}{}}}(\overline{\psi }_\alpha )_{ii}\right)\right],`$ (2.2)
$`\text{d}A`$ $`=`$ $`{\displaystyle \underset{\mu =1}{\overset{4}{}}}\left[{\displaystyle \underset{i<j}{}}\{\text{d}\text{Re}(A_\mu )_{ij}\text{d}\text{Im}(A_\mu )_{ij}\}{\displaystyle \underset{i=1}{\overset{N}{}}}\{\text{d}(A_\mu )_{ii}\}\delta \left({\displaystyle \underset{i=1}{\overset{N}{}}}(A_\mu )_{ii}\right)\right].`$ (2.3)
This model is invariant under 4D Lorentz transformations 333When one defines the IIB matrix model nonperturbatively, a Wick rotation to Euclidean signature in needed. This is also the case for the present model. Hence by Lorentz invariance we actually mean rotational invariance., where $`A_\mu `$ transforms as a vector and $`\psi _\alpha `$ as a Weyl spinor. $`\mathrm{\Gamma }_\mu `$ are 2 $`\times `$ 2 matrices acting on the spinor indices, and they can be given explicitly as

$$\mathrm{\Gamma }_1=i\sigma _1,\mathrm{\Gamma }_2=i\sigma _2,\mathrm{\Gamma }_3=i\sigma _3,\mathrm{\Gamma }_4=\mathrm{𝟏}.$$
(2.4)
The model is manifestly supersymmetric, and it also has a SU($`N`$) symmetry

$$A_\mu VA_\mu V^{};\psi _\alpha V\psi _\alpha V^{};\overline{\psi }_\alpha V\overline{\psi }_\alpha V^{},$$
(2.5)
where $`V\text{SU}(N)`$. All these symmetries are inherited from the super Yang-Mills theory before dimensional reduction. The model can be regarded as the four-dimensional counterpart of the IIB matrix model .

In contrast to unitary matrix models, where the integration domain for the partition function is compact, the first nontrivial question to be addressed in Hermitian matrix models in general, is whether the model is well-defined as it stands. The problem can be most clearly understood by decomposing the Hermitian matrices into eigenvalues and angular variables, where a potential danger of divergence exists in the integration over the eigenvalues, even at finite $`N`$. This issue has been addressed numerically for the supersymmetric case  at $`N=3`$ as well as the bosonic case  up to $`N=6`$. Exact results are available for $`N=2`$ . There is also a perturbative argument which is valid when all the eigenvalues are well separated from each other . This reasoning agrees with the conclusions obtained for small $`N`$ . In particular, supersymmetric models in $`D=4,6,10`$ are expected to be well-defined for arbitrary $`N`$. Our simulations confirm that this is indeed the case for $`D=4`$.

Since the model is well-defined without any cutoff, the parameter $`g`$, which is the only parameter of the model, can be absorbed by rescaling the variables,

$`A_\mu `$ $`=`$ $`g^{1/2}X_\mu ,`$ (2.6)
$`\psi _\alpha `$ $`=`$ $`g^{3/4}\mathrm{\Psi }_\alpha .`$ (2.7)
Therefore, $`g`$ is a scale parameter rather than a coupling constant, i.e. the $`g`$ dependence of physical quantities is completely determined on dimensional grounds. The parameter $`g`$ should be tuned appropriately as one sends $`N`$ to infinity, so that each correlation function of Wilson loops has a finite large $`N`$ limit. Whether such a limit really exists or not is one of the dynamical issues we address in this work.

We now discuss the Eguchi-Kawai equivalence , which is the equivalence between reduced models and the corresponding gauge theories in the large $`N`$ limit. In its proof based on the Schwinger-Dyson equation, one has to assume quantities of the type

$$\text{tr}(\text{e}^{ik_\mu A_\mu })\text{tr}(\text{e}^{ik_\mu A_\mu })(k_\mu 0)$$
(2.8)
to vanish. Assuming in addition large $`N`$ factorization, the vanishing of (2.8) is equivalent to $`\text{tr}(\text{e}^{ik_\mu A_\mu })`$$`=0`$, which is guaranteed if the eigenvalues of $`A_\mu `$ are uniformly distributed on the whole real axis in the large $`N`$ limit. The fact that the present model is well-defined without any cutoff implies that the eigenvalue distribution of $`A_\mu `$ is not uniform, but it has a finite extent for finite $`N`$. Hence, the Eguchi-Kawai equivalence is quite nontrivial even in the supersymmetric case. Here the situation is more subtle than in the case of the unitary matrix model version of a large $`N`$ reduced model . There, the model has the $`(\text{}_N)^D`$ symmetry $`U_\mu \text{e}^{2\pi im_\mu /N}U_\mu `$ ($`m_\mu =0,1,\mathrm{},N1`$), hence quantities like $`\text{tr}(U_\mu )^n`$ vanish, unless the symmetry is spontaneously broken.

One might be tempted to consider a model defined by the partition function (2.1) but without imposing the traceless condition on $`A_\mu `$. We denote such a model as the U($`N`$) model, to be distinguished from the original model, which we call the SU($`N`$) model. The U($`N`$) model has the U(1)4 symmetry

$$A_\mu A_\mu +\alpha _\mu \mathrm{𝟏}_N,$$
(2.9)
where $`\alpha _\mu `$ is a real vector. Note, however, that the trace part of $`A_\mu `$ in the U($`N`$) model simply decouples because $`A_\mu `$ appears in the action only through commutators. The transformation (2.9) acts on the decoupled trace part and hence it cannot play any physical rôle. Indeed the quantity (2.8) calculated with the U($`N`$) model or with the SU($`N`$) model is exactly the same. Thus, considering the U($`N`$) model does not help.

Next we comment on the method we use to study the model. Details can be found in Appendix A. The integration over fermionic variables can be done explicitly and the result is given by $`det`$, $``$ being a $`2(N^21)`$ $`\times `$ $`2(N^21)`$ complex matrix which depends on $`A_\mu `$. Hence the system we want to simulate can be written in terms of bosonic variables as

$$Z=\text{d}A\text{e}^{S_b}det.$$
(2.10)
A crucial point for the present work is that the determinant $`det`$ is actually real positive, as we prove in Appendix A. Due to this property, we can introduce a $`2(N^21)`$ $`\times `$ $`2(N^21)`$ Hermitian positive matrix $`𝒟=^{}`$, so that $`det=\sqrt{det𝒟}`$, and the effective action of the system takes the form

$$S_{\text{eff}}=S_b\frac{1}{2}\mathrm{ln}det𝒟.$$
(2.11)
We apply the Hybrid R algorithm  to simulate this system. In the framework of this algorithm, each update of a configuration is made by solving a Hamiltonian equation for a fixed “time” $`\tau `$. The algorithm is plagued by a systematic error due to the discretization of $`\tau `$ that we used to solve the equation numerically. We performed simulations at three different values of the time step $`\mathrm{\Delta }\tau `$. Except in Fig. 2, we find that the results do not depend much on $`\mathrm{\Delta }\tau `$ (below a certain threshold), so we just present the results for the value $`\mathrm{\Delta }\tau =0.002`$, which appears to be sufficiently small.

We also note that there is an exact result

$$\text{tr}F^2=\text{tr}(\underset{\mu \nu }{}[A_\mu ,A_\nu ]^2)=6g^2(N^21),$$
(2.12)
which can be obtained by a scaling argument, similar to the one used for the bosonic case . We used this exact result to check the code and the numerical accuracy.

## 3 The space-time structure

We first study the space-time structure of the reduced model. In the IIB matrix model, the eigenvalues of the bosonic matrices $`A_\mu `$ are interpreted as the space-time coordinates . However, since the matrices $`A_\mu `$ are not simultaneously diagonizable in general, the space-time is not classical. In order to extract the space-time structure, we first define the space-time uncertainty $`\mathrm{\Delta }`$ by

$$\mathrm{\Delta }^2=\frac{1}{N}\text{tr}(A_\mu ^2)\underset{U\text{SU}(N)}{\mathrm{max}}\frac{1}{N}\underset{i}{}\{(UA_\mu U^{})_{ii}\}^2,$$
(3.1)
which is invariant under Lorentz transformation and SU($`N`$) transformation (2.5) . This formula has been derived in Ref.  based on analogy to quantum mechanics, regarding $`A_\mu `$ as an operator acting on a space of states. It has the natural property that $`\mathrm{\Delta }^2`$ vanishes if and only if the matrices $`A_\mu `$ are diagonalizable simultaneously. For each configuration $`A_\mu `$ generated by a Monte Carlo simulation, we maximize $`_i\{(UA_\mu U^{})_{ii}\}^2`$ with respect to the SU($`N`$) matrix $`U`$. We denote the matrix which yields the maximum as $`U_{\mathrm{max}}`$, and we define $`x_{\mu i}=(U_{\mathrm{max}}A_\mu U_{\mathrm{max}}^{})_{ii}`$ as the space-time coordinates of $`N`$ points ($`i=1,\mathrm{},N`$) in four-dimensional space-time.

Note that $`x_{\mu i}`$ should be identified with the dynamical variables denoted by the same $`x_{\mu i}`$ in Ref. . There, the bosonic matrices $`A_\mu `$ and the fermionic matrices $`\psi _\alpha `$ are decomposed into diagonal and off-diagonal elements as

$`(A_\mu )_{ij}`$ $`=`$ $`x_{\mu i}\delta _{ij}+a_{\mu ij}(a_{\mu ii}=0),`$
$`(\psi _\alpha )_{ij}`$ $`=`$ $`\xi _{\alpha i}\delta _{ij}+\phi _{\alpha ij}(\phi _{\alpha ii}=0).`$ (3.2)
The off-diagonal parts $`a_{\mu ij}`$ and $`\phi _{\alpha ij}`$ are integrated out using the “Lorentz gauge” in the one-loop approximation, which is valid when the points $`x_{\mu i}`$ $`(i=1,\mathrm{},N)`$ are well separated from each other. Thus one obtains the effective action for $`x_{\mu i}`$ and $`\xi _{\alpha i}`$, which can be considered as a low-energy effective action of the supersymmetric large $`N`$ reduced model. In order to get the effective action only for $`x_{\mu i}`$, one still has to integrate over $`\xi _{\alpha i}`$, which cannot be done exactly for $`D=6`$ and $`D=10`$ (IIB matrix model). In $`D=4`$, however, the integration over $`\xi _{\alpha i}`$ can be carried out exactly and the system of $`x_{\mu i}`$ is described by a simple branched polymer with an attractive potential between the points connected by a bond. In $`D=6`$ and $`D=10`$, the system of $`x_{\mu i}`$ is expected to be described by some complicated branched-polymer like structure. Thus, although the one-loop approximation might seem quite drastic, the low energy effective theory of $`x_{\mu i}`$ still has a nontrivial dynamics. In Ref. , some plausible mechanisms for the collapse of the $`x_{\mu i}`$ distribution in IIB matrix model have been discussed. What we have described in the previous paragraph provides a way to extract the low-energy effective theory of $`x_{\mu i}`$ from the full model without perturbative expansions. In particular, we can check explicitly whether the one-loop approximation adopted in Ref.  really captures the low energy dynamics of the supersymmetric large $`N`$ reduced model.

We first look at the distribution $`\rho (r)`$ of the distances $`r`$, where the distance between two arbitrary points $`x_ix_j`$ is measured by $`\sqrt{(x_ix_j)^2}`$. In Fig. 1 we plot the results for $`N=16,\mathrm{\hspace{0.17em}24},\mathrm{\hspace{0.17em}32}`$ and 48. We first note that the distribution at small $`r`$ falls off rapidly below $`r/\sqrt{g}1.5`$, independently of $`N`$. (This behavior is also seen in the bosonic case shown in Fig. 10.) This observation is in agreement with the argument in Ref.  that the ultraviolet behavior of the space-time structure of the model is controlled by the SU(2) matrix model. There, this argument has been used to justify the introduction of a $`N`$-independent ultraviolet cutoff in the low energy effective theory, which otherwise suffers from ultraviolet divergence due to coinciding $`x_{\mu i}`$’s. Our observation confirms that the ultraviolet cutoff is indeed generated dynamically if one treats the full model nonperturbatively instead of making perturbative expansions around diagonal matrices.

In both, the supersymmetric as well as the bosonic case, we observe that the distribution shifts towards larger $`r`$ as one increases $`N`$. In order to quantify this behavior, we define the extent of space-time by

$$R_{\mathrm{new}}=\frac{2}{N(N1)}\underset{i<j}{}\sqrt{(x_ix_j)^2}=_0^{\mathrm{}}\text{d}rr\rho (r).$$
(3.3)
We denote this quantity by $`R_{\mathrm{new}}`$ in order to distinguish it from the definition of the extent of the space-time $`R=\sqrt{\frac{1}{N}\text{tr}(A_\mu ^2)}`$ used in Ref. . $`R`$, which roughly corresponds to $`\sqrt{_0^{\mathrm{}}\text{d}rr^2\rho (r)}`$, is logarithmically divergent in the 4D supersymmetric case due to the asymptotic behavior $`\rho (r)r^3`$ at large $`r`$ . On the other hand, $`R_{\mathrm{new}}`$ does not suffer from this divergence as eq. (3.3) shows. In Fig. 2 we plot the results for the space-time extent $`R_{\mathrm{new}}`$ as well as those for the space-time uncertainty $`\sqrt{\mathrm{\Delta }^2}`$ for $`N=16,24,32`$ and 48. We repeat the same measurements for the bosonic model with $`N`$ up to 256 and include the results in Fig. 2 for comparison. We see that the effect of fermions enhances $`R_{\mathrm{new}}`$ and suppresses $`\sqrt{\mathrm{\Delta }^2}`$ considerably. However, the power of the large $`N`$ behavior does not seem to be affected.

Let us discuss the results for $`R_{\mathrm{new}}`$. In the bosonic case, the data can be nicely fitted to a power behavior with $`R_{\mathrm{new}}/\sqrt{g}=1.56(1)N^{1/4}`$. As expected, the observed large $`N`$ behavior of $`R_{\mathrm{new}}`$ is the same as the one obtained for $`R`$ in Ref. . In the supersymmetric case, the large $`N`$ behavior of $`R_{\mathrm{new}}`$ can be predicted by the branched polymer picture based on the one-loop approximation . Since the Hausdorff dimension of branched polymers is four, $`d_\text{H}=4`$, the number of points $`N`$ grows as the extent $`R_{\mathrm{new}}`$ of the branched polymer, $`N(R_{\mathrm{new}}/\mathrm{})^{d_\text{H}}`$. Here $`\mathrm{}`$ is the minimum length of the bond, which is of $`O(\sqrt{g})`$ as we have already discussed. Thus one obtains $`R_{\mathrm{new}}\sqrt{g}N^{1/4}`$. The data in Fig. 2 seem to be consistent with this prediction. Fitting the data to this power behavior, we obtain $`R_{\mathrm{new}}/\sqrt{g}=3.30(1)N^{1/4}`$.

One might be surprised that supersymmetry does not affect the power of the large $`N`$ behavior of the space-time extent $`R_{\mathrm{new}}`$. We recall, however, that in the bosonic case the explanation is completely different — although the power is the same . There the one-loop perturbative expansion around diagonal matrices yields a logarithmic attractive potential between all the pairs of eigenvalues. The one-loop effective potential is dominant as far as the extent of the eigenvalue distribution is larger than $`\sqrt{g}N^{1/4}`$. One can therefore put an upper bound on the space-time extent $`R\sqrt{g}N^{1/4}`$. What happens actually is that this upper bound is saturated. The behavior $`R\sqrt{g}N^{1/4}`$ can also be shown to all orders in the $`1/D`$ expansion .

Let us turn to the results for the space-time uncertainty. Using the one-loop perturbative expansion, $`\mathrm{\Delta }^2`$ can be roughly estimated as

$$\mathrm{\Delta }^2=\frac{1}{N}\underset{ij}{}a_{\mu ij}a_{\mu ji}=\frac{1}{N}\underset{ij}{}\frac{g^2}{(x_ix_j)^2}\frac{g^2N}{R_{\mathrm{new}}^2}.$$
(3.4)
The powers of $`R_{\mathrm{new}}`$ and $`\sqrt{\mathrm{\Delta }^2}`$, as well as the coefficients we observe, are in qualitative agreement with this estimation. The bosonic case has been studied before in Ref. . The data in Fig. 2 can be nicely fitted to a power behavior with $`\sqrt{\mathrm{\Delta }^2/g}=0.907(3)N^{1/4}`$. Thus, in the bosonic case we obtain $`\sqrt{\mathrm{\Delta }^2}0.58R_{\mathrm{new}}`$ , which indicates a signficant deviation from the classical space-time picture. On the other hand, in the supersymmetric case we obtain $`\sqrt{\mathrm{\Delta }^2/g}=0.730(2)N^{1/4}`$, hence our result amounts to $`\sqrt{\mathrm{\Delta }^2}0.22R_{\mathrm{new}}`$, coming closer to the classical space-time picture.

We will see in the next section that the scale parameter $`g`$ should be taken to be $`O(1/\sqrt{N})`$ in order to obtain a universal scaling behavior for the Wilson loop correlators. This means that the space-time uncertainty in the physical scale remains finite, rather than vanishing, in the large $`N`$ limit. Therefore the present model satisfies the space-time uncertainty principle proposed for nonperturbative definitions of string theories .

## 4 Wilson loop correlation functions

In the interpretation of a large $`N`$ reduced model as a string theory, Wilson loop operators correspond to string creation operators . Therefore, the existence of a non-trivial large $`N`$ limit of the Wilson loop correlators is an absolutely crucial issue. It has been addressed before in the 2D Eguchi-Kawai model, where non-trivial large $`N`$ scaling has indeed been observed .

We define the “Wilson loop” and the “Polyakov line” operators as

$`W(k)`$ $`=`$ $`{\displaystyle \frac{1}{N}}\text{tr}(\text{e}^{ikX_1}\text{e}^{ikX_2}\text{e}^{ikX_1}\text{e}^{ikX_2}),`$
$`P(k)`$ $`=`$ $`{\displaystyle \frac{1}{N}}\text{tr}(\text{e}^{ikX_1}),`$ (4.1)
where $`X_\mu `$ are dimensionless matrices defined in eq. (2.6). For convenience we have chosen particular components of $`X_\mu `$ in the above definitions, but the choice of the directions becomes irrelevant when taking the vacuum expectation value, due to Lorentz symmetry and parity invariance. In the actual calculations we take an average over all possible choices of the components in order to enhance the statistics.

The real parameter $`k`$ represents the dimensionless “momentum” that characterizes the momentum density distributed along the string. The physical (dimensionful) momentum variable is given by $`k_{\text{phys}}=k/\sqrt{g}`$. We have to tune $`g`$ depending on $`N`$, so that the correlation functions of the above operators have definite large $`N`$ limits as functions of $`k_{\text{phys}}`$. In the following, we always set $`g=1`$ for $`N=48`$ without loss of generality.

In all plots except for Fig. 4, we further assume $`g`$ to be proportional to $`1/\sqrt{N}`$. This turns out to be consistent with large $`N`$ scaling, hence $`g1/\sqrt{N}`$ can be regarded as one of our observations.

### 4.1 One-point function and Eguchi-Kawai equivalence

In this subsection we discuss the one-point functions, and we start with the Wilson loop $`W(k)`$. Also Ref.  presents some recent results about this quantity.

In the small $`k`$ regime it can be expanded as

$`W(k)`$ $`=`$ $`1+{\displaystyle \frac{1}{2N}}k^4\text{tr}([X_1,X_2]^2)+O(k^6)`$ (4.2)
$`=`$ $`1{\displaystyle \frac{1}{4}}k^4\left(N{\displaystyle \frac{1}{N}}\right)+O(k^6),`$
where we have used the exact result (2.12). Therefore, in order to make the small $`k`$ regime scale, we have to take $`g1/\sqrt{N}`$, as we mentioned above. In Fig. 4 we plot $`W(k)`$ against $`k/\sqrt{g}`$. The small $`k`$ region scales as it should, and the results agree with the analytical prediction (4.2). Moreover the scaling extends up to $`k/\sqrt{g}=O(1)`$.

If the model is equivalent to ordinary gauge theory — namely to 4D pure super Yang-Mills theory with four supercharges — which is confining, then the Wilson loop should exhibit an area law behavior. In order to illustrate this behavior, we show a logarithmic plot of $`W(k)`$ versus the area $`k^2/g`$ in Fig. 4. In this figure only, we fine-tune $`g`$ as a function of $`N`$ so that the scaling in the intermediate regime of $`k`$ becomes even better. We stay with the convention $`g(48)=1`$ and use the optimal values $`g(32)=1.291`$, $`g(24)=1.563`$, $`g(16)=1.929`$, which is not far from $`g1/\sqrt{N}`$. The small deviation can be understood as a manifestation of finite $`N`$ effects. Fig. 4 shows indeed a region of $`k`$ that corresponds to the area law behavior $`W(k)\mathrm{exp}(\text{const.}k^2)`$. Surprisingly, the area law behavior is also observed in the bosonic model, as Fig. 12 shows, which is quite contrary to what one might have expected . In both cases, supersymmetric and bosonic, it is not clear from the data whether the area law extends to $`k=\mathrm{}`$ in the large $`N`$ limit. We will discuss the observed area law behavior from a theoretical point of view later.

We now proceed to the one-point function of the Polyakov line. In the 2D Eguchi-Kawai model  this quantity vanishes due to $`(\text{}_N)^D`$ symmetry. In the present model, however, there is no exact symmetry that could make such a quantity vanish, as we explained in Section 2. Note for instance that $`P(k=0)=1`$ for any configuration. Fig. 5 shows the results for $`P(k)`$. It falls off rapidly as $`k`$ increases 444One may consider the small $`k`$ expansion here, as in eq. (4.2). The result is $`P(k)=1\frac{1}{2N}k^2\text{tr}(X_1^2)+\mathrm{}`$. The fact that $`\text{tr}(A_\mu ^2)`$ is logarithmically divergent means that actually $`P(k)`$ has a non-analytic behavior $`1+\mathrm{const}.k^2\mathrm{ln}|k|`$ around $`k=0`$.. Again we observe a good scaling with $`g1/\sqrt{N}`$. Remember also that $`P(k)`$ is actually just a Fourier transform of the eigenvalue distribution. Therefore, the value of $`k`$ at which $`P(k)`$ drops to zero, which we denote as $`k_0`$, should be inversely proportional to the space-time extent $`R_{\mathrm{new}}`$. The observed scaling with $`g1/\sqrt{N}`$ is consistent with our result in the previous section, $`R_{\mathrm{new}}\sqrt{g}N^{1/4}`$. The result for the bosonic case is shown in Fig. 13. We obtain a similar behavior except for some oscillations in the large $`k`$ region. In particular, scaling is confirmed with $`g1/\sqrt{N}`$. The value of $`k_0`$ is larger than the supersymmetric $`k_0`$, as expected. The ratio of $`k_0`$ in the two cases is indeed roughly the inverse of the corresponding ratio of $`R_{\mathrm{new}}`$ (the bosonic $`k_0`$ is about twice as large as the supersymmetric one).

The above observations concerning $`P(k)`$ and $`R_{\mathrm{new}}`$ have an interesting implication on the Eguchi-Kawai equivalence. We recall that from the results for $`W(k)`$, we phenomenologically concluded that the Eguchi-Kawai equivalence holds at least in a finite range of scale for both, the supersymmetric and the bosonic case. We would like to understand this from a theoretical point of view. As we mentioned in Section 2, in the proof of Eguchi-Kawai equivalence, $`P(k)`$ is assumed to vanish. We have found that $`P(k)`$ is indeed small for $`k>k_0`$, but not for $`k<k_0`$. This means that the proof works for $`k>k_0`$, but not for small $`k`$, which corresponds to the ultraviolet regime in the corresponding gauge theory. We also observed that $`k_0`$ remains finite with respect to a physical scale in the large $`N`$ limit. A complementary understanding can be obtained by taking Gross-Kitazawa’s point of view . As explained in Ref. , the extent of the eigenvalue distribution of $`A_\mu `$ determines the momentum cutoff of the corresponding gauge theory . The observation in Section 3 that $`R_{\mathrm{new}}\sqrt{g}N^{1/4}`$ implies that the momentum cutoff remains finite in physical scale as $`N\mathrm{}`$. Let us assume that the momentum cutoff is finite, but large enough to attract the renormalization flow to the fixed point which corresponds to the universality class of gauge theory. Then the flow follows closely the renormalization trajectory of the gauge theory, in a certain regime. That would explain why the equivalence holds at least in a finite range of scale. However, since the momentum cutoff does not go to infinity in the large $`N`$ limit, it is conceivable that the renormalization flow will leave the renormalization trajectory of the gauge theory at some low-energy scale eventually. In this case the observed area law would not extend to $`k=\mathrm{}`$ even in the large $`N`$ limit.

### 4.2 Multi-point functions and universal scaling

In this subsection we proceed to the large $`N`$ scaling of multi-point functions of Wilson loops. We first note that in the bosonic case, there are analytical results to all order in the $`1/D`$ expansion . The statement is that

$$𝒪_1𝒪_2\mathrm{}𝒪_n_{con}O\left(\frac{1}{N^{2(n1)}}\right)\text{for the bosonic case},$$
(4.3)
where $`𝒪_i`$ denotes a Wilson loop or a Polyakov line as defined in eq. (4.1), and $`\mathrm{}_{con}`$ means that only the connected part is taken. The correlation functions should be considered as functions of $`k_{\text{phys}}=k/\sqrt{g}`$, where $`g`$ is taken to be proportional to $`1/\sqrt{N}`$. Our results for the bosonic model shown in Figs. 11 to 17 clearly confirm this analytical prediction. Let us consider a wave-function renormalization for each operator, $`𝒪_i^{(\mathrm{ren})}=Z𝒪_i`$, so that connected correlation functions of the renormalized operators $`𝒪_i^{(\mathrm{ren})}`$ become finite in the large $`N`$ limit. Relation (4.3) means, however, that we cannot make all the multi-point functions finite. If we make the two-point functions finite by choosing $`ZO(N)`$, then all the higher-point functions vanish in the large $`N`$ limit. In the supersymmetric case, we will see that scaling is observed again with $`g1/\sqrt{N}`$, but in contrast to the bosonic case a universal $`Z`$ that makes all the correlators finite seems to exist. In the following, we always set $`Z(N=48)=1`$, without loss of generality.

Let us start with the two-point functions, for which we measure the following two correlation functions,

$`G_2^{(W)}(k)`$ $`=`$ $`\{\text{Im}W(k)\}^2`$
$`G_2^{(P)}(k)`$ $`=`$ $`\{\text{Im}P(k)\}^2.`$ (4.4)
We take the imaginary part in order to avoid subtraction of a disconnected part. 555We also measured a number of multi-point functions, which are not presented here since the relative errors are rather large. (Note in this regard that since $`\text{Im}W(k)`$ and $`\text{Im}P(k)`$ are parity odd, the one-point functions $`\text{Im}W(k)`$ and $`\text{Im}P(k)`$ vanish due to parity invariance of the model.) The results are shown in Figs. 6 and 7, respectively. If we multiply the data by $`(N/48)^2`$, they scale nicely with $`g1/\sqrt{N}`$.

As a three-point function, we measure

$$G_3^{(W)}(k)=(\text{Im}W(k))^2\text{Re}W(k)(\text{Im}W(k))^2\text{Re}W(k).$$
(4.5)
We multiply the data either by $`(N/48)^3`$, which is required for the universal scaling of all the multi-point correlation functions, or by $`(N/48)^4`$, which is the factor predicted for the bosonic model. The results are compared in Fig. 8. We do observe a nice scaling behavior with a factor of $`(N/48)^3`$, but the scaling becomes worse for a factor of $`(N/48)^4`$.

Similarly, as a four-point function we measure

$$G_4^{(W)}(k)=(\text{Im}W(k))^43(\text{Im}W(k))^2.$$
(4.6)
We multiply the data either by $`(N/48)^4`$, which is required for the universal scaling of all the multi-point correlation functions, or by $`(N/48)^6`$, which is the factor predicted for the bosonic model. The results are compared in Fig. 9. Again the scaling behavior obtained with the factor for universal scaling is superior over the behavior with the bosonic factor.

To summarize our results concerning Wilson loop correlators, we observe that

$`𝒪O(1)`$ (4.7)
$`𝒪_1𝒪_2\mathrm{}𝒪_n_{con}O\left({\displaystyle \frac{1}{N^n}}\right)\text{for}n2.`$
These correlators scale as functions of $`k_{\text{phys}}=k/\sqrt{g}`$, where $`g`$ is taken to be proportional to $`1/\sqrt{N}`$. This means that all the multi-point functions of the renormalized operators $`𝒪_i^{(\mathrm{ren})}=Z𝒪_i`$ become finite in the large $`N`$ limit if we set $`ZO(N)`$, in contrast to the bosonic case. We will discuss further the implications of this universal scaling behavior in the next subsection.

Finally we comment on large the $`N`$ factorization. In ordinary gauge theory, large $`N`$ factorization can be shown by weak-coupling expansion as well as strong-coupling expansion. In a large $`N`$ reduced model with Hermitian matrices, one cannot do a weak-coupling or a strong-coupling expansion, because $`g`$ is not a coupling constant but a scale parameter, as we have mentioned. Hence large $`N`$ factorization is nontrivial. In the bosonic case, large $`N`$ factorization holds to all orders of the $`1/D`$ expansion . Our observation (4.7) implies

$$𝒪_1𝒪_2\mathrm{}𝒪_n=𝒪_1𝒪_2\mathrm{}𝒪_n+O\left(\frac{1}{N^2}\right),$$
(4.8)
where the $`O\left(1/N^2\right)`$ contributions are due to $`𝒪_1𝒪_2_{con}𝒪_3\mathrm{}𝒪_n`$, etc. Therefore the large $`N`$ factorization is also valid in the supersymmetric case.

### 4.3 Interpretation of the large $`N`$ scaling

In this subsection, we further discuss the physical significance of the large $`N`$ scaling (4.7) we observed.

If one views the supersymmetric matrix model considered here as a non-perturbative definition of type IIB string theory, and the Wilson loops as fundamental strings, string unitarity requires a large $`N`$ behavior of the form $`N^{a\chi _n}`$ for the connected correlators of $`n`$ Wilson loops, where $`\chi _n=2n`$ is the Euler characteristic of the worldsheet. In order to compare our results for the supersymmetric case (4.7) as well as for the bosonic case (4.3) to this behavior, we first drop the extra $`1/N^n`$ factor in (4.3) and (4.7), which is due to the chosen normalization (4.1) of the operators $`𝒪_i`$. Then we find that the connected correlators of Wilson loops change from an $`O(N^{\chi _n})`$ behavior to an $`O(1)`$ behavior by the introduction of supersymmetry. Our results for the supersymmetric case indicate $`a=0`$, which is not in contradiction to the above requirement of string unitarity, but it is an extreme case where one is far away from a perturbative expansion in genus. It would be interesting to understand the result $`a=0`$ analytically, directly from the matrix model. From the string theoretical point of view this indicates that the supersymmetric matrix model might automatically realize a kind of double scaling limit, as it has been known from (ordinary) matrix models of two-dimensional quantum gravity. In these models contributions from all genera are important.

While we do not presently have an analytic understanding of the difference between the connected correlators in the supersymmetric and the bosonic case, we can try to make an educated guess, based on the perturbative expansion (3.2).

Let us consider the bosonic case. After integrating out the off-diagonal elements perturbatively, schematically each diagram gives a contribution

$$\underset{i_1,\mathrm{},i_F}{}\left(\frac{g^2}{(x_ix_j)^2}\right)^L\left((x_kx_l)\frac{1}{g^2}\right)^{V_3}\left(\frac{1}{g^2}\right)^{V_4},$$
(4.9)
where $`x_i`$ denotes the diagonal elements as usual, and $`F`$, $`L`$, $`V_3`$, $`V_4`$ are the number of index loops (faces), propagators (links), 3-point and 4-point vertices, respectively. They obey the relations

$$F+VL=\chi ,4V_4+3V_3=2L,V=V_3+V_4.$$

Here $`\chi `$ is the Euler characteristic of the diagram given as $`\chi =22hn`$, where $`h`$ is the genus (the number of handles in the diagram) and $`n`$ is the number of the operators.

Let us now integrate over the diagonal elements $`x_i`$, under the assumption that the infrared region dominates. This assumption implies that the $`x_i`$’s are not close. In fact we assume that they are generically separated by some scale $`R`$, which is a measure of the extension of our universe. Viewing the $`x_i`$’s as the coordinates of the worldsheet in target space, the above assumption implies that the worldsheet is rough: points are scattered quite randomly. We can now estimate the value for the diagram considered by simply replacing the integration over the $`x_i`$’s by their typical separation, assuming the existence of a suitable measure which implements the above hypothesis. We obtain

$$N^F\left(\frac{g^2}{R^2}\right)^L\left(\frac{R}{g^2}\right)^{V_3}\left(\frac{1}{g^2}\right)^{V_4}N^\chi ,$$
(4.10)
where we have used $`R\sqrt{g}N^{1/4}`$ . So the assumption of infrared dominance reproduces the observed large $`N`$ behavior.

In the supersymmetric case, the fermion diagonal elements make the power-counting more complicated. However, let us naively consider the contribution (4.9). When we integrate over the diagonal elements $`x_i`$, let us assume that the ultraviolet rather than the infrared region dominates. Namely, we assume that smooth worldsheets are favoured dynamically due to supersymmetry 666It should be remarked here that such an increased smoothness of the bosonic part in a supersymmetric worldsheet has been observed in toy models for superstrings .. More precisely we assume that the dominant $`x_i`$ configurations are such that $`x_i`$’s, which are connected in a given diagram, are as close together as the actual distribution of $`x_i`$’s allows. Above we have seen that there seems to be such a minimal length $`\mathrm{}`$, which is of the order $`\sqrt{g}`$, characterizing the distribution of $`|x_ix_j|`$. Replacing $`|x_ix_j|`$ with this minimal length in (4.9) we arrive at

$$\left(\frac{g^2}{\mathrm{}^2}\right)^L\left(\frac{\mathrm{}}{g^2}\right)^{V_3}\left(\frac{1}{g^2}\right)^{V_4}O(1).$$
(4.11)
Thus the assumption of ultraviolet dominance leads to our observed results in the supersymmetric case. Of course it remains to be understood why there is a difference between infrared and ultraviolet dominance in the bosonic and supersymmetric cases.

If the above naive argument is true, it also implies that the contributions in the supersymmetric case do not depend on the genus $`h`$ either, and consequently that diagrams of all topologies contribute with equal weight. This is reminiscent of the double scaling limit of matrix models. For example, let us consider a Hermitian one-matrix model with the partition function

$$Z=\text{d}\varphi \text{e}^{N(\text{tr}\varphi ^2\lambda \text{tr}\varphi ^3)},$$
(4.12)
$`\varphi `$ being a $`N\times N`$ Hermitian matrix. The $`n`$-point correlation functions of loop operators behave as

$$\text{tr}\varphi ^{l_1}\mathrm{}\text{tr}\varphi ^{l_n}(Nϵ^{5/2})^{22hn}ϵ^n,$$
(4.13)
where $`ϵ`$ is the parameter related to the $`l_i`$ and to the coupling constant $`\lambda `$ as

$$l_i\frac{1}{ϵ},\lambda \lambda _cϵ^2,$$
(4.14)
and $`\lambda _c`$ is the critical coupling constant. Note that the large $`N`$ behavior for fixed $`ϵ`$ is $`N^\chi `$, which we obtain for the bosonic large $`N`$ reduced model. When one takes the large $`N`$ limit with $`Nϵ^{5/2}=g_{\mathrm{str}}^1`$ fixed, the dependence on $`h`$ disappears and the remaining power behavior $`ϵ^n`$ can be absorbed into the wavefunction renormalization of each operator. In the large $`N`$ reduced model, we do not have the parameter corresponding to $`\lambda `$. The large $`N`$ scaling we observed for the supersymmetric case is formally the same as one obtains in the double scaling limit of the one-matrix model. In this sense, one might say that the double scaling limit is taken automatically in the supersymmetric large $`N`$ reduced model, and that the string coupling constant $`g_{\mathrm{str}}`$ is not a tunable parameter but is fixed dynamically. If this is the case, it should be considered as a very satisfactory feature of the supersymmetric large $`N`$ reduced model as a nonperturbative definition of superstring theory, since we do expect the string coupling constant $`g_{\mathrm{str}}`$, which is related to the vacuum expectation value of the dilaton field, to be fixed dynamically (if superstring theory is treated nonperturbatively). It is also interesting that the qualitative difference of the large $`N`$ behavior for the bosonic and the supersymmetric case might be traced back to the smoothness of the worldsheet, which is itself a dynamical question to be addressed. This point needs further clarification.

## 5 Summary and discussion

In this paper, we have studied the large $`N`$ dynamics of a supersymmetric large $`N`$ reduced model by means of Monte Carlo simulations.

We studied the space-time structure represented by the eigenvalues of the bosonic matrices. In particular, we found that the large $`N`$ power behavior of the space-time extent is consistent with the branched-polymer picture based on the one-loop perturbative expansion around diagonal matrices. The effect of fermions in the space-time extent was observed by the enhancement of the coefficient in the power behavior, but not in the power itself. The power appears to be the same for the bosonic and supersymmetric case. We emphasized, however, that the theoretical explanation is completely different. We also found that the space-time uncertainty is clearly reduced for the supersymmtric case, which means that space-time comes closer to the classical behavior. Even in the supersymmetric case, the space-time uncertainty is found to be finite in the physical scale in the large $`N`$ limit. We argued that this implies that the model satisfies the uncertainty principle for the nonperturbative definition of superstring theory.

The large $`N`$ scaling behavior of Wilson loop correlators is observed at fixed $`g^2N`$. Although this scaling of $`g`$ is the same as in the bosonic model, there is a striking difference from the bosonic case in the wave-function renormalization with the multi-point functions. In the bosonic case, there was no universal scaling behavior: keeping two-point functions finite, all the higher-point functions vanish. In the supersymmetric case, we observed a clear trend for all the higher-point functions to become finite in the large $`N`$ limit. We gave a perturbative argument that this result for the supersymmetric case might be understood if we assume smooth worldsheets to dominate. This argument also implies that all the topologies of the worldsheet contribute with equal weight to the amplitude. All these features are reminiscent of the double scaling limit of matrix models.

We also addressed the issue of Eguchi-Kawai equivalence. By searching for the area law behavior in the one-point function of the Wilson loop, we concluded that the equivalence does hold at least in a finite region of scale. What is rather surprising is that the area law behavior has been observed also for the bosonic model. This suggests that the bosonic model is also equivalent to ordinary large $`N`$ Yang-Mills theory at least in a finite region of scale, which is contrary to what has been generally believed. We argued, however, that this conclusion can be understood from a more theoretical point of view based on the large $`N`$ behavior obtained for $`R_{\mathrm{new}}`$ and the one-point function of the Polyakov line. It is an open question whether this equivalence extends to the far infrared regime.

To summarize, we have gained new insight into the dynamical properties of the large $`N`$ behavior of a supersymmetric large $`N`$ reduced model. We hope that our findings shed light on the dynamical aspects of the most interesting 10D version of our model, i.e. the IIB matrix model. In this respect, it is encouraging that the large $`N`$ scaling of Wilson loop correlators in the present model has been observed at fixed $`g^2N`$, which coincides with the result obtained by requiring that the loop equations of the IIB matrix model should reproduce the string field Hamiltonian. We presume that a large $`N`$ scaling of Wilson loop correlators — like the one we observed — also holds for the IIB matrix model; then the only difference would be the spontaneous breakdown of Lorentz symmetry. One of the good news revealed in the present work is that low energy effective theory, based on the one-loop approximation, does already capture the low energy dynamics of the supersymmetric matrix model. We therefore hope to address the most interesting issue of spontaneous breakdown of Lorentz invariance by using the low energy effective theory — which is in 10D far more complicated than in the 4D case. We are going to report on Monte Carlo studies of IIB matrix model along these lines in the near future.

## Acknowledgment

We thank T. Ishikawa, C.F. Kristjansen, Y.M. Makeenko, T. Nakajima, M. Staudacher, A. Tsuchiya and G. Vernizzi for valuable discussions. J. N. is supported by the Japan Society for the Promotion of Science as a Postdoctoral Fellow for Research Abroad. The computation has been done partly on Fujitsu VPP500 at High Energy Accelerator Research Organization (KEK), Fujitsu VPP700E at The Institute of Physical and Chemical Research (RIKEN), and NEC SX4 at Research Center for Nuclear Physics (RCNP) of Osaka University. This work is supported by the Supercomputer Project (No.99-53) of KEK.

## Appendix A The algorithm for the Monte Carlo simulation

In this appendix, we explain the algorithm we use for the Monte Carlo simulation of the supersymmetric matrix model. Only in this appendix we set $`g=1`$ for simplicity.

We first carry out the integration over fermionic matrices to obtain the explicit formula for the fermion determinant. We calculate

$$Z_f[A]=\text{d}\psi \text{d}\overline{\psi }\text{e}^{S_f},$$
(A.1)
where we use the notation introduced in eq. (2.1). We define a set of generators $`t^a`$ gl($`N`$,) by

$$(t^a)_{ij}=\delta _{ii_a}\delta _{jj_a}(a=1,\mathrm{},N^2),$$
(A.2)
where $`i_a`$ and $`j_a`$ are integers running from 1 to $`N`$, specified uniquely by

$$a=N(i_a1)+j_a.$$
(A.3)
We also introduce the notation $`\overline{a}=N(j_a1)+i_a`$. The fermionic matrix $`\psi _\alpha `$ can be expanded in terms of $`t^a`$ as

$$(\psi _\alpha )_{ij}=\underset{a=1}{\overset{N^2}{}}\psi _{a\alpha }(t^a)_{ij},$$
(A.4)
where $`\psi _{a\alpha }=(\psi _\alpha )_{i_aj_a}`$. $`\overline{\psi }_\alpha `$ and $`A_\mu `$ can be expanded similarly with the coefficients $`\overline{\psi }_{a\alpha }=(\overline{\psi }_\alpha )_{i_aj_a}`$ and $`A_{a\mu }=(A_\mu )_{i_aj_a}`$. Note also that $`A_{\overline{a}\mu }=(A_{a\mu })^{}`$ due to the Hermiticity of $`A_\mu `$.

We define the structure constants $`g_{abc}`$ of gl($`N`$,) by

$`g_{abc}`$ $`=`$ $`\text{tr}(t^c[t^a,t^b])`$ (A.5)
$`=`$ $`\delta _{j_ai_b}\delta _{j_bi_c}\delta _{j_ci_a}\delta _{j_ci_b}\delta _{j_bi_a}\delta _{j_ai_c}.`$
The fermionic action then reads

$`S_f`$ $`=`$ $`g_{abc}\overline{\psi }_{c\alpha }(\mathrm{\Gamma }_\mu )_{\alpha \beta }A_{a\mu }\psi _{b\beta }`$ (A.6)
$`=`$ $`\overline{\psi }_{a\alpha }_{a\alpha ,b\beta }^{}\psi _{b\beta },`$
where

$$_{a\alpha ,b\beta }^{}=g_{abc}(\mathrm{\Gamma }_\mu )_{\alpha \beta }A_{c\mu }.$$
(A.7)
We first integrate out $`(\psi _\alpha )_{NN}`$ and $`(\overline{\psi }_\alpha )_{NN}`$ using the $`\delta `$ functions in the measure (2.2). We get a factor of $`1/N^4`$ followed by the replacements

$$(\psi _\alpha )_{NN}\underset{j=1}{\overset{N1}{}}(\psi _\alpha )_{jj};(\overline{\psi }_\alpha )_{NN}\underset{j=1}{\overset{N1}{}}(\overline{\psi }_\alpha )_{jj}$$
(A.8)
in the fermionic action. The integration over the remaining Grassmann variables yields $`det`$, where $``$ is the $`\mathrm{\hspace{0.17em}2}(N^21)`$ $`\times `$ $`2(N^21)`$ complex matrix defined by

$$_{a\alpha ,b\beta }=_{a\alpha ,b\beta }^{}_{N^2\alpha ,b\beta }^{}\delta _{i_aj_a}_{a\alpha ,N^2\beta }^{}\delta _{i_bj_b}$$
(A.9)
(the indices $`a`$ and $`b`$ run from 1 to $`N^21`$). Thus, we obtain

$$Z_f[A]=\frac{1}{N^4}det.$$
(A.10)
We first want to show that the determinant $`det`$ is real positive 777This has been already reported in Ref.  as a numerical observation. For related work, see Ref. .. For this purpose we note that the matrix $``$ satisfies the identity $`\sigma _2\sigma _2=^{}`$. Hence if $`\phi _{a\alpha }`$ is an eigenvector of $``$ with an eigenvalue $`\lambda `$, then $`\psi _{a\alpha }=(\sigma _2)_{\alpha \beta }(\phi _{a\beta })^{}`$ is an eigenvector of $``$ with an eigenvalue $`\lambda ^{}`$. It is important that the two vectors $`\phi _{a\alpha }`$, $`\psi _{a\alpha }`$ are linearly independent. The determinant, which is the product of all the eigenvalues of $``$, should therefore be real and positive semi-definite. In the case of 6D or 10D (IIB matrix model) versions of the supersymmetric large $`N`$ reduced model, the fermion integral yields a complex effective action in general. This causes the notorious sign problem, which makes standard Monte Carlo simulations practically inapplicable for large $`N`$. In the present case, since the determinant $`det`$ is real positive, we can introduce a $`\mathrm{\hspace{0.17em}2}(N^21)`$ $`\times `$ $`2(N^21)`$ Hermitian matrix $`𝒟=^{}`$, which has real positive eigenvalues, and $`det=\sqrt{det𝒟}`$. Therefore we have written the effective action for the bosonic matrices $`A_\mu `$ in eq. (2.11) as

$$S_{\text{eff}}=S_b\frac{1}{2}\mathrm{ln}det𝒟.$$
(A.11)
We apply the Hybrid R algorithm  to simulate this system 888Ref.  gives an overview of effective algorithms for dynamical fermions, including the Hybrid R algorithm..

The first step of the Hybrid R algorithm is to apply the molecular dynamics method . We introduce a conjugate momentum for $`A_{a\mu }`$ as $`X_{a\mu }`$, which satisfies $`X_{\overline{a}\mu }=(X_{a\mu })^{}`$. The partition function can be re-written as

$$Z=\text{d}X\text{d}A\text{e}^H,$$
(A.12)
where $`H`$ is the “Hamiltonian” defined by

$$H=\frac{1}{2}\underset{\mu a}{}X_{\overline{a}\mu }X_{a\mu }+S_b[A]\frac{1}{2}\mathrm{ln}det𝒟.$$
(A.13)
The update of $`X_{a\mu }`$ can be done by just generating $`X_{a\mu }`$ with the probability distribution $`\mathrm{exp}(\frac{1}{2}|X_{a\mu }|^2)`$. In order to update $`A_{a\mu }`$, we use the Hamiltonian equations

$`{\displaystyle \frac{\text{d}A_{a\mu }(\tau )}{\text{d}\tau }}`$ $`=`$ $`{\displaystyle \frac{H}{X_{a\mu }}}=X_{\overline{a}\mu },`$ (A.14)
$`{\displaystyle \frac{\text{d}X_{a\mu }(\tau )}{\text{d}\tau }}`$ $`=`$ $`{\displaystyle \frac{H}{A_{a\mu }}}={\displaystyle \frac{1}{2}}\text{tr}\left({\displaystyle \frac{𝒟}{A_{a\mu }}}𝒟^1\right){\displaystyle \frac{S_b}{A_{a\mu }}}.`$ (A.15)
Along the “classical trajectory” given by the Hamiltonian equation,

(i) $`H`$ is invariant,

(ii) the motion is reversible,

(iii) the phase-volume is preserved,

$$\frac{(A(\tau ),X(\tau ))}{(A(0),X(0))}=1,$$
(A.16)
where $`(A(\tau ),X(\tau ))`$ is a point on the trajectory after evolution from $`(A(0),X(0))`$. Therefore, generating a new set of $`(A,X)`$ by solving the Hamiltonian equation for a fixed “time” interval $`\tau `$ satisfies detailed balance. This procedure — together with the proceeding generation of $`X_{a\mu }`$ with the Gaussian distribution — is called “one trajectory”, which corresponds to “one sweep” in ordinary Monte Carlo simulations.

In order to solve the Hamiltonian equation numerically, we have to discretize the “time” $`\tau `$. A discretization which maintains the properties (ii) and (iii) is known. The slight violation of (i) for finite $`\mathrm{\Delta }\tau `$ causes systematic errors. One can in principle eliminate the systematic error completely, by making a Metropolis accept/reject decision at the end of each trajectory. But in the present case, the overhead for this procedure is rather large. We therefore decided to omit that step, and just use a sufficiently small $`\mathrm{\Delta }\tau `$. Still we can use the specific discretization of Ref. , which we explain below, to minimize the systematic error. As we explain later, we do find a good convergence in small $`\mathrm{\Delta }\tau `$, and the systematic error is well under control.

We introduce a short-hand notation for the discretized $`X_{a\mu }(\tau )`$ and $`A_{a\mu }(\tau )`$,

$$X_{a\mu }^{(r)}=X_{a\mu }(r\mathrm{\Delta }\tau );A_{a\mu }^{(s)}=A_{a\mu }(s\mathrm{\Delta }\tau ).$$
(A.17)
The Hamiltonian equations are discretized as

$`A_{a\mu }^{(\frac{1}{4})}`$ $`=`$ $`A_{a\mu }^{(0)}+{\displaystyle \frac{\mathrm{\Delta }\tau }{4}}X_{\overline{a}\mu }^{(0)}`$
$`A_{a\mu }^{(n+\frac{1}{2})}`$ $`=`$ $`A_{a\mu }^{(n+\frac{1}{4})}+{\displaystyle \frac{\mathrm{\Delta }\tau }{4}}X_{\overline{a}\mu }^{(n)}`$
$`A_{a\mu }^{(m+\frac{1}{4})}`$ $`=`$ $`A_{a\mu }^{(m\frac{1}{2})}+{\displaystyle \frac{3\mathrm{\Delta }\tau }{4}}X_{\overline{a}\mu }^{(m)}`$
$`A_{a\mu }^{(\nu )}`$ $`=`$ $`A_{a\mu }^{(\nu \frac{1}{2})}+{\displaystyle \frac{\mathrm{\Delta }\tau }{2}}X_{\overline{a}\mu }^{(\nu )}`$
$`X_{a\mu }^{(n+1)}`$ $`=`$ $`X_{a\mu }^{(n)}+\mathrm{\Delta }\tau \left\{{\displaystyle \frac{1}{2}}R_{a\mu }^{(n+\frac{1}{2})}{\displaystyle \frac{S_b}{A_{a\mu }}}\left(A_{a\mu }^{(n+\frac{1}{2})}\right)\right\},`$ (A.18)
where $`n=0,1,\mathrm{},\nu 1`$, $`m=1,\mathrm{},\nu 1`$, and $`R_{a\mu }^{(n+\frac{1}{2})}`$ is defined by

$`R_{c\mu }^{(n+\frac{1}{2})}`$ $`=`$ $`\mathrm{\Phi }_{a\alpha }^{}\left({\displaystyle \frac{𝒟(A_{a\mu }^{(n+\frac{1}{2})})}{A_{c\mu }}}\right)_{a\alpha b\beta }\mathrm{\Phi }_{b\beta },`$ (A.19)
$`𝒟(A_{a\mu }^{(n+\frac{1}{2})})\mathrm{\Phi }`$ $`=`$ $`^{}(A_{a\mu }^{(n+\frac{1}{4})})\eta .`$ (A.20)
Here $`\eta _{a\alpha }`$ are complex variables generated by the Gaussian distribution $`\mathrm{exp}(_{a\alpha }|\eta _{a\alpha }|^2)`$. The judicious choice of the argument of $`^{}`$ is the tool to reduce the systematic error .

We solve eq. (A.20) with respect to $`\mathrm{\Phi }`$ by means of the conjugate gradient method , which is iterative. Each iteration involves a multiplication of the matrix $`𝒟`$ with some vector $`v`$. Since $`𝒟`$ is a $`\mathrm{\hspace{0.17em}2}(N^21)`$ $`\times `$ $`2(N^21)`$ matrix, storing $`𝒟`$ requires $`O(N^4)`$ memory, and multiplying $`𝒟`$ with $`v`$ naively involves $`O(N^4)`$ arithmetic operations. Actually we can do much better than this. We first recall that $`𝒟=^{}`$, where $``$ is the $`\mathrm{\hspace{0.17em}2}(N^21)`$ $`\times `$ $`2(N^21)`$ matrix defined in eq. (A.9). The point is that the number of nonzero elements of $``$ is only $`O(N^3)`$ (not $`O(N^4)`$). Indeed, the multiplication $`v`$ can be done economically as follows.

We consider

$$w_{a\alpha }=_{a\alpha b\beta }v_{b\beta },$$
(A.21)
and define the quantities $`w_{a\alpha }^{}`$ and $`v_{a\alpha }^{}`$, where $`a`$ runs from 1 to $`N^2`$ as in $`^{}`$, by

$`v_{a\alpha }^{}`$ $`=`$ $`v_{a\alpha }\text{for}a=1,\mathrm{},N^21,`$
$`v_{N^2\alpha }^{}`$ $`=`$ $`{\displaystyle \underset{i_a=j_a}{}}v_{a\alpha },`$ (A.22)
$`w_{a\alpha }^{}`$ $`=`$ $`_{a\alpha b\beta }^{}v_{b\beta }^{}.`$ (A.23)
Now $`w_{a\alpha }`$ can be written as

$$w_{a\alpha }=\{\begin{array}{cc}w_{a\alpha }^{}w_{N^2\alpha }^{}\hfill & \text{for}i_a=j_a\hfill \\ w_{a\alpha }^{}\hfill & \text{otherwise.}\hfill \end{array}$$
(A.24)
Thus the problem reduces to calculating the matrix-vector product in eq. (A.23). Using definition (A.7), we obtain

$$(w_\alpha ^{})_{ij}=(\mathrm{\Gamma }^\mu )_{\alpha \beta }[A_\mu ,v_\beta ^{}]_{ji},$$
(A.25)
where $`w_\alpha ^{}`$ and $`v_\alpha ^{}`$ are $`N\times N`$ matrices associated with $`w_{a\alpha }^{}`$ and $`v_{a\alpha }^{}`$, respectively, as in eq. (A.4). The commutator in eq. (A.25) requires $`O(N^3)`$ arithmetic operations. Thus we save $`O(N)`$ operations. In addition, we do not have to store neither $`g_{abc}`$ nor $``$. Multiplication of $`^{}`$ with some vector $`v`$ is done in the same way.

A similar technique should be used to calculate $`R_{a\mu }`$ in eq. (A.19). Note first that it can be written as $`R_{c\mu }=T_{c\mu }+(T_{\overline{c}\mu })^{}`$, where $`T_{c\mu }`$ is given by

$`T_{c\mu }`$ $`=`$ $`\mathrm{\Psi }_{a\alpha }\left({\displaystyle \frac{}{A_{c\mu }}}\right)_{a\alpha b\beta }\mathrm{\Phi }_{b\beta },`$
$`\mathrm{\Psi }_{a\alpha }`$ $`=`$ $`(_{a\alpha b\beta }\mathrm{\Phi }_{b\beta })^{}.`$ (A.26)
We define $`\mathrm{\Phi }^{}`$ and $`\mathrm{\Psi }^{}`$ in terms of $`\mathrm{\Phi }`$ and $`\mathrm{\Psi }`$, as we defined $`v^{}`$ in terms of $`v`$ before in eq. (A). Now we can re-write $`T_{c\mu }`$ as

$$T_{c\mu }=\frac{}{A_{c\mu }}\left(\mathrm{\Psi }_{a\alpha }^{^{}}(^{})_{a\alpha b\beta }\mathrm{\Phi }_{b\beta }^{}\right).$$
(A.27)
Using again eq. (A.7), we obtain

$$(T_\mu )_{ij}=(\mathrm{\Gamma }_\mu )_{\alpha \beta }[\mathrm{\Psi }_\alpha ^{^{}},\mathrm{\Phi }_\beta ^{}]_{ji},$$
(A.28)
where $`\mathrm{\Phi }_\alpha ^{}`$, $`\mathrm{\Psi }_\alpha ^{}`$ and $`T_\mu ^{}`$ are $`N\times N`$ matrices associated with $`\mathrm{\Phi }_{a\alpha }^{^{}}`$, $`\mathrm{\Psi }_{a\alpha }^{^{}}`$ and $`T_{a\mu }^{^{}}`$, respectively, as in eq. (A.4).

There are two parameters $`\nu `$ and $`\mathrm{\Delta }\tau `$ in this algorithm. We can choose $`\nu \mathrm{\Delta }\tau `$ so that a typical autocorrelation time is minimized. We have taken $`\nu \mathrm{\Delta }\tau =1`$ throughout the present work, and $`\nu =200,\mathrm{\hspace{0.17em}280},\mathrm{\hspace{0.17em}500}`$ for each of the cases $`N=16,\mathrm{\hspace{0.17em}24},\mathrm{\hspace{0.17em}32}`$, and $`\nu =500,\mathrm{\hspace{0.17em}600}`$ for $`N=48`$. Except in Fig. 2, we observed that the results are reasonably well converged at $`\nu =500`$, $`\mathrm{\Delta }\tau =0.002`$, so we just present those results. For Fig. 2 we carried out an extrapolation to $`\mathrm{\Delta }\tau =0`$ by assuming the $`\mathrm{\Delta }\tau `$ dependence of some observables $`Q(\mathrm{\Delta }\tau )`$ to be

$$Q(\mathrm{\Delta }\tau )Q(\mathrm{\Delta }\tau =0)(\mathrm{\Delta }\tau )^2\text{tr}(A_\mu ^2)_{\mathrm{\Delta }\tau }.$$
(A.29)
This assumption has been checked for $`\text{tr}F^2`$ with the exact result (2.12). We also observed that $`\text{tr}(A_\mu ^2)_{\mathrm{\Delta }\tau }`$ behaves as

$$\text{tr}(A_\mu ^2)_{\mathrm{\Delta }\tau }c_1c_2\mathrm{log}\mathrm{\Delta }\tau ,$$
(A.30)
for small $`\mathrm{\Delta }\tau `$, where $`c_1`$ and $`c_2`$ are constants depending on $`N`$. 999In QCD the $`\mathrm{\Delta }\tau `$ dependence of the systematic error is $`O(\mathrm{\Delta }\tau ^2)`$ . A similar argument leads to the assumption (A.29). Due to eq. (A.30), the $`\mathrm{\Delta }\tau `$ dependence of the systematic error in our case is expected to be $`(\mathrm{\Delta }\tau )^2\mathrm{log}\mathrm{\Delta }\tau `$. This implies that it diverges logarithmically for $`\mathrm{\Delta }\tau 0`$, which is consistent with the theoretical prediction discussed below eq. (3.3).

Let us comment on the required computational effort of our algorithm. The dominant part comes from solving the linear system (A.20) using the conjugate gradient method. First of all, we find that the number of iterations necessary for the convergence of the method seems to grow linearly with the size of the matrix $`𝒟`$, namely as $`O(N^2)`$. This is much worse than the full QCD case with a fixed quark mass, where the number of iterations does not depend on the system size. We may interpret this phenomenon as a sort of “critical slowing down”, since the present system corresponds to QCD in the chiral limit. As we have seen, the number of arithmetic operations for each iteration is of order $`N^3`$. Therefore, the required computational effort of our algorithm is estimated to be $`O(N^5)`$.

For the bosonic case, we use the heat bath algorithm in the way proposed in Ref. , which requires an effort of $`O(N^4)`$. We note, however, that application of a Hybrid Monte Carlo algorithm  allows for an $`O(N^3)`$ algorithm for the bosonic case, which might be useful for proceeding to much larger $`N`$.

Finally, we comment on the numbers of configurations used for the measurements. For the supersymmetric case, they are 3060, 1508, 1296, 436 for $`N=16,24,32,48`$, respectively. For the bosonic case, we used 1000 configurations for each $`N`$.

## Appendix B Results for the bosonic case

For comparison we show in this appendix the results for the bosonic case. By the bosonic case we mean a model obtained by just dropping the fermions from the supersymmetric matrix model described by the partition function (2.1). Fig. 10 shows the distribution $`\rho (r)`$ defined in Section 3. Figs. 11 to 17 show the Wilson loop and Polyakov line correlators defined in Section 4. We take $`g1/\sqrt{N}`$ ($`g=1`$ for $`N=48`$) and plot the results against $`k_{\text{phys}}=k/\sqrt{g}`$, as in the supersymmetric case. We multiply the results by $`(N/48)^{2(n1)}`$ for $`n`$-point functions.

The data scale nicely in agreement with the theoretical prediction for large $`N`$ given by eq. (4.3). For comparison we also show the 3-point and the 4-point Wilson loop correlators with the renormalization factors, which were used successfully in Sec. 4 for the supersymmetric case. We see very clearly that the bosonic prediction is the correct one in this case.


format: text

text:
# 1 Introduction

## 1 Introduction

Since the NMC experiment in 1991, which precisely measured the structure functions of the proton and neutron, $`F_2^p(x)`$ and $`F_2^n(x)`$, for a wide region of Bjorken’s $`x`$, it has been realized that the Gottfried sum rule is violeted and the light sea-quark distributions, $`\overline{u}(x)`$ and $`\overline{d}(x)`$, are asymmetric. A considerable excess of the $`\overline{d}`$ quark density relative to the $`\overline{u}`$ quark density was seen. The same result was confirmed by the E866 experiment which measured the cross section ratio of the Drell-Yang processes, $`\sigma (p+d)/\sigma (p+p)`$, though the violation of the Gottfried sum rule is smaller than reported by the NMC. Now, study on the origin of the flavor asymmetry of light sea–quarks has been a challenging subject in particle and nuclear physics because it is closely related to the dynamics of nonperturbative QCD. Although the most widely accepted idea to understand it is the meson cloud model, many discussions are still under going with several approaches such as chiral quark model, Skyrme model, Pauli blocking effects, etc..

However, what is going on with the polarized case, $`\mathrm{\Delta }\overline{u}(x)`$ and $`\mathrm{\Delta }\overline{d}(x)`$? In these years, measurement of the polarized structure function of the nucleon in polarized deep–inelastic scatterings have shown that the nucleon spin is carried by quarks a little and the strange sea–quark is negatively polarized in quite large. The results were not anticipated by conventional theories and often referred as ‘the proton spin crisis. By using many data with high precision on the polarized structure functions of the proton, neutron and deuteron accumulated so far, good parametrization models of polarized parton distribution functions have been proposed at the next–to–leading order(NLO) of QCD. The behavior of polarized valence $`u`$ and $`d`$ quarks has been well–known from such analyses. However, the knowledge of polarized sea–quarks and gluons is still poor. Although people usually assume the symmetric light sea–quark polarized distribution, i.e. $`\mathrm{\Delta }\overline{u}(x)=\mathrm{\Delta }\overline{d}(x)`$, in analyzing the polarized structure functions of nucleons, there is no physical ground of such an assumption. In order to understand the nucleon spin structure, it is very important to know if the light sea–quark flavor symmetry is broken even for polarized distributions and to determine how $`\mathrm{\Delta }\overline{u}(x)`$ and $`\mathrm{\Delta }\overline{d}(x)`$ behave in the nucleon. Related to these subjects, it is interesting to know that even if we start with the symmetric distributions for the polarized light sea–quarks, $`\mathrm{\Delta }\overline{u}=\mathrm{\Delta }\overline{d}`$, at an initial $`Q_0^2`$, the symmetry can be violated for higher $`Q^2`$ regions, if the polarized distributions are perturbatively evolved in NLO calculations of QCD. In addition, some people have estimated the amount of its violation at an initial $`Q_0^2`$ using some effective models. However, their results do not agree with each other. Therefore, it is interesting to extract the value of $`\mathrm{\Delta }\overline{d}(x)\mathrm{\Delta }\overline{u}(x)`$ from the experimental data and test the flavor symmetry of $`\mathrm{\Delta }\overline{d}(x)`$ and $`\mathrm{\Delta }\overline{u}(x)`$ experimentally.

Recently, using longitudinal polarized lepton beams and longitudinal polarized fixed targets, SMC group at CERN and HERMES group at DESY observed the cross sections of the following semi–inclusive processes,

$$\stackrel{}{l}+\stackrel{}{N}l^{}+h+X,$$
(1)
and obtained the data on spin asymmetries for proton, deuteron and 3He targets, where $`h`$ is a created charged hadron or one of $`\pi ^\pm `$, $`K^\pm `$, $`p`$ and $`\overline{p}`$. A created hadron depends on the flavor of a parent quark and thus properly combining these data it is possible to decompose polarized quark distributions into the ones with individual flavor. These data provide a good material to test the light flavor symmetry of polarized sea–quark distributions and it might be timely to test the symmetry by using the present data.

## 2 Simple formulas of $`\mathrm{\Delta }\overline{d}(x)\mathrm{\Delta }\overline{u}(x)`$

In this letter, we propose new formulas for extracting a difference, $`\mathrm{\Delta }\overline{d}\mathrm{\Delta }\overline{u}`$, from the data of the above–mentioned semi–inclusive processes and estimate the value of it from the present data in order to test if the light flavor symmetry of polarized sea–quark distributions is originally violated.

Let us start with the semi–inclusive asymmetry for the process of eq.(1) with proton targets, which is written by

$$A_{1p}^h(x,Q^2)=\frac{_{q,H}e_q^2\{\mathrm{\Delta }q(x,Q^2)D_q^H(Q^2)+\mathrm{\Delta }\overline{q}(x,Q^2)D_{\overline{q}}^H(Q^2)\}}{_{q,H}e_q^2\{q(x,Q^2)D_q^H(Q^2)+\overline{q}(x,Q^2)D_{\overline{q}}^H(Q^2)\}}\times \{1+R(x,Q^2)\},$$
(2)
in the leading order(LO) of QCD, where $`\mathrm{\Delta }q(x,Q^2)`$ ($`\mathrm{\Delta }\overline{q}(x,Q^2)`$) and $`q(x,Q^2)`$ ($`\overline{q}(x,Q^2)`$) are the spin–dependent and spin–independent quark distribution functions at some values of $`x`$ and $`Q^2`$, respectively, and $`R(x,Q^2)`$ is a ratio of the absorption cross section of longitudinally and transversely polarized virtual photons by the nucleon, $`R(x,Q^2)=\sigma _L/\sigma _T`$. $`D_q^H(Q^2)`$ is given by integration of the fragmentation function, $`D_q^H(z,Q^2)`$, over the measured kinematical region of $`z`$, i.e. $`D_q^H(Q^2)=_{z_{min}}^1𝑑zD_q^H(z,Q^2)`$, where $`D_q^H(z,Q^2)`$ represents the probability of producing a hadron $`H`$ carrying momentum fraction $`z`$ at some $`Q^2`$ from a struck quark with flavor $`q`$. $`h`$ is the observed hadron concerned with here. When $`h`$ is $`h^+`$, the fragmentation function of, for example, $`u`$–quark decaying into $`h^+`$ is given by

$$D_u^{h^+}(z,Q^2)=D_u^{\pi ^+}(z,Q^2)+D_u^{K^+}(z,Q^2)+D_u^p(z,Q^2),$$
(3)
because $`h^+`$ is dominantly composed of $`\pi ^+`$, $`K^+`$ and $`p`$. Assuming the reflection symmetry along the V–spin axis, the isospin symmetry and charge conjugation invariance of the fragmentation functions, many fragmentation functions can be classified into the following 6 functions,

$`DD_u^{\pi ^+}=D_{\overline{d}}^{\pi ^+}=D_d^\pi ^{}=D_{\overline{u}}^\pi ^{},`$
$`\stackrel{~}{D}D_d^{\pi ^+}=D_{\overline{u}}^{\pi ^+}=D_u^\pi ^{}=D_{\overline{d}}^\pi ^{}=D_s^{\pi ^+}=D_{\overline{s}}^{\pi ^+}=D_s^\pi ^{}=D_{\overline{s}}^\pi ^{},`$
$`D^KD_u^{K^+}=D_{\overline{s}}^{K^+}=D_{\overline{u}}^K^{}=D_s^K^{},`$ (4)
$`\stackrel{~}{D^K}D_d^{K^+}=D_s^{K^+}=D_{\overline{u}}^{K^+}=D_{\overline{d}}^{K^+}=D_u^K^{}=D_d^K^{}=D_{\overline{d}}^K^{}=D_{\overline{s}}^K^{},`$
$`D^pD_u^p=D_d^p=D_{\overline{u}}^{\overline{p}}=D_{\overline{d}}^{\overline{p}},`$
$`\stackrel{~}{D^p}D_s^p=D_{\overline{u}}^p=D_{\overline{d}}^p=D_{\overline{s}}^p=D_u^{\overline{p}}=D_d^{\overline{p}}=D_s^{\overline{p}}=D_{\overline{s}}^{\overline{p}},`$
where $`D^H`$ and $`\stackrel{~}{D^H}`$ are called favored and unfavored fragmentation functions, respectively. Here we follow the commonly taken assumption on the fragmentation functions, for simplicity.

Now, we can rewrite eq.(2) as

$`{\displaystyle \underset{q,H}{}}e_q^2\{\mathrm{\Delta }q(x,Q^2)D_q^H(Q^2)+\mathrm{\Delta }\overline{q}(x,Q^2)D_{\overline{q}}^H(Q^2)\}`$
$`={\displaystyle \frac{A_{1p}^h(x,Q^2)[_{q,H}e_q^2\{q(x,Q^2)D_q^H(Q^2)+\overline{q}(x,Q^2)D_{\overline{q}}^H(Q^2)\}]}{\{1+R(x,Q^2)\}}}`$
$`=\mathrm{\Delta }N_p^h(x,Q^2),`$ (5)
where $`\mathrm{\Delta }N_p^h(x,Q^2)`$ is reffered to the spin–dependent production processes of charged hadrons with proton targets. From a combination of $`\mathrm{\Delta }N_{p,n}^{h^+,h^{}}(x,Q^2)`$ for proton and neutron targets, we can obtain the following formula,

$`\mathrm{\Delta }\overline{d}(x,Q^2)\mathrm{\Delta }\overline{u}(x,Q^2)`$
$`={\displaystyle \frac{\mathrm{\Delta }N_p^{h^+}(x,Q^2)\mathrm{\Delta }N_n^{h^+}(x,Q^2)\mathrm{\Delta }N_p^h^{}(x,Q^2)+\mathrm{\Delta }N_n^h^{}(x,Q^2)}{2I_1(Q^2)}}`$
$`{\displaystyle \frac{\mathrm{\Delta }N_p^{h^+}(x,Q^2)\mathrm{\Delta }N_n^{h^+}(x,Q^2)+\mathrm{\Delta }N_p^h^{}(x,Q^2)\mathrm{\Delta }N_n^h^{}(x,Q^2)}{2I_2(Q^2)}},`$ (6)
where

$`I_1(Q^2)=5D(Q^2)+4D^K(Q^2)+3D^p(Q^2)5\stackrel{~}{D}(Q^2)4\stackrel{~}{D^K}(Q^2)3\stackrel{~}{D^p}(Q^2),`$
$`I_2(Q^2)=3D(Q^2)+4D^K(Q^2)+3D^p(Q^2)+3\stackrel{~}{D}(Q^2)+2\stackrel{~}{D^K}(Q^2)+3\stackrel{~}{D^p}(Q^2).`$
Furthermore, if one can specify the detected charged hadron in experiment, one can obtain more simplified formulas for the difference of polarized light sea–quark densities. For the case of semi–inclusive $`\pi ^\pm `$–productions with proton and neutron targets, the difference can be written by

$`\mathrm{\Delta }\overline{d}(x,Q^2)\mathrm{\Delta }\overline{u}(x,Q^2)`$
$`={\displaystyle \frac{1}{6\{D(Q^2)+\stackrel{~}{D}(Q^2)\}}}\times [\{J(Q^2)1\}\{\mathrm{\Delta }N_p^{\pi ^+}(x,Q^2)\mathrm{\Delta }N_n^{\pi ^+}(x,Q^2)\}`$ (8)
$`\{J(Q^2)+1\}\{\mathrm{\Delta }N_p^\pi ^{}(x,Q^2)\mathrm{\Delta }N_n^\pi ^{}(x,Q^2)\}],`$
where $`J(Q^2)=\frac{3(D(Q^2)+\stackrel{~}{D}(Q^2))}{5(D(Q^2)\stackrel{~}{D}(Q^2))}`$. Eqs.(6) and (8) are main results of this work. Based on these formulas, one can extract $`\mathrm{\Delta }\overline{d}(x,Q^2)\mathrm{\Delta }\overline{u}(x,Q^2)`$ by using the values of $`\mathrm{\Delta }N_N^h(x,Q^2)`$ which can be derived from experimental data of spin asymmetries $`A_{1N}^h(x,Q^2)`$, if the spin–independent quark distribution functions and fragmentation functions are well known.

## 3 Extraction of $`\mathrm{\Delta }\overline{d}(x)\mathrm{\Delta }\overline{u}(x)`$ from semi-inclusive data

The remaining task is to numerically estimate the value of $`\mathrm{\Delta }\overline{d}(x,Q^2)\mathrm{\Delta }\overline{u}(x,Q^2)`$ from the present semi–inclusive data in order to examine how these formulas are effective for testing the light flavor symmetry of polarized distributions. In this analysis, we use the parametrization of GRV98(LO) for the unpolarized parton distribution being the $`\overline{u}/\overline{d}`$ asymmetric and the R1990 parametrization for the ratio $`R`$ in eq.(2). The fragmentation functions of eq.(4) are determined so as to fit well the EMC data and by integrating them from $`z_{min}=0.2`$ to $`1`$, we have obtained $`D^H(Q^2)`$ and $`\stackrel{~}{D}^H(Q^2)`$. At present, we have some data of $`A_{1p}^{h^\pm }`$ and $`A_{1d}^{h^\pm }`$ measured by the SMC group and also some data of $`A_{1p}^{h^\pm }`$, $`A_{1^3He}^{h^\pm }`$, $`A_{1p}^{\pi ^\pm }`$ and $`A_{1^3He}^{\pi ^\pm }`$ by the HERMES group. From these data, we can estimate the values of $`\mathrm{\Delta }\overline{d}(x,Q^2)\mathrm{\Delta }\overline{u}(x,Q^2)`$ from eqs.(6) and (8) by using $`\mathrm{\Delta }N_N^h`$ calculated from the data set of ($`A_{1p}^{h^\pm }`$, $`A_{1d}^{h^\pm }`$) by SMC and ($`A_{1p}^{h^\pm }`$, $`A_{1^3He}^{h^\pm }`$) and ($`A_{1p}^{\pi ^\pm }`$, $`A_{1^3He}^{\pi ^\pm }`$) by HERMES. Here, for the data of 3He targets, the values of $`A_{1n}^h`$ were derived from the data of $`A_{1^3He}^h`$ according to the way in ref.. In the present analysis, we have neglected the $`Q^2`$ dependence being fixed as $`Q_0^2=4`$GeV2 because no significant $`Q^2`$ dependence has been observed in this region in the spin asymmetry $`A_{1N}`$ for inclusive data. The results calculated from eqs.(6) and (8) are presented in fig.1. We have checked the model dependence of unpolarized quark distribution functions and found that the results are not sensitive to those models.

To examine the behavior of $`\mathrm{\Delta }\overline{d}(x)\mathrm{\Delta }\overline{u}(x)`$ in more detail and to test the light flavor asymmetry of $`\mathrm{\Delta }\overline{d}(x)`$ and $`\mathrm{\Delta }\overline{u}(x)`$, we have parametrized it as

$$\mathrm{\Delta }\overline{d}(x)\mathrm{\Delta }\overline{u}(x)=Cx^\alpha (\overline{d}(x)\overline{u}(x)),$$
(9)
and determined the values of $`C`$ and $`\alpha `$ from the $`\chi ^2`$–fit to the results presented in fig.1. The results were $`C=3.40(3.87)`$ and $`\alpha =0.567(0.525)`$ for the GRV98(LO)(MRST98(LO)) unpolarized distributions, while the values of $`\chi ^2`$/d.o.f. were 0.91(0.90) for GRV98(LO)(MRST98(LO)). $`C<0(0)`$ is a remarkable result, suggesting an asymmetry of $`\mathrm{\Delta }\overline{d}(x)`$ and $`\mathrm{\Delta }\overline{u}(x)`$. It is interesting to note that the negative value of $`C`$ is consistent with instanton interaction predictions. Also, the similar result is indicated from the chiral quark soliton model. However, it must be premature to lay stress on this result because of too large errors of the present data, though this result might suggest a violation of the polarized light flavor sea–quark symmetry. We urge to have more data with high precision to confirm this result.

Some comments are in order for the usefulness of our formulas: (i) Our formulas depend on the unpolarized parton distribution functions and the fragmentation functions. Unfortunately, some of them are poorly known at present. In addition, $`\mathrm{\Delta }N_{p(n)}^h`$ depends on the semi–inclusive asymmetry, $`A_{1p(n)}^h`$, and contains some experimental errors. Therefore, it might be rather difficult to extract the exact value of $`\mathrm{\Delta }\overline{d}(x)\mathrm{\Delta }\overline{u}(x)`$ from the present data. However, we believe that they must be quite useful for future experimental test of the polarized light sea-quark asymmetry if we have more precise data and good information on these functions. Our formulas are simple and can be easily tested in experiment. (ii) At present we see only asymmetries, $`A_{1p(n)}^h`$, in literature. However, if the precise experimental data on the polarized cross sections will be presented, then our formulas make more sense by replacing $`A_{1p(n)}^h`$ by the polarized cross sections themselves, where the unpolarized parton distributions and $`R(x,Q^2)`$ do not come in and we do not need to worry about their uncertainty.

## 4 Summary and discussion

In conclusion, we have proposed simple new formulas for extracting a difference of the polarized light sea–quark density, $`\mathrm{\Delta }\overline{d}(x)\mathrm{\Delta }\overline{u}(x)`$, from polarized deep–inelastic semi–inclusive data and numerically estimated it using these formulas from the present experimental data for semi–inclusive processes. Unfortunately, the precision of the present data is not enough for extracting an exact value of the difference, $`\mathrm{\Delta }\overline{d}(x)\mathrm{\Delta }\overline{u}(x)`$, and unambiguously testing the polarized light flavor sea–quark asymmetry. However, the HERMES group is now measuring semi–inclusive processes by using a new detector called RICH which can identify each of charged particles over a wide kinematical range and these data of charged pions with high statistics are expected to allow us to test more clearly the asymmetry of polarized light flavor sea–quark densities.

Another interesting way to study the asymmetry of $`\mathrm{\Delta }\overline{u}(x)`$ and $`\mathrm{\Delta }\overline{d}(x)`$ is the Drell–Yan process for polarized proton/deuteron–polarized proton collisions. The process provides informations on the raito of $`\mathrm{\Delta }\overline{u}(x)/\mathrm{\Delta }\overline{d}(x)`$ and thus it is complementary to our processes.

## Acknowledgements

One of us (T. Y.) thanks S. Kumano, H. Kitagawa and Y. Sakemi for valuable discussions.