text:
# Critical behavior in reaction-diffusion systems exhibiting absorbing phase transition

## I Introduction

The classification of universality classes of second order phase transitions is still one of the most important uncompleted task of statistical physics. Recently phase transitions of genuine nonequilibrium systems have intensively been investigated among reaction-diffusion (RD) type of models exhibiting absorbing states . There has been a hope that in such homogeneous systems symmetries and spatial dimensions are the most significant factors like in equilibrium ones, but gradually it turned out that there may be more relevant constituents as well. The most well known example is the parity conserving class (PC), which differs from the robust universality class of directed percolation (DP). The DP hypothesis stated by Janssen and Grassberger , according to which in one component systems exhibiting continuous phase transitions to single absorbing state (without extra symmetry and inhomogeneity or disorder) short ranged interactions can generate DP class transition only. However parity conservation itself proved to be an insufficient condition in many cases  and rather an underlying $`A3A`$, $`2A\mathrm{}`$ (BARW2) process  of particles and the $`Z_2`$ symmetry of absorbing states is necessary for this class . On the other hand parity conservation in N-component branching and annihilating random walk (N-BARW) systems , or by triplet production models  was found to be responsible for novel classes again. In one dimensional, multi-component reaction-diffusion systems site restriction turned out to be a relevant, new factor . Global conservation laws by directed percolation and lattice gas models were shown to be irrelevant , while systems with multiple absorbing states  or with multi-components also exhibit DP class scaling behavior .

An other important puzzle was being investigated intensively during the past three years that emerges at phase transitions of binary production reaction-diffusion systems  (PCPD). In these systems particle production by pairs competes with pair annihilation and single particle diffusion. If production wins steady states with finite particle density appear in (site restricted), while in unrestricted (bosonic) models the density diverges. By lowering the production/annihilation rate a doublet of absorbing states without symmetries emerges. One of such states is completely empty, the other possesses a single wandering particle. In case of site restricted systems the transition to absorbing state is continuous. It is important to note that these models do not break the DP hypothesis, because they exhibit multiple absorbing states which are not frozen, lonely particle(s) may diffuse in them. However no corresponding symmetry or conservation law has been found yet. Non-DP type of phase transition in a binary production system was already mentioned in the early work of Grassberger . A corresponding bosonic field theoretical model the annihilation fission (AF) process was introduced and studied by Howard and Täuber . These authors claim a non-DP type of transition in AF, because the action does not contain linear mass term and the theory is non-renormalizible perturbatively unlike the Reggeon field theory of DP. In field theories of models exhibiting DP class transition the canonical $`A\mathrm{}`$, $`A2A`$ reactions are generated by the renormalization transformation unlike here. Further facts opposing DP criticality are the set of different mean-field exponents and the different upper critical dimension of binary production PCPD like models ($`d_c=2`$ vs. $`d_c=4`$).

Forthcoming numerical studies reported somewhat different critical exponents, but there has been a consensus for about two years that this model should possess novel, non-DP type of transition. The first density matrix study by Carlon et al.  did not support a DP transition, but reported exponents near to those of the PC class. Since the PCPD does not conserve the particle number modulo 2, neither exhibits $`Z_2`$ symmetric absorbing states the PC criticality was unfavored. Simulation studies by Hinrichsen  and Ódor  and coherent anomaly calculations by Ódor  resulted in novel kind of critical behavior, although there was an uncertainty in the precise values of critical exponents. Exponent estimates showed diffusion ($`D`$) dependence that was under-pinned by pair mean-field results , possessing two distinct classes as the function of $`D`$. Recently Park et al. reported well defined set of critical exponents in different versions of binary production PCPD-like processes . However these simulations were done at a fixed, high diffusion/annihilation rate and agree with Ódor’s corresponding results . Kockelkoren and Chaté on the other hand claim an other set of critical exponents  that agrees with Ódor’s low diffusion/annihilation data.

The PCPD model can be mapped onto a two-component model  in which pairs are identified as a particle species following DP process and single particles as an other, coupled species following annihilating random walk. Simulations of such a two-component system at $`D=0.5`$ showed a continuous phase transition with exponents agreeing with those of the PCPD for high diffusions. This model is similar to an other one , which exhibits global particle number conservation as well. Field theory  and simulations  for the latter model reported two different universality classes as the function of $`D`$. It would be interesting to see if this conservation law is relevant or not like in case of the DP .

Interestingly, higher level cluster mean-field approximations result in a single class behavior by varying $`D`$ and it turned out that by assuming logarithmic corrections the single class scenario can be supported by simulations too . The origin of such logarithmic corrections may be a marginal perturbation between pairs and single particles in a coupled system description. A filed theoretical explanation would be necessary.

Two more recent studies  reported non-universality in the dynamical behavior of the PCPD. While in the former one Dickman and Menezes explored different sectors (a reactive and a diffusive one) in the time evolution and gave non-DP exponent estimates, in the latter one Hinrichsen set afloat a speculative conjecture that the ultimate long time behavior might be characterized by DP scaling behavior. In a forthcoming preprint  Hinrichsen provided a discussion about the possibility of the DP transition based on a series of assumptions. His starting point is a Langevin equation that is mapped onto a wetting process by Cole-Hopf transformation. By analyzing this process within a certain potential he gave arguments for a DP transition. While this Langevin equation with real noise is valid for the bosonic version of PCPD at and above the critical point, its usage in case of site occupancy restricted models is hypothetic, the noise can be complex at the transition point and may even change sign by the transformation. Furthermore the diffusive field of solitary particles is neglected.

In a very recent preprint  Barkema and Carlon continue this line and show that some simulation and density matrix renormalization results may also be interpreted as a signal of a phase transition belonging to the DP class. By assuming correction to scaling exponents that are equal to DP exponents and relevant up to quadratic or 3-rd order in the asymptotic limit they fitted their numerical results in case of two independent exponents. The extrapolations resulted is close to DP values for $`D=0.5`$. However for smaller $`D`$-s and by surface critical exponents this technique gave exponent estimates which are out of the error margin of DP.

An other novel class that may appear in triplet production systems was proposed in  (TCPD). This reaction-diffusion model differs from the PCPD that for a new particle generation at least three particles have to meet. For such generalizations Park et al. proposed a phenomenological Langevin equation that exhibits real, multiplicative noise . By simple power counting they found that the triplet model exhibits distinct mean-field exponents and upper critical dimension $`4/3d_c8/3`$. The simulations in 1d  indeed showed non-trivial critical exponents, which do not seem to correspond to any known universality classes. Kockelkoren and Chaté reported similar results in stochastic cellular automata (SCA) versions of general $`nA(n+k)A`$, $`mA(ml)A`$ type of models , where multiple particle creation on a given site is suppressed by an exponentially decreasing creation probability ($`p^{N/2}`$) of the particle number. They claim that their simulation results in 1d are in agreement with the fully occupation number restriction counterparts and set up a general table of universality classes, where as the function of $`n`$ and $`m`$ only 4 non-mean-field classes exist, namely the DP class, the PC class, the PCPD and TCPD classes. However more extensive simulations of 1 and 2 dimensional site restricted lattice models  do not support some of these results in case of different triplet and quadruplet models. In the $`3A4A`$ $`3A\mathrm{}`$ triplet model 1d numerical data can be interpreted as mean-field behavior with logarithmic corrections and in two dimensions clear mean-field exponents appear, hence the upper critical dimension is $`d_c=1`$, which contradicts the Langevin equation prediction. Surprisingly other non-trivial critical behavior were also detected in the $`3A5A`$ $`3AA`$ parity conserving triplet model and in some quadruplet models . The cause of differences between the results of these studies is subject of further investigations. Again proper field theoretical treatment would be important.

The classification of universality classes of nonequilibrium systems by the exponent $`\mu `$ of a multiplicative noise in the Langevin equation was suggested some time ago by Grinstein et al. . However it turned out that there may not be corresponding particle systems to real multiplicative noise cases  and an imaginary part appears as well if one derives the Langevin equation of a RD system starting from the master equation in a proper way. Furthermore for higher-order processes the emerging nonlinearities in the master equation action do not allow a rewriting in terms of Langevin-type stochastic equations of motion, hence for high-order processes like those of the TCPD a Langevin representation may not exist.

This situation resembles to some extent to a decade long debate over the critical phase transition of driven diffusive systems . The latest papers in this topic suggest that the phenomenological Langevin equation originally set up for such systems do not correspond exactly to the lattice models investigated. Simulations for different lattice models show, that instead of an external current the anisotropy is the real cause of the critical behavior observed in simulations .

## II Mean-field classes

In this section I show that mean-field classes of site restricted lattice models with general microscopic processes of the form

$$nA\stackrel{\sigma }{}(n+k)A,mA\stackrel{\lambda }{}(ml)A,$$
(1)
with $`n>1`$, $`m>1`$, $`k>0`$, $`l>0`$ and $`ml0`$ are different from those of the DP and PC processes backing numerical results which claim novel type of criticality below $`d_c`$. The mean-field equation that can be set up for the lattice version of these processes (with creation probability $`\sigma `$ and annihilation or coagulation probability $`\lambda =1\sigma `$) is

$$\frac{\rho }{t}=ak\sigma \rho ^n(1\rho )^kal(1\sigma )\rho ^m,$$
(2)
where $`\rho `$ denotes the site occupancy probability and $`a`$ is a dimension dependent coordination number. Each empty site has a probability (1-$`\rho `$) in mean-field approximation, hence the need for $`k`$ empty sites at a creation brings in a $`(1\rho )^k`$ probability factor. By expanding $`(1\rho )^k`$ and keeping the lowest order contribution one can see that for site restricted lattice systems a $`\rho ^{n+1}`$-th order term appears with negative coefficient that regulates eq.(2) in the active phase. In the inactive phase one expects a dynamical behavior dominated by the $`mA\mathrm{}`$ process, for which the particle density decay law is known $`\rho (t)t^{1/(m1)}`$ . The steady state solutions were determined in  analytically and one can distinguish three different situations at the phase transition: (a) $`n=m`$, (b) $`n>m`$ and (c) $`n<m`$ .

### A The $`n=m`$ symmetric case

As discussed in  the leading order singularities of steady state solution can be obtained. By approaching the the critical point $`\sigma _c=\frac{l}{k+l}`$ in the active phase the steady state density vanishes continuously as

$$\rho |\sigma \sigma _c|^{\beta ^{MF}},$$
(3)
with the order parameter exponent exponent $`\beta ^{MF}=1`$. At the critical point the density decays with a power-law

$$\rho t^{\alpha _{MF}},$$
(4)
with $`\alpha _{MF}=\beta ^{MF}/\nu _{||}^{MF}=1/n`$, hence $`\nu _{||}^{MF}=n`$, providing a series of mean-field universality classes for $`n>1`$ (besides DP an PC where $`\nu _{||}^{MF}=1`$) and backing the results, which claim novel type of non-trivial transitions below the critical dimension. Unfortunately determining the the value of $`d_c`$ is a non-trivial task without a proper Langevin equation. These scaling exponents can be obtained from bosonic, coarse grained formulation too , where a $`\rho ^{n+1}`$-th order term, with negative coefficient had to be added by hand to suppress multiple site occupancy. It is known however that hard-core particle blocking may result in relevant perturbation in $`d=1`$ dimension , so for cases where the upper critical dimension is $`d_c1`$ the site restricted, $`N>1`$ cluster mean-field approximation that takes into account diffusion would be a more adequate description of the model (see ).

### B The $`n>m`$ case

In this case the mean-field solution provides first order transition (see ), hence it does not imply anything with respect to possible classes for models below the critical dimension ($`d<d_c`$). Note however, that by higher order cluster mean-field approximations, where the diffusion plays a role the transition may turn into continuous one (see for example ). The simulations by Kockelkoren and Chaté report DP class transition for such models in one dimension .

### C The $`n<m`$ case

In this case the critical point is at zero branching rate $`\sigma _c=0`$, where the density decays with $`\alpha ^{MF}=1/(m1)`$ as in case of the $`n=1`$ branching and $`m=l`$ annihilating models showed by Cardy and Täuber  (BkARW classes). However the steady state solution for particle production with $`n>1`$ parents gives different $`\beta `$ exponents than those of BkARW classes, namely $`\beta ^{MF}=1/(mn)`$, defining a whole new series of mean-field classes , for a simplicity I shall call them PkARW classes. It is important to note, that one has not found a corresponding symmetry or conservation law to these classes. These mean-field classes imply novel kind of critical behavior for $`d<d_c`$. For $`n=2`$, $`m=3`$ the mean-field exponents are $`\beta =1`$ and $`\alpha =1/2`$ agreeing with the mean-field exponents of the PCPD class. Indeed for $`n=2`$, $`m=3`$  reports PCPD class dynamical criticality. This supports the expectation that non-mean-field classes follow the distinctions observed in the corresponding mean-field classes. To go further in sections III and IV I investigate the phase transition of the simplest unexplored PkARW classes, in the $`n=2`$, $`m=4`$ model.

### D The role of $`k`$ and $`l`$

In the mean-field approximation $`k`$ and $`l`$ do not affect the universal properties, however simulations in one dimension showed  that in case of the $`m=n=k=l=3`$ model the critical point was shifted to zero branching rate and a BkARW class transition emerged there contrary to what was expected for $`n=m`$. For a stochastic cellular automaton version of these reactions  reported a non-trivial critical transition. In general one may expect such effects for large $`k`$ and $`l`$ values, for which N-cluster mean-field approximation – that takes into account diffusion – would give a better description.

## III Simulations of the $`2A4A`$, $`4A2A`$ model in two dimensions

In the II B section I introduced PkARW mean-field classes for $`n<m`$. Here I explore the phase transition in the simplest model from this class, in the $`2A4A`$, $`4A2A`$ model with $`D=0.5`$ diffusion rate. Two dimensional simulations were performed on $`L=1000`$ linear sized lattices with periodic boundary conditions. One Monte Carlo step (MCS) — corresponding to $`dt=1/P`$ (where $`P`$ is the number of particles) — is built up from the following processes. A particle and a number $`x(0,1)`$ are selected randomly; if $`x<D=0.5`$ a site exchange is attempted with one of the randomly selected empty nearest neighbors (nn); if $`xD=0.5`$ two particles are created with probability $`\sigma `$ at randomly selected empty nn sites provided the number of nn particles was greater than or equal $`2`$; or if $`x0.5`$ two particles are removed with probability $`1\sigma `$. The simulations were started from fully occupied lattices and the particle density ($`\rho (t)`$) decay was followed up to $`4\times 10^5`$ MCS.

First the critical point was located by measuring the dynamic behavior of $`\rho (t)`$. It turned out that the transition is at zero branching rate ($`\sigma _c=0`$). The density decay was analyzed by the local slopes defined as

$$\alpha _{eff}(t)=\frac{\mathrm{ln}\left[\rho (t)/\rho (t/m)\right]}{\mathrm{ln}(m)}$$
(5)
where I used $`m=4`$. As Fig.1 shows the local slopes curve for $`t>10^5`$ MCS extrapolates to the mean-field value $`\alpha =0.334(1)`$. This value agrees with the mean-field value $`\alpha ^{MF}=1/3`$.

Density decays for several $`\sigma `$-s in the active phase ($`0.0002\sigma 0.05`$) were followed on logarithmic time scales and averaging was done over $`100`$ independent runs in a time window, which exceeds the level-off time by a decade. The steady state density in the active phase near the critical phase transition point is expected to scale as

$$\rho (\mathrm{},\sigma )|\sigma \sigma _c|^\beta .$$
(6)
Using the local slopes method one can get a precise estimate for $`\beta `$ as well as for the corrections to scaling

$$\beta _{eff}(\sigma _i)=\frac{\mathrm{ln}\rho (\mathrm{},\sigma _i)\mathrm{ln}\rho (\mathrm{},\sigma _{i1})}{\mathrm{ln}(\sigma _i)\mathrm{ln}(\sigma _{i1})}.$$
(7)
As one can see on Fig.2 the effective exponent clearly tends to the expected mean-field value $`\beta =0.5`$ as $`\sigma 0`$. Assuming a correction to scaling of the form

$$\beta _{eff}=\beta at^\delta $$
(8)
non-linear fitting results in $`\delta =0.44(1)`$ correction to scaling exponent.

Besides these scaling correction assumptions I also tried to apply different, lowest order logarithmic corrections to the data, but these fittings gave exponents slightly away from mean-field values and the corresponding coefficients proved to be very small, therefore I concluded that $`d_c<2`$. In the next section I do the same analysis in $`d=1`$.

## IV Simulations of the $`2A4A`$, $`4A2A`$ model in one dimensions

The simulations in one dimension were carried out on $`L=20000`$ sized systems with periodic boundary conditions. The initial states were again fully occupied lattices, and the density of particles is followed up to $`4\times 10^6`$ MCS. An elementary MCS consists of the following processes:

$`A\mathrm{}\mathrm{}A`$ with probability D,

$`AAAA\mathrm{}AA\mathrm{}`$ with probability $`(1\sigma )(1D)`$,

$`AA\mathrm{}\mathrm{}AAAA`$ or $`\mathrm{}\mathrm{}AAAAAA`$ with probability $`\sigma (1D)`$,

The critical point was located at $`\sigma _c=0`$ again. As one can see on Figure 3 there is a crossover of the local slopes for $`t>5\times 10^5`$ MCS and a linear extrapolation for this region results in $`\alpha =0.329(5)`$ agreeing with the mean-field value.

The steady state data were analyzed in the active region for $`0.0003\sigma 0.5`$ as in two dimensions, by the local slopes method (eq. 7) and by assuming correction to scaling of the form (8). This resulted in $`\delta ^{}=0.332`$ correction to scaling exponent. The local slopes plotted as the function $`\sigma ^\delta `$ shown on Fig. 4 extrapolates to $`\beta =0.49(1)`$ in agreement with the mean-field exponent.

## V Conclusions

In this paper I reviewed and discussed the state of the art of the phase transitions of reaction-diffusion systems exhibiting explicit diffusion and production by $`n>1`$ parents. Arguments are given against DP criticality that has recently been suggested in some papers. These are supported by a series of mean-field classes that can be classified by the existence of a $`n=m`$ symmetry in the system and by the $`n`$ and $`m`$ values. Especially the need for a proper field theoretical treatment is emphasized. The upper critical dimension in these models is not known.

I determined by simulations the $`\alpha `$ and $`\beta `$ exponents of the $`2A4A`$, $`4A2A`$ model in one an two dimensions. These results indicate that for this binary production system the upper critical dimension is $`d_c<1`$. This model conserves the parity of particles still its transition does not belong to the PC class. Note that while the density decay results in one dimension are in agreement with that of Kockelkoren and Chaté-s simulations , the off-critical order parameter exponent is $`\beta ^{MF}=1/(mn)`$ which shows that there are more classes exist at zero branching rate besides the BkARW universality classes. It is still an open question if there is any variant of PkARW models that exhibits non-mean-field criticality in physical dimensions ($`d1`$).

Acknowledgements:

Support from Hungarian research funds OTKA (Grant No. T-25286) is acknowledged. The author thanks the trial access to the NIIFI Cluster-GRID of Hungary.


format: text

text:
# References

Over the decades matrix models have contributed to our understanding of gauge fields and strings and their synthesis in different contexts. The recent proposal of Dijkgraaf and Vafa , (see  for prehistory to this), has enabled us to relate exact results on the superpotential and gauge couplings for a wide class of $`𝒩=1`$ supersymmetric gauge theories with the planar limit of a matrix model in which the superpotential is taken as an ordinary potential. One may regard this relation as that connecting two faces of the notion of integrability of these theories in their infrared limit. The one face is seen as the integrability of Seiberg-Witten, $`i.e.`$, $`𝒩=2`$ theories   with soft breaking superpotential as Whitham flows . Here the aspects associated with gaugino condensates are not revealed immediately. The other face is, of course, the integrability of matrix models where the (planar) free energy provides an integral representation of the prepotential. Here, bosonic integrations of matrices have no way to capture the properties associated with supersymmetry and we have to input a given breaking pattern of the original gauge theory. It seems that we need to regard both gaugino condensates and vacuum values as moduli in order to gain unified understanding of this subject . Studies of the Dijkgraaf-Vafa proposal from these perspectives have been continuing both on its structure and explicit calculations . For the other aspects on this subject, see an extensive list of references seen, for instance, in  and the last reference of .

The proposal of Dijkgraaf and Vafa has been put and consolidated in a field theoretic vein by Cachazo, Douglas, Seiberg and Witten  in the case of an adjoint chiral superfield and an arbitrary superpotential. In particular, they have derived an equation of the anomalous Ward-Takahashi identity associated with the Konishi anomaly  which takes the same form as the loop equation of the ordinary hermitian one-matrix model in the planar limit. In their discussion, there is actually a system of three Ward-Takahashi identities involving generating functions $`R(z),w_\alpha (z),T(z)`$, which are all chiral operators related by constant fermionic shift of the gaugino field strengths.111 This fermionic shift is allowed by the decoupling of the $`U(1)`$ factor in the adjoint representation. The first operator $`R(z)`$ realizes the proposal of  while the third operator $`T(z)`$ is closely related to the Seiberg-Witten differential. 222This latter statement holds even in those cases in which unbroken nonabelian gauge symmetry survive. The system of three Ward-Takahashi identities appears to offer a unified understanding of  and $`𝒩=2`$ Seiberg-Witten theory. See  for a recent progress. In the ordinary bosonic one-matrix model, on the other hand, there is only one kind of resolvent whose Laplace transform is a loop operator, and is no hope within this model to obtain the second, which is fermionic, and the third resolvents carrying lower dimensions than the first one. A model which derives all of the three equations of  and which provides a constructive definition serving both for the planar processes and for the non-planar (and therefore) gravitational ones has remained elusive.

To incorporate the counterpart of the fermionic shift operation, fermionic degrees of freedom are required. The supereigenvalue model proposed sometime ago in , which materializes discretized two-dimensional supergravity, (see  for some of the subsequent developments), has come to our attention.333Supermatrices have been known not to serve well in the old context of $`2d`$ supergravity. (See, for instance,  for this). For discussions in the present context, see . There is a subtle issue of a decoupling of the superunitary group $`U(M|N)`$ in the integrations over the supermatrix. These integrations in general do not reduce to those of the (super)eigenvalues. Furthermore, supertrace requires us to introduce oppositely charged particles in the Dyson gas description of the (super)eigenvalues. The supermatrices appear to us more appropriate to the quiver gauge theories and not to the pure Yang-Mills.  In this letter, we will present a variant of the supereigenvalue model in which the ordinary grassmann $`𝒩=2`$ superspace and its supercoordinates $`(\lambda ,\theta ,\overline{\theta })`$ 444 $`\theta `$ and $`\overline{\theta }`$ are two independent grassmann variables and no complex conjugation is involved. are introduced to label a set of eigenvalues $`\lambda _i`$ and their grassmann partners $`\theta _i,\overline{\theta }_i`$. We will show that this model captures the three equations of  as planar loop equations.

Let us first consider the following partition function:

$`Z_{\stackrel{ˇ}{N}}={\displaystyle \mathrm{}\underset{i=1}{\overset{\stackrel{ˇ}{N}}{}}d\lambda _id\theta _id\overline{\theta }_i\underset{i<j}{}(\lambda _i\lambda _j+\theta _i\overline{\theta }_i\theta _j\overline{\theta }_j)^2\mathrm{exp}\left(\frac{1}{g_s}\underset{i}{}𝐕(\lambda _i,\theta _i,\overline{\theta }_i)\right)},`$ (1)
where $`g_s\frac{S}{\stackrel{ˇ}{N}}`$ and

$`𝐕(\lambda _i,\theta _i,\overline{\theta }_i)=W(\lambda _i;g_{\mathrm{}})+\overline{W}^{(s)}(\lambda _i;\overline{\xi }_{\mathrm{}})\theta _i+\overline{\theta }_iW^{(s)}(\lambda _i;\xi _{\mathrm{}})+\stackrel{~}{W}(\lambda _i;\stackrel{~}{g}_{\mathrm{}})\theta _i\overline{\theta }_i`$ (2)
is an $`𝒩=2`$ superfield in the terminology of the two-dimensional superconformal field theory in which the coefficients of the Taylor expansion in $`\lambda _i`$ of the components $`W,\overline{W}^{(s)},W^{(s)},\stackrel{~}{W}`$ are realized as the couplings $`g_{\mathrm{}},\overline{\xi }_{\mathrm{}},\xi _{\mathrm{}},\stackrel{~}{g}_{\mathrm{}}`$. All of the singlet observables are obtained from this partition function by acting upon these couplings. We have fixed the normalization of the grassmann variable, so that the relative coefficient between $`\lambda _i\lambda _j`$ and $`\theta _i\overline{\theta }_i\theta _j\overline{\theta }_j`$ in the measure factor is $`1`$. 555 In this letter, we set the value of the exponent of the measure factor to $`2`$. For a general value $`\beta `$, the coefficient of an order $`g_s`$ term in the loop equation below gets modified after a rescaling of the superpotential and the resolvent.

We bring this model into more direct contact to the proposal of  and its description as gauge theory in  by reducing this partition function. Let us set

$$\overline{W}^{(s)}=W^{(s)}=0,\stackrel{~}{W}=W^{},$$
(3)
which is, of course, a nonsupersymmetric reduction. We regard the function $`W`$ an $`n`$-th order polynomial in its argument. The model obtained this way is

$`Z_{\stackrel{ˇ}{N}}={\displaystyle \mathrm{}\underset{i=1}{\overset{\stackrel{ˇ}{N}}{}}d\lambda _id\theta _id\overline{\theta }_i\underset{i<j}{}(\lambda _i\lambda _j+\theta _i\overline{\theta }_i\theta _j\overline{\theta }_j)^2\mathrm{exp}\left(\frac{1}{g_s}\underset{i}{}W(\lambda _i+\theta _i\overline{\theta }_i)\right)}.`$ (4)
There is another supereigenvalue model which possesses $`𝒩=2`$ supersymmetry, which we will present briefly in the end of this letter. In contract with the model (4), the model with $`𝒩=2`$ supersymmetry is obtained by a different choice of the measure factor $`_{i<j}(\lambda _i\lambda _j(\theta _i\overline{\theta }_j+\overline{\theta }_i\theta _j))^2`$. 666These two alternative choices may be characterized as follows: by requiring the antisymmetry in $`(i,j)`$ and $`𝒩=2`$ neutralness, possible bilinear combinations are either $`\theta _i\overline{\theta }_i\theta _j\overline{\theta }_j`$ or $`\theta _i\overline{\theta }_j+\overline{\theta }_i\theta _j`$. The difference is that the former is odd under the exchange of $`\theta `$ and $`\overline{\theta }`$ while the latter is even. Thus the former may be called chiral choice and the latter non-chiral. In fact the combinations $`\lambda \pm \theta \overline{\theta }`$ are the natural coordinates for chiral and anti-chiral superfields respectively.

Let us briefly recall the logic which connects the Virasoro constraints at finite $`\stackrel{ˇ}{N}`$ with its loop equation in the ordinary hermitian one-matrix model. Consider the partition function with all grassmann degrees of freedom removed in eq.(1). Setting up the system of Schwinger-Dyson equations associated with most general (polynomial) variations of the matrix is equivalent to inserting the Virasoro operators (which is a total derivative) $`\mathrm{}_m^{(B)}{\displaystyle \underset{i=1}{\overset{\stackrel{ˇ}{N}}{}}}\frac{d}{d\lambda _i}\lambda _i^{m+1},m1`$ in the partition function and equating the resulting expressions to zero. Summing over $`m`$ by introducing an expansion parameter $`x`$, we obtain

$`0`$ $`=`$ $`{\displaystyle \mathrm{}\underset{i}{}d\lambda _i\left(\underset{j}{}\frac{d}{d\lambda _j}\frac{1}{x\lambda _j}\right)\mathrm{exp}\left(S_{\mathrm{eff}}^{(B)}\right)},`$
$`S_{\mathrm{eff}}^{(B)}`$ $``$ $`2{\displaystyle \underset{i<j}{}}\mathrm{log}(\lambda _i\lambda _j)+{\displaystyle \frac{1}{g_s}}{\displaystyle \underset{i}{}}W\left(\lambda _i\right).`$
Carrying out the differentiations and straightforward algebras, we obtain the finite $`\stackrel{ˇ}{N}`$ loop equation for the resolvent operator $`\omega _{\stackrel{ˇ}{N}}^{(B)}(x)g_s{\displaystyle \underset{i=1}{\overset{\stackrel{ˇ}{N}}{}}}\frac{1}{x\lambda _i}`$ :

$`0=\omega _{\stackrel{ˇ}{N}}^{(B)}(x)^2+W^{}(x)\omega _{\stackrel{ˇ}{N}}^{(B)}(x)+{\displaystyle \frac{S}{\stackrel{ˇ}{N}}}{\displaystyle \underset{i}{}}{\displaystyle \frac{W^{}(x)W^{}(\lambda _i)}{x\lambda _i}}.`$
Here $`\mathrm{}`$ denotes the averaging with respect to the partition function. The absence of the $`g_s\omega _{\stackrel{ˇ}{N}}^{(B)}(x)`$ term is related to $`c=1`$ nature of the Virasoro constraints. In the planar limit in which the correlators of the singlet operators factorize, we obtain $`\omega ^{(B)}(x)\underset{\stackrel{ˇ}{N}\mathrm{}}{lim}\omega _{\stackrel{ˇ}{N}}^{(B)}(x)`$:

$`\omega ^{(B)}(x)^2+W^{}(x)\omega ^{(B)}(x)+\underset{\stackrel{ˇ}{N}\mathrm{}}{lim}{\displaystyle \frac{S}{\stackrel{ˇ}{N}}}{\displaystyle \underset{i}{}}{\displaystyle \frac{W^{}(x)W^{}(\lambda _i)}{x\lambda _i}}=0.`$
We now apply the above procedure to our model given by eq. (4). We find that the appropriate realization of the Virasoro operators which involve the grassmann eigenvalues as well is

$`\mathrm{}_m`$ $``$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{d}{d\lambda _i}}\lambda _i^{m+1}+\gamma (m+1){\displaystyle \underset{i}{}}\left({\displaystyle \frac{d}{d\theta _i}}\theta _i+{\displaystyle \frac{d}{d\overline{\theta }_i}}\overline{\theta }_i\right)\lambda _i^m,m1`$
$`\widehat{\mathrm{}}(x)`$ $`=`$ $`{\displaystyle \underset{j}{}}{\displaystyle \frac{d}{d\lambda _j}}{\displaystyle \frac{1}{x\lambda _j}}+\gamma {\displaystyle \underset{j}{}}\left({\displaystyle \frac{d}{d\theta _j}}\theta _j+{\displaystyle \frac{d}{d\overline{\theta }_j}}\overline{\theta }_j\right){\displaystyle \frac{1}{(x\lambda _j)^2}}.`$ (5)
This form can be motivated by realizing the $`𝒩=2`$ superconformal algebra in terms of (super-)differential operators. The coefficient $`\gamma `$ is then determined:

$$\gamma =\frac{1}{2}.$$
(6)
Inserting $`\widehat{\mathrm{}}(x)`$, which is a total derivative, in $`Z_{\stackrel{ˇ}{N}}`$, carrying out algebras, and equating the resulting expression to zero, we obtain

$`0`$ $`=(12\gamma ){\displaystyle \underset{j}{}}{\displaystyle \frac{1}{(x\lambda _j)^2}}`$
$``$ $`{\displaystyle \frac{1}{g_s}}{\displaystyle \underset{j}{}}\left({\displaystyle \frac{1}{x\lambda _j}}{\displaystyle \frac{d}{d\lambda _j}}+\gamma {\displaystyle \frac{1}{(x\lambda _j)^2}}\left(\theta _j{\displaystyle \frac{d}{d\theta _j}}+\overline{\theta }_j{\displaystyle \frac{d}{d\overline{\theta }_j}}\right)\right)W(\lambda _j+\theta _j\overline{\theta }_j)`$
$`+`$ $`2{\displaystyle \underset{j}{}}\left({\displaystyle \frac{1}{x\lambda _j}}{\displaystyle \frac{d}{d\lambda _j}}+\gamma {\displaystyle \frac{1}{(x\lambda _j)^2}}\left(\theta _j{\displaystyle \frac{d}{d\theta _j}}+\overline{\theta }_j{\displaystyle \frac{d}{d\overline{\theta }_j}}\right)\right){\displaystyle \underset{i<k}{}}\mathrm{log}\left(\lambda _i\lambda _k+\theta _i\overline{\theta }_i\theta _k\overline{\theta }_k\right).`$
The first line of this equation vanishes with our choice $`\gamma =\frac{1}{2}`$. The second line with this choice is

$`{\displaystyle \frac{1}{g_s}}{\displaystyle \underset{j}{}}{\displaystyle \frac{W^{}\left(\lambda _j+\theta _j\overline{\theta }_j\right)}{x\lambda _j\theta _j\overline{\theta }_j}}.`$
The calculation of the third line is a little less trivial to see. After antisymmetrization and isolation of purely bosonic term, we obtain

$`{\displaystyle \underset{ji}{}}{\displaystyle \frac{1}{(x\lambda _j)(x\lambda _i)}}{\displaystyle \underset{ji}{}}{\displaystyle \frac{\theta _j\overline{\theta }_j\theta _i\overline{\theta }_i}{(x\lambda _j)(x\lambda _i)(\lambda _j\lambda _i+\theta _j\overline{\theta }_j\theta _i\overline{\theta }_i)}}`$
$`+2\gamma {\displaystyle \underset{ji}{}}\left({\displaystyle \frac{\theta _j\overline{\theta }_j}{(x\lambda _j)^2}}{\displaystyle \frac{\theta _i\overline{\theta }_i}{(x\lambda _i)^2}}\right){\displaystyle \frac{1}{\lambda _j\lambda _i+\theta _j\overline{\theta }_j\theta _i\overline{\theta }_i}}.`$
The cancellation of the $`\frac{1}{x^2}`$ contribution from the second and the third terms of this equation in the asymptotic expansion again requires $`\gamma =\frac{1}{2}`$. With this choice, the third line of eq.(S0.Ex7) is

$$\underset{ij}{}\frac{1}{(x\lambda _j\theta _j\overline{\theta }_j)(x\lambda _i\theta _i\overline{\theta }_i)}.$$
(8)
Introducing the resolvent operator $`\omega _{\stackrel{ˇ}{N}}(x)g_s{\displaystyle \underset{i=1}{\overset{\stackrel{ˇ}{N}}{}}}\frac{1}{x\lambda _i\theta _i\overline{\theta }_i}`$, we obtain

$`0`$ $`=`$ $`\omega _{\stackrel{ˇ}{N}}(x)^2g_s\omega _{\stackrel{ˇ}{N}}^{}(x)+W^{}(x)\omega _{\stackrel{ˇ}{N}}(x)+{\displaystyle \frac{S}{\stackrel{ˇ}{N}}}{\displaystyle \underset{i}{}}{\displaystyle \frac{W^{}(x)W^{}(\lambda _i+\theta _i\overline{\theta }_i)}{x\lambda _i\theta _i\overline{\theta }_i}}`$ (9)
$``$ $`_1(x).`$
The fourth term is a polynomial in $`x`$. In the planar limit $`\omega (x)lim\omega _{\stackrel{ˇ}{N}}(x)`$, this equation reduces to

$$\omega (x)^2+W^{}(x)\omega (x)+lim\frac{S}{\stackrel{ˇ}{N}}\underset{i}{}\frac{W^{}(x)W^{}(\lambda _i+\theta _i\overline{\theta }_i)}{x\lambda _i\theta _i\overline{\theta }_i}=0.$$
(10)
So far we have managed to obtain a single loop equation for the resolvent. The observables of our model are, on the other hand,

$`{\displaystyle \underset{i}{}}\lambda _i^{\mathrm{}},{\displaystyle \underset{i}{}}\lambda _i^{\mathrm{}}\theta _i,{\displaystyle \underset{i}{}}\lambda _i^{\mathrm{}}\overline{\theta }_i,{\displaystyle \underset{i}{}}\lambda _i^{\mathrm{}}\overline{\theta }_i\theta _i,\mathrm{}0.`$
With the single expansion parameter $`x`$, we are not able to generate all of these which are independent. In order to overcome this point, we insert

$`{\displaystyle \frac{1}{2}}\left(1+\eta {\displaystyle \underset{i}{}}{\displaystyle \frac{d}{d\theta _i}}+{\displaystyle \underset{i}{}}{\displaystyle \frac{d}{d\overline{\theta }_i}}\overline{\eta }\right)^2,`$
following the insertion of the Virasoro operators in eq.(S0.Ex7) and carry out the remaining procedure. The grassmann expansion parameters $`\eta ,\overline{\eta }`$ have now been introduced. Equivalently, we can make a fermionic shift of the grassmann eigenvalues in eq.(9):

$`\theta _i\theta _i+\eta ,\overline{\theta }_i\overline{\theta }_i+\overline{\eta }.`$
Both the resolvent operator $`\omega _{\stackrel{ˇ}{N}}(x)`$ and

$$S_{\mathrm{eff}}\frac{1}{g_s}\underset{i}{}W(\lambda _i+\theta _i\overline{\theta }_i)2\underset{i<j}{}\mathrm{log}(\lambda _i\lambda _j+\theta _i\overline{\theta }_i\theta _j\overline{\theta }_j)$$

undergo this shift while $`d(\theta _i+\eta )=d\theta _i,d(\overline{\theta }_i+\overline{\eta })=d\overline{\theta }_i`$. In either way, we obtain another set of Schwinger-Dyson equations as equations multiplying $`\eta ,\overline{\eta },\eta \overline{\eta }`$;

$`0`$ $`=`$ $`_2(x)+𝒮_\theta ^{\mathrm{sub}}_1(x),`$
$`0`$ $`=`$ $`_{\overline{2}}(x)+𝒮_{\overline{\theta }}^{\mathrm{sub}}_1(x),`$ (11)
$`0`$ $`=`$ $`_3(x)+\left(𝒮_{\theta \overline{\theta }}^{\mathrm{sub}}+𝒮_\theta ^{\mathrm{sub}}𝒮_{\overline{\theta }}^{\mathrm{sub}}𝒮_\theta ^{\mathrm{sub}}𝒮_{\overline{\theta }}^{\mathrm{sub}}\right)_1(x)+_{\overline{2}}(x)𝒮_{\overline{\theta }}^{\mathrm{sub}}+𝒮_\theta ^{\mathrm{sub}}_2(x),`$
where

$`_2(x)`$ $``$ $`2\psi _{\stackrel{ˇ}{N}}(x)\omega _{\stackrel{ˇ}{N}}(x)g_s\psi _{\stackrel{ˇ}{N}}^{}(x)+W^{}(x)\psi _{\stackrel{ˇ}{N}}(x)+\rho _{\stackrel{ˇ}{N}}(x),`$
$`_{\overline{2}}(x)`$ $``$ $`2\overline{\psi }_{\stackrel{ˇ}{N}}(x)\omega _{\stackrel{ˇ}{N}}(x)g_s\overline{\psi }_{\stackrel{ˇ}{N}}^{}(x)+W^{}(x)\overline{\psi }_{\stackrel{ˇ}{N}}(x)+\overline{\rho }_{\stackrel{ˇ}{N}}(x),`$ (12)
$`_3(x)`$ $``$ $`2t_{\stackrel{ˇ}{N}}(x)\omega _{\stackrel{ˇ}{N}}(x)+2\overline{\psi }_{\stackrel{ˇ}{N}}(x)\psi _{\stackrel{ˇ}{N}}(x)g_st_{\stackrel{ˇ}{N}}^{}(x)+W^{}(x)t_{\stackrel{ˇ}{N}}(x)+c_{\stackrel{ˇ}{N}}(x),`$
and

$`\rho _{\stackrel{ˇ}{N}}(x)`$ $`=`$ $`{\displaystyle \frac{S}{\stackrel{ˇ}{N}}}{\displaystyle \underset{i}{}}\theta _i{\displaystyle \frac{W_{\mathrm{sub}}^{(2)}(x;\lambda _i+\theta _i\overline{\theta }_i)}{\left(x\lambda _i\theta _i\overline{\theta }_i\right)^2}},\overline{\rho }_{\stackrel{ˇ}{N}}(x)={\displaystyle \frac{S}{\stackrel{ˇ}{N}}}{\displaystyle \underset{i}{}}\overline{\theta }_i{\displaystyle \frac{W_{\mathrm{sub}}^{(2)}(x;\lambda _i+\theta _i\overline{\theta }_i)}{\left(x\lambda _i\theta _i\overline{\theta }_i\right)^2}},`$
$`c_{\stackrel{ˇ}{N}}(x)`$ $`=`$ $`{\displaystyle \frac{S}{\stackrel{ˇ}{N}}}{\displaystyle \underset{i}{}}{\displaystyle \frac{W_{\mathrm{sub}}^{(2)}(x;\lambda _i+\theta _i\overline{\theta }_i)}{\left(x\lambda _i\theta _i\overline{\theta }_i\right)^2}}+{\displaystyle \frac{S}{\stackrel{ˇ}{N}}}{\displaystyle \underset{i}{}}2\overline{\theta }_i\theta _i{\displaystyle \frac{W_{\mathrm{sub}}^{(3)}(x;\lambda _i+\theta _i\overline{\theta }_i)}{\left(x\lambda _i\theta _i\overline{\theta }_i\right)^3}}.`$ (13)
We have introduced

$$W_{\mathrm{sub}}^{(m+1)}(x;\lambda _i+\theta _i\overline{\theta }_i)W^{}(x)\underset{\mathrm{}=0}{\overset{m}{}}\frac{(x\lambda _i\theta _i\overline{\theta }_i)^{\mathrm{}}}{\mathrm{}!}W^{(\mathrm{}+1)}(\lambda _i+\theta _i\overline{\theta }_i).$$

We see that the quantities seen in eq.(S0.Ex20) are polynomials in $`x`$. Another bosonic resolvent and two fermionic resolvents are defined respectively by

$`t_{\stackrel{ˇ}{N}}(x)`$ $``$ $`g_s{\displaystyle \underset{i}{}}\left({\displaystyle \frac{1}{\left(x\lambda _i\theta _i\overline{\theta }_i\right)^2}}+{\displaystyle \frac{2\overline{\theta }_i\theta _i}{\left(x\lambda _i\theta _i\overline{\theta }_i\right)^3}}\right),`$
$`\psi _{\stackrel{ˇ}{N}}(x)`$ $``$ $`g_s{\displaystyle \underset{i}{}}{\displaystyle \frac{\theta _i}{\left(x\lambda _i\theta _i\overline{\theta }_i\right)^2}},\overline{\psi }_{\stackrel{ˇ}{N}}(x)g_s{\displaystyle \underset{i}{}}{\displaystyle \frac{\overline{\theta _i}}{(x\lambda _i\theta _i\overline{\theta }_i)^2}}.`$
In the planar limit, the terms involving the singlet operators coming from the variation of the action which we have denoted by

$`𝒮_\theta {\displaystyle \underset{i}{}}{\displaystyle \frac{S_{\mathrm{eff}}}{\theta _i}},𝒮_{\overline{\theta }}{\displaystyle \underset{i}{}}{\displaystyle \frac{S_{\mathrm{eff}}}{\overline{\theta _i}}},𝒮_{\theta \overline{\theta }}{\displaystyle \underset{i,j}{}}{\displaystyle \frac{^2S_{\mathrm{eff}}}{\theta _j\overline{\theta _i}}},𝒮_{}^{\mathrm{sub}}𝒮_{}𝒮_{}`$ (14)
factorize and become ignorable. We obtain

$`2\psi (x)\omega (x)+W^{}(x)\psi (x)+lim\rho _{\stackrel{ˇ}{N}}(x)=0,`$
$`2\overline{\psi }(x)\omega (x)+W^{}(x)\overline{\psi }(x)+lim\overline{\rho }_{\stackrel{ˇ}{N}}(x)=0,`$ (15)
$`2t(x)\omega (x)+2\overline{\psi }(x)\psi (x)+W^{}(x)t(x)+limc_{\stackrel{ˇ}{N}}(x)=0.`$
As is announced, the system of three planar loop obtained in eq.(S0.Ex24) takes the same form as that of the equations of Konishi anomaly in .

So far, we have dealt with the case in which the underlying gauge group is $`U(N)`$; this has led us to the supereigenvalue generalization of the hermitian one-matrix model. It is easy to generalize to the cases in which the underlying gauge group is $`SO(2N)/Sp(2N),SO(2N+1)`$. (See the recent discussion ). The measure factor is changed into

$`{\displaystyle \underset{i<j}{\overset{\stackrel{ˇ}{N}}{}}}(\lambda _i\lambda _j+\theta _i\overline{\theta }_i\theta _j\overline{\theta }_j)^2(\lambda _i+\lambda _j+\theta _i\overline{\theta }_i+\theta _j\overline{\theta }_j)^2,`$ $`\mathrm{for}SO(2\stackrel{ˇ}{N})`$
$`{\displaystyle \underset{i<j}{\overset{\stackrel{ˇ}{N}}{}}}(\lambda _i\lambda _j+\theta _i\overline{\theta }_i\theta _j\overline{\theta }_j)^2(\lambda _i+\lambda _j+\theta _i\overline{\theta }_i+\theta _j\overline{\theta }_j)^2{\displaystyle \underset{i}{\overset{\stackrel{ˇ}{N}}{}}}(\lambda _i+\theta _i\overline{\theta }_i)^2,`$ $`\mathrm{for}Sp(2\stackrel{ˇ}{N}),SO(2\stackrel{ˇ}{N}+1).`$
These cases can be handled by changing the set which integer labels $`i,j`$ have belonged to into the new ones:

$`i,j\{1,2,\mathrm{},\stackrel{ˇ}{N}\}`$ $``$ $`i,j\{\stackrel{ˇ}{N},\mathrm{},2,1,1,2,\mathrm{},\stackrel{ˇ}{N}\},\mathrm{for}SO(2\stackrel{ˇ}{N})`$
$`i,j\{1,2,\mathrm{},\stackrel{ˇ}{N}\}`$ $``$ $`i,j\{\stackrel{ˇ}{N},\mathrm{},2,1,0,1,2,\mathrm{},\stackrel{ˇ}{N}\},\mathrm{for}Sp(2\stackrel{ˇ}{N}),SO(2\stackrel{ˇ}{N})`$
and by demanding

$`\lambda _i=\lambda _i,\lambda _0=0.`$
The summations and the products in the previous equations are now taken with respect the integers belonging to these new sets. The potential effectively becomes a polynomial consisting only of even powers. The remainder of the discussions are the same.

Finally, let us discuss the supereigenvalue model which implements $`𝒩=2`$ super Virasoro constraints. The action is given by the $`𝒩=2`$ superfield (2). In the partition function (1), the measure factor is replaced by

$$\underset{i<j}{}(\lambda _i\lambda _j(\theta _i\overline{\theta }_j+\overline{\theta }_i\theta _j))^2,$$
(16)
where each factor is the super-distance invariant under the supertranslation generated by the supercovariant derivatives. The (super)-differential operators to be inserted in the partition function are

$`\widehat{j}(x)`$ $`=`$ $`{\displaystyle \underset{j}{}}\left({\displaystyle \frac{}{\theta _j}}\theta _j{\displaystyle \frac{}{\overline{\theta }_j}}\overline{\theta }_j\right){\displaystyle \frac{1}{x\lambda _j}},`$
$`\widehat{g}(x)`$ $`=`$ $`{\displaystyle \underset{j}{}}\left({\displaystyle \frac{}{\theta _j}}\overline{\theta }_j{\displaystyle \frac{}{\lambda _j}}\right){\displaystyle \frac{1}{x\lambda _j\theta _j\overline{\theta }_j}},`$ (17)
$`\widehat{\overline{g}}(x)`$ $`=`$ $`{\displaystyle \underset{j}{}}\left({\displaystyle \frac{}{\overline{\theta }_j}}\theta _j{\displaystyle \frac{}{\lambda _j}}\right){\displaystyle \frac{1}{x\lambda _j\overline{\theta }_j\theta _j}},`$
and $`\widehat{\mathrm{}}(x)`$ (the case $`\gamma =\frac{1}{2}`$) defined in eq. (5). They can be organized into a current superfield. The same is true for the resolvents of this model:

$$\mathrm{\Omega }(x,\eta ,\overline{\eta })g_s\underset{j}{}\frac{(\eta \theta _j)(\overline{\eta }\overline{\theta }_j)}{x\lambda _j}.$$
(18)
Introducing the supercovariant derivatives $`D\frac{}{\eta }+\overline{\eta }\frac{}{x},\overline{D}\frac{}{\overline{\eta }}+\eta \frac{}{x}`$, we have been able to derive a manifestly $`𝒩=2`$ supersymmetric loop equation. Space only permits us to present its final form:

$$0=D\mathrm{\Omega }\overline{D}\mathrm{\Omega }g_s\mathrm{\Omega }^{}+\left(D𝐕\overline{D}\mathrm{\Omega }+\overline{D}𝐕D\mathrm{\Omega }\right)+𝐅.$$
(19)
Here $`𝐅(x,\eta ,\overline{\eta })`$ is a polynomial in $`x`$ given explicitly in terms of the components of $`𝐕`$. Factorization of the two-point functions in the large $`\stackrel{ˇ}{N}`$ limit ($`g_s0`$) gives us an $`𝒩=2`$ planar loop equation.

We plan to give full details of the contents of this letter together with other issues in a future publication. We thank Alyosha Morozov for helpful discussion on this subject.