text:
# Scaling Limits for the 2D Metal-Insulator Transition at 𝐵=0 in Si-MOSFETs

## I Introduction

Several years ago the observation of a conductivity phase transition in the two-dimensional electron system (2DES) in Si metal-oxide-semiconductor field-effect-transistors (MOSFETs) was reported. This observation was unexpected and controversial because it was contrary to the widely accepted scaling theory for a noninteracting 2DES. Further experimental studies on MOSFETs, SiGe heterostructures, and GaAs heterostructures indicate that a conductivity transition is possible as a function of carrier density in a 2DES at $`B=0`$. The observed conductivity transitions at $`B=0`$ bear a marked similarity to transitions between states with different Hall resistances in the quantum Hall effect (QHE). This unexpected behavior and similarity to QHE transitions has also sparked considerable theoretical interest. Furthermore, these parallels make it advisable to consider developments in the QHE when studying the $`B=0`$ conductivity transition in 2DES. In particular, a recent publication by Shahar et al.reports that the behavior of the resistivity as a function of temperature and electron density, $`\rho (T,n_s)`$, is inconsistent with scaling at temperatures above 0.5 K. Thus, they question the use of scaling analysis for the study of the QHE in general. We feel that a serious consideration of the appropriateness of the scaling analysis for the $`B=0`$ conductivity transition in 2DES is warranted. We have data available that can address this question using a set of data manipulation routines in Mathematicatm especially developed for this purpose. Below we briefly describe the data used in this study and the data manipulation routines that allow us to consistently propagate errors and eventually place confidence limits on the applicability of this scaling analysis to the case of the $`B=0`$ conductivity transition in 2DES. We feel that these methods have applicability beyond the particular system to the QHE and other quantum phase transitions.

Our aim has been to develop an objective method that gives a quantitative answer to the question of how well the data actually scales. We achieve this by analyzing our ability to fit the data to a general scaling model that assumes as little as possible about the functional form of the scaling function. We use the variance-weighted mean-square deviation from this fit, $`\chi ^2`$, as a quantitative measure of the ’goodness’ of the fit. Such a properly weighted and normalized $`\chi ^2`$ close to one indicates good agreement of the data with the model. An objective estimate of $`\chi ^2`$ requires a quantitative knowledge of experimental uncertainties. Therefore, great care has been taken to estimate the magnitude of these sources of uncertainty and to propagate them consistently to the final fit.

Using this method to investigate $`E`$-field scaling of the lowest temperature data available, we find rather good agreement with scaling in the best sample with a $`\chi ^2`$ varying from 0.4 to 1.9 depending on our selection of the data. We find that $`\chi ^2`$ gets smaller as more data at low power is discarded. Excluding more data at lower powers also tends to increase the consistency among the different samples.

## II Description of the data including sources of error

The data used in this study was amassed at the University of Oklahoma while studying the temperature $`T`$, electron density $`n_s`$, and excitation dependence of the resistivity $`\rho `$ of a number of high-mobility $`(\mu _{\mathrm{max}}25\mathrm{`}000\mathrm{cm}^2/\mathrm{Vs})`$ metal-oxide-semiconductor field-effect transistors (MOSFETs). These devices are described in detail elsewhere.

This data includes large sequences of currents $`I`$ with their corresponding probe voltages $`V`$ recorded as a function of gate voltage $`V_gn_s`$ and $`T`$ using a four-terminal technique. The resulting $`I`$$`V`$curves at low temperatures are highly nonlinear with distinctly different behaviors above and below a critical density $`n_c`$ or resistivity $`\rho _c`$. In this article we consider the nonlinear $`I`$$`V`$curves taken at the lowest temperature ($`0.2`$ K) and present a scaling analysis of the resistivity as a function of the applied electric field, $`E`$. This $`E`$ is taken to be the measured $`V`$ divided by the distance between the respective probes.

We find three sources of randomness that we assume are independent:

* The uncertainty in the measured probe voltage characterized by a root mean square (RMS) uncertainty $`\sigma _V`$.
* The uncertainty in the measured excitation current, $`\sigma _I`$, which we find is negligible.
* The uncertainty in $`n_s`$ due to the applied $`V_g`$ similarly characterized by $`\sigma _n`$.

For these samples at this $`T`$ we also have a source of systematic error that we associate with the well-known difficulty in producing ohmic contacts to a 2DES. The non-ohmic contacts caused the measured $`I`$$`V`$curves to be asymmetric with respect to $`I=0`$ and to show a sharp anomaly centered at $`I=0`$. By carefully accounting for all the voltage drops in the measurement circuit, we were able to determine the contribution due to the contacts. We found that a percentage of this contact potential (proportional to the resistivity of the sample) added to the measured $`V`$ essentially eliminated the observable effects of the contacts, see Fig. 1.

Furthermore, we found that these contact potentials did not effect our analysis because they occurred at very low excitation currents where the physics of the system is still dominated by the finite $`T`$. The results presented here are independent of these contact anomalies and our treatment of them.

In order to measure the effective $`\sigma _V`$ for our data, we produced a modified data set that was equally spaced in excitation $`I`$ by using a spline interpolation scheme. The Fourier transform of the data was then examined. The transformed data is consistent with a low frequency signal superimposed on uniform white noise of constant amplitude. This white noise amplitude varies from about 2 $`\mu `$V in the most metallic curves up to 90$`\mu `$V in the most insulating ones. We used this amplitude as a measure of $`\sigma _V`$ for each particular $`V_g`$ in the following analysis.

The uncertainty in the electron density, $`\sigma _n`$, is more difficult to estimate. To calibrate $`n_s`$ to $`V_g`$, we sweept the gate voltage at a constant magnetic field to obtain $`\rho _{xx}(V_g)`$. Minima in $`\rho _{xx}`$ that correspond to known integral filling factors, and therefore known $`n_s`$, were then determined. The linear fit to the resulting $`(n_s,V_g)`$ pairs that gives $`V_g(n_s)`$ with a statistical error of $`\sigma _n=110^9`$ cm-2. This corresponds to a resolution of 0.8 millivolt in $`V_g`$. This determination is not obviously a measure of the density fluctuations for constant $`V_g`$. Furthermore, an effective gradient in the gate voltage will be imposed on the 2DES by the source-drain excitation. We made every attempt to make this potential drop symmetric with respect to the $`V`$ probes, but the $`n_s`$ gradient could not be fully eliminated. Thus, we estimate the density uncertainty, $`\sigma _n`$, by taking $`\sigma _{V_g}1\text{m}V`$. This estimate is compatible with our observations of the precision and reproducibility of producing $`n_s`$ with $`V_g`$ in our samples as evidenced by a repeatability of a given $`\rho _{xx}(V_g)`$ where $`\rho _{xx}(n_s)`$ is particularly steep.

Qualitatively, as has been reported previously,  two different nonlinear $`V(I)`$ behaviors can be seen. At low densities the resistivity has a maximum at zero electric field and decreases for stronger fields. This is similar to a normal insulator where a larger electric field can excite more carriers and hence increase the conductivity. At high electron densities the opposite behavior is seen; the resistivity has a minimum at $`E=0`$ and increases for stronger fields. Between these behaviors, a linear $`V(I)`$ is observed with a constant $`\rho `$, see Fig. 2.

This is similar to the nonlinear region of a superconductor. However, the nonlinearities seen here are anomalously strong, they occur at very small fields, and they are unexplained by a clear microscopic model.

Because the $`I`$$`V`$curves are nonlinear, there is a difference between the normal resistivity given by the ratio $`V/I`$, and the differential resistivity, $`\frac{dV}{dI}`$. We have analyzed both, and they have qualitatively similar behaviors. The differential resistivity has the advantage of being insensitive to experimental voltage and current offsets that must be considered when calculating ratios. The $`\frac{dV}{dI}`$ is also more sensitive both to structure in the $`I`$$`V`$and to noise. We generate differential resistivities by doing a n-point differentiating convolution of the $`I`$$`V`$data;

$$\rho _i=\underset{j=1}{\overset{n}{}}c_j(V_{i+j}V_{ij}),$$
(1)
where the data points, $`\{V_i\}`$, are made equidistant in current excitation by spline interpolation. Choosing different Super-Lancos differentiation filters, $`\{c_n\}`$, allows us to optimize the data smoothing to reduce noise while retaining the signal. We then propage $`\sigma _V`$ in the original data through the convolution to obtain the variance in the resistance, $`\sigma _\rho `$.

## III Scaling model

The question at issue is whether or not a quantum ($`T=0`$) phase transition between two different conductivity states takes place at $`B=0`$ in these systems. As mentioned in the introduction this is still highly debated in spite of large experimental and theoretical efforts. Our aim here is to develop a more rigorous and objective way to determine the consistency of experimental data with such a phase transition.

A quantum phase transition occurs at zero temperature whereas all conductivity measurements are taken at both finite bath temperature and nonzero excitation voltage to drive a current through the sample. Theoretical arguments suggest that an underlying quantum phase transition should give rise to scaling at finite temperatures and excitation $`E`$-fields,  i.e. that the resistivity of the sample should have a specific temperature and $`E`$-field dependence;

$$\rho (n_s,T,E)=\rho (x,y)=f(\frac{\delta }{E^{1/\beta }},\frac{\delta }{T^{1/\nu z}}),$$
(2)
where $`\delta =\frac{n_sn_c}{n_c}`$ is the relative deviation of the electron density from a critical density $`n_c`$, and $`\nu `$, $`z`$ and $`\beta `$ are critical exponents associated with the phase transition. Assuming the simplest $`E`$-field heating mechanism, $`\beta `$ would be related to the temperature exponents as $`\beta =\nu (z+1)`$. Here $`x=\delta /E^{1/\beta }`$ and $`y=\delta /T^{1/\nu z}`$ are called the scaled variables.

In cases where either the $`E`$-field or the temperature has been judged negligible compared to the other, data has been examined for one-parameter scaling, i.e. whether or not it collapses to a function of only one scaled variable. At the core of all arguments for a phase transition lies the fact that one in many cases can find a critical density $`n_c`$ and exponent, $`\nu z`$ or $`\beta `$, so that data falls onto a single curve. This is usually done by trying out different critical densities and exponents. The best fit is then determined by eye as the parameter choice that gives the best data collapse. Because the question of a phase transition in our system still is controversial, we developed an objective method to determine the compatibility of the data with scaling.

Our goal is now to give a quantitative measure of the consistency of the data with a single curve when we plot the resistivities taken at the lowest temperature as a function of the scaled $`E`$-field, $`x=\frac{(n_sn_c)}{n_c}E^{1/\beta }`$.

We want to assume as little as possible about the form of the scaling function, $`f(x)`$, but previous investigations suggest that $`Log[f(x)]`$ is a reasonably smooth curve without too much curvature. We therefore investigate a model where $`Log[f(x)]`$ is a polynomial with only a few terms, $`f(x)=\mathrm{exp}[a_0+a_1x+a_2x^2+a_3x^3+\mathrm{}]`$. Because theory does not tell us the form of $`f(x)`$, any functional form for $`f`$ that is general enough is acceptable, provided it contains a small number of free parameters compared to the number of data points being fitted. If we truncate this polynomial at say third order, we obtain a scaling model for $`\rho (n_s,E)`$ that contains the six parameters $`n_c,\beta ,a_0,a_1,a_2`$ and $`a_3`$;

$`\rho _{model}(n_s,E)`$ $`=`$ $`e^{a_0+a_1x+a_2x^2+a_3x^3}\text{where}`$ (3)
$`x`$ $`=`$ $`{\displaystyle \frac{(n_sn_c)/n_c}{E^{1/\beta }}}.`$ (4)
We can now address the question of the consistency of the data with scaling by investigating our ability to fit the measured data to this scaling model. We do this with a nonlinear regression routine in MathematicaTM that fits all parameters at once and hence obtain values for $`n_c`$ and $`\beta `$ with statistical error estimates. At the same time we get an objective measurement of the ’goodness of fit’ via the $`\chi ^2`$ of the fit.

In order to have an unbiased fit and get $`\chi ^2`$ in the correct units we need to keep track of the estimated uncertainty for each experimental point. In this case we need to consider not only the error in $`\rho `$ itself but also the uncertainty in $`n_s`$ and in $`E`$. We find that the uncertainty in density, $`\sigma _n`$ normally dominates the uncertainty of points in the $`f(x)`$-plot. To take this into account we propagate all errors to a single effective error in the resistivity, $`\sigma _{\rho ,\text{eff}}`$:

$`\sigma _{\rho ,\text{eff}}^2`$ $`=`$ $`\left[\sigma _\rho \right]^2+\left[\sigma _n{\displaystyle \frac{\rho _{model}}{n_s}}\right]^2+\left[\sigma _V{\displaystyle \frac{\rho _{model}}{V}}\right]^2.`$ (5)
We estimate the error in the resistivity $`\sigma _\rho `$ by propagating the errors from the probe voltage $`\sigma _V`$ either via the convolution (1) or the ratio $`V/I`$. When calculating the total effective error $`\sigma _{\rho ,\text{eff}}`$, we still assume $`\sigma _\rho `$ and $`\sigma _V`$ to be independent. We then weigh the points in the fit with a weight $`1/\sigma _{\rho ,\text{eff}}^2`$. Because the weights depend on the fitted polynomial, we have to start with an unweighted fit, then propagate the errors and iterate the fit until a self-consistent solution is found.

## IV Results

Using either the derivative $`dV/dI`$ or the ratio $`V/I`$ as the resistivity $`\rho `$, we have used the scaling model described in the previous section to investigate $`E`$-field scaling in three different samples. Here we used a fifth order polynomial (3). We also achieved reasonable results with lower order polynomials. Because data is taken at a finite temperature, ($`0.2`$K) we know that the resistivity is dominated by the temperature for small enough excitation $`E`$-fields; hence, we introduce a low power cut-off to look for scaling in the $`E`$-field only.

Using the ratio $`V/I`$ for the resistivity $`\rho `$ we obtain a $`\chi ^2`$ reasonably close to one. Figure 3 shows the result from sample one,

using a low power cut-off $`P_c=50`$ pW. The critical exponent we find in this way is comparable to previous investigations, even though the exponent and $`\chi ^2`$ do depend somewhat on the particular cut-offs used.

We have most data from sample one. It also shows the best behavior; $`\chi ^2`$ is about one and the physically interesting fit parameters, in particular the exponent $`\beta `$, are relatively stable as the low-power cut-off is increased, see Table I.

The results in the other samples show stronger dependence on the cut-off. Only for quite large low power cut-offs ($`500`$pW) are the exponents from different samples compatible with each other, see Table II

for results from sample two. We note that at an excitation power of 1 nW temperature increases in the sample were observed via the temperature monitoring RuOx resistor. Thus, data at these high excitations may correspond to the hydrodynamic regime discussed in this context earlier.

Sample one shows rather good consistency with the scaling model. Here all the fitted parameters are reasonably independent of the low power cut-off, $`P_c`$. The number of points in a fit varies from 864 for the smallest cut-off to 467 for the largest. Actual variations in the fitted parameters are consistent with the uncertainty estimates based on the statistical deviations from the model as the low power cut-off is changed. For example, the critical density $`n_c`$ is virtually unchanged for all cut-offs. This is also reflected in very small statistical uncertainties for $`n_c`$ in each fit. The exponent $`\beta `$ (=$`\nu (z+1)`$) is consistent with $`\beta =3.8\pm 0.05`$, and the critical resistivity $`e^{a_0}`$ is approximately $`2.9h/e^2`$. The mean square deviation per point, $`\chi ^2`$, however, does vary substantially with the cut-off, decreasing as the cut-off is increased. This could be understood as an effect of the finite temperature the data is taken at; at low powers the temperature is not negligible and the model (3) is too crude to give a good fit. It is also possible that there are two regimes of nonlinearity as mentioned above. The fact that $`\chi ^2`$ is less than one for the best fits may indicate that we have overestimated the statistical uncertainties in the data.

Sample two shows a larger dependence on the low power cut-off, and only at large cut-offs ($`400`$pW) is the exponent $`\beta `$ compatible with the results in sample one. Clearly the exponent varies much more than the statistical error estimates in each of the fits with lower cut-offs. Again $`\chi ^2`$ varies in the same way as in sample one. At the largest cut-offs, however, both $`\beta `$ and $`\chi ^2`$ are close to the values found in sample one. The critical resistivity is, however, clearly not the same as in sample one, $`e^{a_0}7.5h/e^2`$.

Sample three shows a behavior similar to sample one (exponent $`\beta 3.7`$ and $`e^{a_0}3.9`$) but with much poorer statistics.

In order to gain a further understanding of the data, and our fitting procedure, we have also studied the effect of introducing a cut-off in the scaled variable $`x`$. By keeping data where $`|x|x_{\text{cut}}`$ we effectively cut away points that are taken far from the critical density. Note that the criteria for a specific density to be “near” or “far” from the critical density is determined also by the excitation field; a small field corresponds to a more narrow range of densities. Since $`x`$ depends also on the critical density and exponent given by the fit, this cut has to be made self-consistently. Ideally one would expect fitted parameters to be independent of the actual value of $`x_{\text{cut}}`$, provided $`x_{\text{cut}}`$ is not too large. Figure 4

shows how the fitted exponent $`\beta `$ varies with $`x_{\text{cut}}`$ in sample one. Each point in the figure corresponds to a complete fit using both the $`|x|x_{\text{cut}}`$ cut-off and a low power cut-off ($`P_c`$). The statistical uncertainty in each fit, shown as errorbars on the points, are clearly smaller than the actual variation in $`\beta `$. This behavior may be an indication that the idea of multiple regimes is appropriate, but further investigations are clearly necessary to settle the question. The large uncertainties in the points with very small $`x_{\text{cut}}`$ reflect the fact that there are no longer enough data left for a good fit.

We have also investigated the derivative $`dV/dI`$ for scaling. In this case we find similar but more varied results with generally higher $`\chi ^2`$. The $`\chi ^2`$’s determined for these fits also depend on the specific Super-Lancos filter used (1). The derivative scaling also requires the introduction of cut-offs for high power data. This was done by discarding data outside three times the FWHM (Full Width Half Maximum) of the central structure. Results depend strongly on this FWHM-cut. The results, presented above, based on ratios show very little dependence on such high power cut-offs and we have used all data consistent with other cut-offs.

## V Conclusions

We have developed and used a general scaling model to fit all parameters in $`E`$-field scaling at once. We find some results that are consistent with scaling, $`\chi ^2`$ of about one, but also surprisingly large variations in the fitted parameters when different physically appropriate cut-offs are introduced.

It would be desirable to generalize and use this method for a full two-parameter scaling test, treating the $`E`$-field and temperature at equal footing. Any one-parameter scaling suffers from the fact that the other (unscaled) parameter is not zero. As we have discussed above, any $`E`$-field scaling will always suffer from the finite temperature of the surrounding bath. In the same way, temperature scaling is sensitive to the nonzero excitation field used to measure resistivities, and heating of the 2DES will always be a problem for low enough temperatures. Both these problems are reduced in a two-parameter scaling test where the scaling model (2) is supposedly valid for any $`E`$ and $`T`$ such that $`x`$ and $`y`$ are close enough to the transition.

## VI Acknowledgements

The data used in this study was amassed at the University of Oklahoma between 1993 and 1996 under the auspices of the Laboratory for Electronic Properties of Materials, which was supported by NSF grants EPSCoR-92-OU-NSF and EPSCoR-95-OUMATERIALS, and The Department of Physics and Astronomy, which was supported by NSF grants DMR-89-22222 and DMR-96-264699. Whitney Mason, Sheena Q. Murphy, and Sergey V. Kravchenko were active collaborators with JEF in collecting this data. We wish to thank Kieran Mullen, Bruce Mason, Anders Karlhede and Christian Walck for many useful discussions. We also wish to thank Shankar Das Sarma and the other participants of the 1998 Program on Disorder and Interactions in Quantum Hall and MesoScopic Systems at the University of California, Santa Barbara for suggesting this project and providing much of its initial impetus.


format: text

text:
# Kinetic roughening with anisotropic growth rules

## I Introduction

In the last few years the study of physical phenomena characterized by a degree of self-organization , has attracted a lot of interest. These models are usually cellular automata models defined on a discretized lattice, with a growth rule that can be either stochastic, when the inhomogeneities in the system change with a time scale smaller than the characteristic time scale of the dynamical evolution (noise), or deterministic with a quenched disorder which accounts for the effect of inhomogeneities inside a solid medium. Both kind of dynamical rules are characterized by an evolution towards an attractive fixed point in which scale free fluctuations in time and space are present .

The problem of kinetic roughening belongs to this class of models. It received recently an increasing interest in relation with non-equilibrium growth models  and in view of its practical applications: Chemical Vapor Deposition (CVD)  and electro-chemical deposition  are just two examples.

In this perspective, one is interested in identifying the dynamic universality classes of kinetic roughening processes and several models has been defined starting from the models falling in the universality class of the Kardar-Parisi-Zhang (KPZ) equation . This equation describes the properties of an interface $`h(x,t)`$ driven by a stochastic noise and gives a roughness exponent $`\chi =0.5`$. Other models are more suitable to describe the propagation of interfaces in random media, i.e. with a quenched disorder. These models are driven by an extremal dynamics. In this class fall the so-called Sneppen model  (in  referred as model B) and the pinning model by directed percolation , which predict a roughness exponent equal to $`\chi =0.63`$. These models produce self-affine surfaces. Recently, a model has been introduced to describe some etching experiments, which leads to the formation of self-similar (fractal) structures, and which has been shown to fall in the percolation universality class . Many experiments on surface roughening , however, as well as experiments on chemical etching  produce self-affine surfaces instead of self-similar ones. In this paper, we will focus on kinetic roughening phenomena leading to self-affine ($`\chi <1`$) surfaces.

We recall that the roughness exponent is defined by the ensemble averaged width of the interface as $`W(l,t)=<(h(x,t)<h(x,t)>)^2>^{1/2}l^\chi f(t^{1/z}/l)`$ where $`z`$ is the so-called dynamical exponent, the angular brackets denote the average over all segments of the interface of length $`l`$ and over all different realizations. $`f`$ is a scaling function such that $`f(y)y^\chi `$ for $`y<<1`$ and $`f(y)=cost.`$ for $`y>>1`$. The exponent $`\beta =\chi /z`$ describes the transient roughening, during which the surface evolves from the initial condition toward the final self-affine structure.

In spite of the strong universality exhibited by the KPZ and Sneppen models, in many experimental studies one measures values of $`\chi `$ which are above the ones predicted by both the KPZ and Sneppen universality classes. To give some examples, we remember the experimental studies of paper burning for which one gets $`\chi =0.70\pm 0.03`$ , or the propagation of a forced fluid front in a porous medium, which exhibits a roughness exponent $`\chi =0.73\pm 0.03`$  and $`\chi =0.88\pm 0.08`$ .

In this paper we propose a generalized model for kinetic roughening characterized by anisotropic growth rules and, as a consequence, separated time scales for the dynamics. The existence of this time-scale separation induces a cross-over effect in the roughness properties, which could erroneously appear as a genuine non-universal critical behavior, and could give an explanation for the above cited scattered experimental results.

Some results presented here have been already briefly reported on in a letter . In this long paper we give a detailed, complete description of our previous work. Moreover, we present a set of new numerical results, which allow us to reach different and better founded conclusions with respect to .

The idea underlying the model is that some experimental parameters can introduce a characteristic scale in the system, separating different scaling behaviors. In particular we consider a model which includes explicitly an anisotropy factor, say a growth rule dependent on the local environment of the growing site. The model thus presents a complex interplay between a global equilibrium and the conditions of a local dynamics. This choice is motivated by the observation of roughening phenomena occurring in etching processes which represent an important tool either in academic research or in device technology. Their importance is related to the preparation of single-crystal samples of desired dimensions, shapes and orientations. Etching is usually applied to obtain desired mesas and grooves in semiconductors wafers and multilayers .

In the same field, although in a different context, etching processes are used to produce textured optical sheets, which allow to exploit the light trapping by total internal reflection to increase the effective absorption in the indirect-gap semiconductors crystalline silicon. Light trapping, originally suggested to increase the response speed of silicon photo-diodes while maintaining high quantum efficiency in the near-infrared, was later indicated as an important benefit for solar cells .

The general suffix etching indicates the ensemble of operations which involve the removal of materials by expending energy either by mechanical, thermal or chemical means. In ref.  the authors focused their attention on chemical etching processes as a reference point to formulate the model. One of the most important properties of these processes is represented by the intrinsic anisotropy  of the etch rates. For instance in samples of crystalline silicon etched in solutions of aqueous potassium-hydroxide (K-OH) with isopropil alcohol (IPA), depending on the concentration of the etchant and the temperature, the $`(111)`$ direction etches slower than the others by a factor which can be of order $`100`$ or more . The degree of anisotropy affects the properties of the surface, which turns out to be rough with an apparently non-universal roughness exponent.

Although the definition of the model is very general we will briefly consider the chemistry of the etching process in order to exhibit a physical framework that allows to understand the meaning of the definitions and their interpretation.

The disorder in the etching process is related to the impurities in the lattice. Such impurities, e.g. vacant atoms, reduce the binding energy of atoms nearby the vacancy. By assigning to each site (atom) of our lattice a random number $`x_i`$ we assume that a distribution of vacancies, or other kinds of impurities, is present in the system, and this induces fluctuations in the binding energy of atoms due to this disorder. If we assume to be in a condition of slow dynamics, that is to say the driving field (which in our case in represented by the concentration of etchant) tends to zero , we can look at the etching as an extremal process, where the etchant dissolves the atom with the smallest binding energy. This is correct for low etchant concentrations and corresponds actually to the situation experimentally more interesting, in which rough surfaces are produced.

In order to reproduce the experimental conditions (type of etchant, concentration, temperature) a microscopic model for the physical process should contain some tunable parameters (at least one). In our model, the anisotropy is introduced by a phenomenological tunable parameter, $`p`$, which distinguishes sites with a different local environment.

The introduction of the parameter $`p`$ defines a characteristic scale in the problem. As a result the critical properties of the model are characterized by a continuous crossover between two universality classes corresponding to the roughness exponents $`\chi =1`$ and $`\chi =0.63`$ (Sneppen models $`A`$ and $`B`$ respectively ). In particular one can define a parameter $`r=\frac{p}{1p}`$ which measures the time-scale separation between the dynamics of the different classes of sites. For lengths $`l>l^{}\frac{1/2p}{p}`$ one observes a behavior characteristic of the Sneppen model $`B`$ universality class ($`\chi =0.63`$), while for $`l<l^{}\frac{1/2p}{p}`$ one observes a behavior characteristic of the Sneppen model $`A`$ universality class. The existence of this crossover is difficult to detect directly on the plot for the scaling of the surface width $`W^2(l)`$ (especially for finite sets of data) and it becomes evident looking at the power spectra. This explain why large-scale experiments could give the impression of non-universality in the critical properties of rough surfaces.

Let us look at the meaning of $`p`$ in the case of etching. If we represent the crystalline lattice of the silicon on a two-dimensional plane we can imagine a square lattice where (see Fig.(1)) the atoms can be found in each of the four positions marked in figure by the letters $`ad`$. The four positions correspond to different oxidation states: from the situation $`(c)`$ (oxidation number $`0`$) which occurs only in the bulk, to the situation $`(d)`$ (oxidation number $`3`$). Note that all the surface atoms are passivated by hydrogen atoms. The atoms in the positions $`(a)`$ and $`(b)`$, corresponding respectively to the oxidation numbers $`2`$ and $`1`$ (two and one heteropolar bonds, i.e. Si-H bonds), play an important role in explaining, at least from a heuristic point of view, the origin of the anisotropy in the etched rates . The parameter $`p`$ quantifies the ratio of the etch rates between the sites in the positions $`(a)`$ and $`(b)`$. The basic idea is that in the $`(111)`$plane of silicon, there is only one heteropolar bond per silicon atom. Therefore there are three bonds to break for dissolution, while other planes (except the $`(110)`$) have more than one heteropolar bonds and accordingly a smaller number of bonds must be broken.

The paper is organized as follows. In section II we describe in detail the model and the setup of numerical simulations. In section III we present and discuss the numerical results in relation with the Sneppen model A and B universality classes. In section IV a discussion of the results and some conclusions are drawn together with a planning of future researches.

## II The model

We give now a detailed definition of the model. The model is defined on a square $`2D`$ lattice tilted at $`45^{}`$ (see Fig. 1). We consider a $`1+1`$ dimensional interface $`h(x)=h(x,t)`$ with $`x=1,2,\mathrm{},L`$, where $`L`$ is the linear extension of the interface in the $`x`$ direction. The initial condition for the dynamical evolution of this interface is given by: $`h(2x,0)=1x[1,L/2]`$ and $`h(2x1,0)=0x[1,L/2],`$ in order to have both classes of variables (lattice planes) participating to the dynamics from the beginning, but different initial conditions do not change the properties of the model. The interface, which satisfies locally the conditions

$$|h(x,t)+1h(x1,t)|1,$$

$$|h(x,t)+1h(x+1,t)|1,$$
(1)
contains two classes of random variables that correspond to two separate classes of sites. The sites $`(M)`$ for which it holds $`^2h>0`$ (called minimum sites) which are, microscopically, the atoms with two heteropolar bonds, and the sites on a slope (slope $`(S)`$ sites) for which one has $`^2h=0`$. These last sites correspond microscopically to atoms with one heteropolar bond. To each class of sites is assigned a class of Gaussian distributed uncorrelated random variables which mimic the disorder, and represents physically, for the case of etching, the binding energy of atoms:

$$\eta (x,h)\{\begin{array}{cc}\left[0:0.5\right]\hfill & \text{if x is such that}^2h=0\\ \left(0.5:1\right]\hfill & \text{if x is such that}^2h>0.\end{array}$$
(2)
The sites with $`^2h<0`$, for which all the chemical bonds are homopolar, i.e. Si-Si, do not take part to the dynamics and they have assigned a zero value of the random variable. Periodic boundary conditions are assumed along the $`x`$ direction.

The system evolves by updating the site $`i^{}`$ with the largest r.v. in one of the two classes of sites chosen, at its turn, with a probability $`p`$. One thus updates with probability $`p`$ a site $`(S)`$ and with probability $`1p`$ a site $`(M)`$ according to the rules (see Fig. (2)):

* $`h(i^{},t+1)=h(i^{},t)+2`$, $`\eta (i^{},h(i^{},t+1),t+1)=0`$;
* Updating all the sites necessary to make satisfied the conditions 1 (this phase is assumed to be instantaneous with respect to extremal dynamics);
* Updating of the random variables for the sites which changed their class of belonging. In particular $`\eta (x,h,t+1)=1/2RAN`$ if $`\eta (x,h,t)=0`$ and $`\eta (x,h,t+1)=\eta (x,h,t)+1/2`$ if $`\eta (x,h,t)0`$, where $`RAN`$ is a random value between $`0`$ and $`1`$;
* Updating of the random variables of the sites which have changed their height but which did not change their class of belonging, (sites S only): $`\eta (x,h,t+1)=1/2RAN`$.

The parameter $`p`$ can vary in the range $`[0:1/2]`$. If we define $`t_S`$ as the characteristic time scale for $`S`$ variables and $`t_M`$ the characteristic time scale for $`M`$ variables one has:

$$r=\frac{t_M}{t_S}=\frac{p}{1p}$$
(3)
The growth of the interface, in the etching process, represent the invasion of the etchants into the silicon wafer. ¿From this point of view the updating of the sites $`(S)`$ mimics the etching of the $`(111)`$ planes and the updating of the $`(M)`$ sites the etching in the $`(100)`$ direction. For $`p=1/2`$ all the sites which take part to the dynamics are equivalent and there is no anisotropy ($`r=1`$), whereas the case $`p=0`$ corresponds to the maximal anisotropy in which $`v_{(111)}/v_{(100)}=r=0`$, where $`v_{(111)}`$ and $`v_{(100)}`$ are the etch rates in the corresponding directions.

Our model can be viewed as a variation of the Sneppen model for quenched surface growth, where two important elements are added: 1) The anisotropy in the distribution of the quenched random field, depending on the local characteristics of the growing surface: 2) a time scale separation for the dynamical evolution of the two classes of variables ($`S`$ sites and $`M`$ sites), which is tuned by the parameter $`p`$.

## III Numerical results

We have studied this cellular automata by numerical simulations in order to analyze its dynamical roughening properties. The sizes we have chosen for the numerical simulations range from $`L=2048`$ to $`L=8192`$. For each value of $`p`$(we have considered $`p=0.0,0.02,0.2,0.5`$) , $`10^2`$ simulations lasting $`10^7`$ time steps have been performed and we have computed the growth exponent $`\beta `$, which rules the time evolution of the width $`W(t)`$ of the surface ($`W(t)t^\beta `$) before the stationary state is reached, and the roughness exponent $`\chi `$, which gives the scaling of the width of the surface ($`W(l)l^\chi `$), in the stationary state. The stationary state is called self-organized in that it is reached spontaneously by the system independently of the initial conditions. This self-organization is confirmed by an analysis of the temporal evolution of the distribution of quenched variables (the histogram $`\mathrm{\Phi }_t(\eta )`$). To characterize temporal correlations in the dynamics and check that the asymptotic state is critical, we studied the distribution of the avalanches in the asymptotic state. As an independent check about the universality of the roughness properties of the model, we have studied numerically the power spectrum $`S(k)`$ of the height profile.

To ensure that the system is in the stationary state, we studied the behavior of the $`n`$-th moments (for $`n=3,4,5`$) of $`h(x,t)`$ normalized to second momentum:

$$m_n(t)=\frac{_i(h(i,t)\overline{h}(t))^n}{\left(_i(h(i,t)\overline{h}(t))^2\right)^{\frac{n}{2}}},$$
(4)
where $`\overline{h}(t)`$ is the mean surface height at time $`t`$. We get that, after a transient, all the odd moments vanish and the even ones tend to constant values (see Fig.s 3-5). In particular the condition for the skewness $`m_3=0`$ (Fig. 3), which characterizes the stationary critical state, is realized after about $`10^3`$ time steps per site, independent of the value of $`p`$. These results imply that the higher moments of the variable $`h(x,t)h`$ scale in a trivial way (they are powers of the second moment), and after the transient the probability distribution of the variable $`h(x,t)`$, which can be viewed as a random variable, is Gaussian. The amplitude of the normalized even moments $`m_{2n}`$, in the asymptotic stationary state, characterizes the roughness properties of the interface.

In Table I we report the measured values for the dynamical exponents $`\beta `$ which turns out to be independent of $`p`$. Fig.s (6-9) show the scaling behavior of $`W^2(l)`$ for different values of $`p`$. For $`p=0.0`$ (i.e. maximal anisotropy) the measured values of $`\chi `$ are affected by a finite-size effect and they tend, in the limit $`L\mathrm{}`$, to the value $`\chi =1.0`$ found for Sneppen model A  (Fig. 6). In this case the surface is composed by very big pyramids (Fig. 10). On the other hand for $`p=0.5`$ one recovers the universality class of Sneppen model $`B`$ with $`\chi =0.63`$ (Fig. 9). For all the other values of $`p`$ between $`0`$ and $`0.5`$ ($`p=0.02`$ in Fig. 7 and $`p=0.2`$ in Fig. 8), trying to perform fits away from the saturation regions, one would be tempted to invoke the existence of a non-universal behavior ruled by the parameter $`p`$. A careful observation puts in evidence that the curves for $`W^2(l)`$ seem to exhibit a crossover between the Sneppen models $`A`$ and $`B`$ universality classes. On the basis of $`p`$ one can define a characteristic length $`l^{}\frac{1/2p}{p}`$ above which one could see the $`\chi =0.63`$ behavior and below which the $`\chi =1`$ behavior. We shall come back on these considerations later on when we shall discuss the power spectra.

We have also studied the time evolution of the distribution of random variables on the invading interface (the histogram $`\mathrm{\Phi }_t(\eta )`$, where $`\eta `$ is a generic value for $`\eta (x,h)`$), which is of great importance for models with extremal dynamics. The results of the simulation are shown in Fig.s (11-14). One can see that $`\mathrm{\Phi }_t(\eta )`$ self-organizes, for $`p0`$, after about $`10^5`$ time steps, into a distribution that is the superposition of two theta functions, one for each class of variables, each one characterized by a critical threshold $`p_c(p)`$ depending on the parameter $`p`$. The meaning of these thresholds is that only $`S`$ variables larger that $`\eta _c^S`$ and $`M`$ variables larger than $`\eta _c^M`$ can grow . For $`p=0`$, instead, the histogram has no self-organized critical state (Fig. 11). Looking carefully at Fig. (11) we can see that, while in the initial transient there are a few $`S`$ sites, in the asymptotic state most sites are $`S`$ sites. In fact nearly all variables larger than $`0.5`$ (the $`M`$ variables) are disappeared. This observation agrees with the actual structure of the surface, which is composed by very big pyramids (Fig. 10), with a roughness exponent $`\chi 1`$. This picture is confirmed by the acceptance profile $`a(\eta )`$, which is shown in Fig. 15. As in the Bak and Sneppen model, the acceptation profile (that is to say the distribution of the values of all updated quenched variables up to the actual time) exhibits correlation properties (it is not flat), reflecting temporal correlations in the dynamics. But, while the acceptation profile for $`S`$ variables is quite similar to that of the BS model, going to zero linearly at $`p_c`$, the acceptation profile for $`M`$ variables has a more complicated behavior. This difference originates maybe from the fact that $`S`$ variables ($`\eta [0.5,1]`$) can turn into $`M`$ ($`\eta [0,0.5]`$) variables during the dynamics of the system, while $`M`$ variables cannot become $`S`$ variables. Moreover, $`S`$ variables can have developed correlations before the transition to $`M`$ variable and this affects the shape of $`a(\eta )`$ for $`\eta [0.5,1]`$. This might account for the linear part of the acceptation profile of $`M`$ variables, around $`\eta =1`$, but the non linear part is more puzzling. The coupling between $`S`$ variables and $`m`$ variables could play a role in this behavior, too, but at the moment we have no clear explanation of it. ¿From the $`a(\eta )`$ we can get a good estimation of the critical thresholds $`\eta _c^S`$ and $`\eta _c^M`$ for different values of $`p`$ (see Table II).

The stationary state is characterized by a constant ratio between $`S`$ sites and $`M`$ sites, that is to say the evolution equation for the densities $`\rho _S`$ and $`\rho _M`$ of sites $`S`$ and $`M`$ respectively, have an attractive fixed point in the stationary state (see Fig. 16a-d), with the asymptotic values of $`\rho _S,\rho _M`$ depending on the parameter $`p`$. One interesting observation is that, even in the case $`p=0`$, that is to say only $`M`$ sites can be selected by the extremal dynamical rule, there is a stationary state for the system with $`\rho _S0`$. This is due to the particular geometry of the lattice, for which the growth of an $`M`$ site implies the creation or annihilation of some $`S`$ sites. In other words, there cannot be surfaces without slopes ($`S`$ sites).

The roughness exponent accounts for scale free spatial fluctuations in the interface profile. In order to characterize the eventual scale free fluctuations in the dynamical evolution of the system at its asymptotic critical state, that is to say long range temporal correlations, we have studied the avalanche distribution. An avalanche is defined as a sequence of causally connected elementary growth events. For the class of models with quenched disorder and an extremal dynamics to which our model belongs, the initiator of a critical, scale invariant, avalanche is identified in the critical state by a site with quenched variable $`\eta _c^M(p)`$ or $`\eta _c^S(p)`$ (respectively for an $`M`$ initiator and for an $`S`$ initiator). The values of $`\eta _c^M`$ and $`\eta _c^S`$ for different values of $`p`$ can be obtained by the asymptotic histogram distributions shown in Fig.s (11-14). In our case there are two classes of variables, the $`S`$ and $`M`$ sites, and two possible initiators for an avalanche. We call the avalanches that start with an $`S`$ site, $`S`$-avalanches, and the avalanches that start with an $`M`$ site, $`M`$-avalanches. An avalanche lasts when a variable which has been updated before the growth of the initiator is selected by the extremal dynamics. The statistics of off-critical avalanches has been shown to have the form :

$$P^X(s;\eta )=s^{\tau _X}f_X(|\eta \eta _c^X|s^{\sigma _X})$$
(5)
where $`X=S,M`$, and $`\eta `$ is the initiator of an $`X`$-avalanche. This distribution becomes a pure power law for $`\eta =\eta _c^X`$. In the limit $`t\mathrm{}`$ the system self-organizes into the critical state $`\eta =\eta _c^X`$, and the (normalized) avalanche size distribution becomes:

$$P^X(s;\eta _c^X)=\frac{s^{\tau ^X}}{_{s=1}^{\mathrm{}}s^{\tau ^X}}$$
(6)
We have performed a set of about $`10^3`$ realizations of size $`L=8192`$, lasting each one $`2\times 10^6`$ time steps, and collected the statistics of $`S`$ and $`M`$ avalanches over the last $`10^6`$ time steps, for $`p=0.02,0.2,0.5`$. These simulations required about $`2`$ months of CPU time on our computers (a network of DEC alpha machines with clocks going from 266MHz to 500MHz), and are at the best of our computation possibilities. To reduce numerical problems connected with the approximation on $`\eta _c^M,\eta _c^S`$, we used an alternative definition of critical avalanches in models with extremal dynamics, which resides on the causal relation between updated sites inside an avalanche (for details on the definition of critical avalanches see ). The results are shown in Figs. (17), (18). Even after this big computational effort, our numerical results are still a bit noisy. In particular the statistic of $`S`$ avalanches for $`p=0.02`$ is really poor. This is due to the fact that for small $`p`$ values most of sites selected by the dynamics are $`M`$ sites. Consequently, it is difficult to observe a quite clear power law behavior for both the $`S`$-avalanches and $`M`$-avalanches distributions. We point out that the presence of long range temporal correlations is not necessary for the model to have self-similar or self-affine spatial properties, as already observed in a different context .

In order to better establish the critical properties of our model we have measured the power spectrum $`S(k)`$ of the equilibrium surface. The model studied here is a discretized cellular automaton which can be thought as a modified version of the Sneppen model for quenched interface growth. The Sneppen model has been shown to be, at least in $`1+1`$ dimensions, in the same universality class of the continuous Kardar Parisi Zhang equation with quenched noise (QKPZ). It is natural, but not necessarily true, to suppose that for our model, too, it is possible to find a formulation as a continuous growth equation. Given a general growth equation for h(x,t) like:

$$\frac{h(x,t)}{t}=A[h(x,t)]+\gamma (x,t)$$
(7)
where $`A[\mathrm{}]`$ is an operator acting on $`h(x,t)`$ and $`\gamma (x,t)`$ is an uncorrelated quenched noise (the “temporal” direction corresponds to the growth direction of the surface), if the operator $`A[\mathrm{}]`$ is linear and local, the equation can be Fourier transformed into

$$i\omega \stackrel{~}{h}(k,\omega )=\stackrel{~}{A}(k)\stackrel{~}{h}(k,\omega )+\stackrel{~}{\gamma }(k,\omega ),$$
(8)
and by introducing the propagator $`G(k,\omega )`$,

$$\stackrel{~}{h}(k,\omega )=G(k,\omega )\stackrel{~}{\gamma }(k,\omega ).$$
(9)
where $`G(k,\omega )=[i\omega \stackrel{~}{A}(k)]^1`$.

The propagator $`G(k,\omega )`$ of eq. 7 is related to the power spectrum $`S(k)`$ of the interface in the asymptotic state. The power spectrum is so defined:

$$S(k)=FT[h(x)h(x^{})]=|\stackrel{~}{h}(k)|^2$$
(10)
where $`FT[\mathrm{}]`$ is the Fourier transform operator, the average is over different realizations of the noise, $`h(x)=h(x,t=\mathrm{})`$, $`\stackrel{~}{h}(k)`$ is the Fourier transform of $`h(x)`$. Eq. 10 is valid is the case the noise is uncorrelated in space and time, which is the case of our model. The relation between $`\stackrel{~}{G}(k,\omega )`$ and $`S(k)`$ is the following ():

$$\stackrel{~}{G}(k,\omega =0)^2=S(k)$$
(11)
Equation (11) tells us that the power spectrum of the interface can give informations on $`k`$-dependent part of the propagator $`G(k,\omega =0)`$ and consequently on the structure of the operator $`\stackrel{~}{A}(k)`$ in Eq. 8. For self-affine surfaces, the power spectrum follows a power law scaling

$$S(k)k^{2\delta }$$
(12)
where $`\delta `$ is related to the global roughness exponent $`\chi _{glob}`$ through the scaling relation $`2\delta =2\chi _{glob}+1`$ .

Fig.s (19-20) report the behavior of the power spectrum $`S(k)`$ of the interface profile in the critical state for system sizes $`L=2048,8192`$ and different values of $`p`$.

For large values of $`k`$ finite size effects connected to the discretized nature of the model become relevant and there is a deviation from the power law behavior. Away from this saturation effect it is evident in this case how $`S(k)`$ is characterized by a clear cross-over between two power law behaviors. In the low $`k`$ region ($`k<k^{}p/(1/2p)=1/l^{}`$) $`S(k)`$ scales with an exponent $`\delta _{low}`$ close to $`1`$ (actually $`1.07(3)`$) for all values of $`p`$. The QKPZ (Sneppen $`B`$ model) universality class is characterized by a global roughness exponent (coinciding with the local exponent we have measured through $`W^2(l)`$) $`\chi _{glob}=\chi =0.63`$, giving $`\delta =1.13`$, quite near to what we find. For intermediate values of $`k`$ ($`k>k^{}p/(1/2p)=1/l^{}`$) one gets an exponent $`\delta _{mid}`$ quite near to the value $`1.7`$ (we find values between $`1.7`$ and $`1.86`$) which corresponds to the QEW universality class ($`\chi _{glob}1.2`$ ), independently of the value of $`p`$. We point out that the QEW model is super rough with a global roughness exponent $`\chi _{glob}1.2`$, and a local exponent $`\chi =1.0`$ (the one we measured through $`W^2(l)`$). For $`p`$ near $`0.5`$ the intermediate $`k`$ region is very small, and it is difficult to distinguish it from the saturation region. The same happens when one tries to fit the low $`k`$ region for $`p`$ close to zero. If $`p=0.0`$ or $`=0.5`$, no cross-over effect is observed, since these two values of $`p`$ correspond to the “pure” QEW and KPZ universality classes, respectively.

From these results we have a confirmation that the model does not exhibit non-universal critical properties. The apparent non-universal roughness exponent is the consequence of a cross-over effect, tuned by the parameter $`p`$. The fact that this cross-over effect is difficult to observe when studying the scaling of the mean square width $`W^2(l)`$ of the sample, could explain the discrepancy between the experimental findings available in the literature.

## IV Conclusions

In this paper we have introduced a model for surface roughening whose main peculiarity is that of taking explicitly into account the anisotropy of the growth process by means of a tunable phenomenological parameter $`p`$ which introduces local, i.e. dependent on the local environment, dynamical rules in the growth. The simple introduction of just one anisotropy parameter $`p`$ is far from being able to capture all the characteristics of etching processes, and in general of surface roughening experiments. In etching experiments, for example, transport phenomena in the solution are likely to be important and both concentration and agitation have strong effects on transport. Nevertheless, our model captures at least some basic elements of the relationship between anisotropy and the apparent non-universality observed experimentally in etching processes. Moreover, the general requirement of a microscopic dynamical rule depending on the local environment could be a key element in the apparently observed non-universality in kinetic roughening phenomena.

As a main outcome, the model exhibits a cross-over behavior in its critical properties. For each value of the anisotropy factor $`p`$ the system reaches a critical stationary state, with a characteristic length separating a KPZ-like (Sneppen B model) behavior from a QEW-like (Sneppen A model) behavior. The cross-over from one scaling behavior to the other is tuned by the anisotropy parameter $`p`$. If one looks at the scaling of $`W^2(l)`$, the cross-over effect cannot be easily discovered, and the system seems to have a non-universal roughness exponent. A careful analysis of the power spectrum $`S(k)`$, however, shows a clear cross-over effect. These results can probably help to explain the relevant discrepancies among experimental results . We believe that this behavior is the outcome of the complex interplay between the global dynamics which selects at each time step the weakest site and the anisotropy effect which takes into account local constraints in the growth.

It is worthwhile to stress how our model suggests the possibility of several analytical approach, from the treatment of the problem in terms of a continuous stochastic dynamical equation, to the single site mean-field approach , or to the application of a method recently proposed for dynamical models driven by an extremal dynamics . Particularly promising, in this respect, is a recently proposed non-perturbative Renormalization Group approach  which allows one to study self-affine problems.

The author acknowledges financial support under the European network project FMRXCT980183.


text:
# Note on Extended Coherent Operators and Some Basic Properties

## 1 Introduction

This is a continuation of the paper . For the aim of this letter see .

We extend coherent operators and study some basic properties such as

1. the disentangling formula,
2. the resolution of unity,
3. the commutation relation,
4. the matrix elements,
5. the trace,
6. the Glauber formula.

These properties are very well known for coherent operators, see . We give explicit forms to (a) $``$ (f) for extended version of coherent operators (extended coherent operators in our terminology $`\mathrm{}`$ easy and bad naming !).

Since our extension is in a certain sense very natural, our results must have been known. But the author could not find such references in spite of his efforts. At any rate let us list our results.

The details and further developments of this letter will be included in .

## 2 Coherent Operators and Basic Properties

We make a brief review of some basic properties of coherent operators within our necessity.

Let $`a(a^{})`$ be the annihilation (creation) operator of the harmonic oscillator. If we set $`Na^{}a`$ (: number operator), then we have

$$[N,a^{}]=a^{},[N,a]=a,[a^{},a]=\mathrm{𝟏}.$$
(1)
Let $``$ be a Fock space generated by $`a`$ and $`a^{}`$, and $`\{|n|n𝐍\{0\}\}`$ be its basis. The actions of $`a`$ and $`a^{}`$ on $``$ are given by

$$a|n=\sqrt{n}|n1,a^{}|n=\sqrt{n+1}|n+1,N|n=n|n$$
(2)
where $`|0`$ is a normalized vacuum ($`a|0=0\mathrm{and}0|0=1`$). From (2) state $`|n`$ for $`n1`$ are given by

$$|n=\frac{(a^{})^n}{\sqrt{n!}}|0.$$
(3)
These states satisfy the orthogonality and completeness conditions

$$m|n=\delta _{mn},\underset{n=0}{\overset{\mathrm{}}{}}|nn|=\mathrm{𝟏}.$$
(4)
Definition The coherent operators and coherent states are defined as

$`U(z)=\text{e}^{za^{}\overline{z}a}\text{for}z𝐂,`$ (5)
$`|z=U(z)|0\text{for}z𝐂.`$ (6)
Let us list the three basic properties of coherent operators and coherent states, see  and  for the proofs.

(a) Disentangling Formula We have

$$U(z)=\text{e}^{|z|^2/2}\text{e}^{za^{}}\text{e}^{\overline{z}a}=\text{e}^{|z|^2/2}\text{e}^{\overline{z}a}\text{e}^{za^{}}.$$
(7)
In the proof we use the elementary Baker-Campbell-Hausdorff formula

$$\text{e}^{A+B}=\text{e}^{\frac{1}{2}[A,B]}\text{e}^A\text{e}^B$$
(8)
whenever $`[A,[A,B]]=[B,[A,B]]=0`$, see . This is the key formula in the theory of coherent states.

The coherent states are usually defined as eigenvectors of annihilation operator $`a|z=z|z\text{for}z𝐂`$. The important feature of coherent states is the following resolution (partition) of unity.

(b) Resolution of Unity We have

$$_𝐂\frac{[d^2z]}{\pi }|zz|=\underset{n=0}{\overset{\mathrm{}}{}}|nn|=\mathrm{𝟏},$$
(9)
where we have set $`[d^2z]=d(\text{Re}z)d(\text{Im}z)`$ for simplicity.

(c) Commutation Relation We have

$$U(z)U(w)=\text{e}^{z\overline{w}\overline{z}w}U(w)U(z)$$
(10)
for $`z,w𝐂`$.

(d) Matrix Elements The matrix elements of $`U(z)`$ are

$`(\text{i})nmn|U(z)|m=\text{e}^{\frac{1}{2}|z|^2}\sqrt{{\displaystyle \frac{n!}{m!}}}(\overline{z})^{mn}L_{n}^{}{}_{}{}^{(mn)}(|z|^2),`$ (11)
$`(\text{ii})nmn|U(z)|m=\text{e}^{\frac{1}{2}|z|^2}\sqrt{{\displaystyle \frac{m!}{n!}}}z^{nm}L_{m}^{}{}_{}{}^{(nm)}(|z|^2),`$ (12)
where $`L_{n}^{}{}_{}{}^{(\alpha )}`$ is the associated Laguerre’s polynomial defined by

$$L_{n}^{}{}_{}{}^{(\alpha )}(x)=\underset{j=0}{\overset{n}{}}(1)^j\left(\genfrac{}{}{0pt}{}{n+\alpha }{nj}\right)\frac{x^j}{j!}.$$
(13)
(e) Trace We have

$$\text{Tr}U(z)=\pi \delta ^2(z)\pi \delta (x)\delta (y)\text{if}z=x+iy.$$
(14)
(f) Glauber Formula Let $`A`$ be any observable. Then we have

$$A=_𝐂\frac{[d^2z]}{\pi }\text{Tr}[AU^{}(z)]U(z)$$
(15)
This formula plays an important role in the field of homodyne tomography.

## 3 Extended Coherent Operators and Basic Properties

We in this section define extended coherent operators and study the basic properties corresponding to ones of coherent operators.

Definition The extended coherent operators and extended coherent states are defined as follows :

$`U(z,t)`$ $`=`$ $`\text{e}^{za^{}\overline{z}a+itN}\text{for}z𝐂,t𝐑,`$ (16)
$`|z,t`$ $`=`$ $`U(z,t)|0\text{for}z𝐂,t𝐑.`$ (17)
In this definition at first sight we would want to transform (16) as

$$za^{}\overline{z}a+itN=it\left(a+\frac{z}{it}\right)^{}\left(a+\frac{z}{it}\right)\frac{i|z|^2}{t}.$$
(18)
But with this form it is impossible to take a limit $`t0`$, so we don’t use this one in this paper.

Let us list the corresponding properties of extended coherent operators (16) and extended coherent states (17). Before stating our result let us prepare some notations.

$`f(t)`$ $`=`$ $`{\displaystyle \frac{\text{e}^{it}1}{it}},`$ (19)
$`g(t)`$ $`=`$ $`{\displaystyle \frac{\text{e}^{it}(1+it)}{t^2}}.`$ (20)
We here note that

$`|f(t)|`$ $`=`$ $`{\displaystyle \frac{\text{sin}(t/2)}{t/2}},`$ (21)
$`|f(t)|^2`$ $`=`$ $`\left\{g(t)+\overline{g(t)}\right\}.`$ (22)
Let us state our result :

(a) Disentangling Formula We have

$`U(z,t)`$ $`=`$ $`\text{e}^{g(t)|z|^2}\text{e}^{f(t)za^{}}\text{e}^{itN}\text{e}^{f(t)\overline{z}a},\text{or}`$ (23)
$`=`$ $`\text{e}^{\overline{g(t)}|z|^2}\text{e}^{\overline{f(t)}\overline{z}a}\text{e}^{itN}\text{e}^{\overline{f(t)}za^{}}.`$ (24)
This becomes the key formula in the theory of extended coherent states.

(b) Resolution of Unity We have

$$_𝐂\frac{|f(t)|^2[d^2z]}{\pi }|z,tz,t|=\underset{n=0}{\overset{\mathrm{}}{}}|nn|=\mathrm{𝟏}.$$
(25)
We note that the measure to satisfy resolution of unity is not unique, so the following measure is not so bad :

$$_𝐑\frac{|f(t)|^2\text{e}^{|t|}}{2}𝑑t_𝐂\frac{[d^2z]}{\pi }|z,tz,t|=\underset{n=0}{\overset{\mathrm{}}{}}|nn|=\mathrm{𝟏}.$$
(26)
(c) Commutation Relation We have

$$U(z,t)U(w,s)=\text{e}^{\left\{\overline{f(t)}\overline{f(s)}z\overline{w}f(t)f(s)\overline{z}w\right\}}U(w\text{e}^{it},s)U(z\text{e}^{is},t)$$
(27)
for $`z,w𝐂\text{and}t,s𝐑`$.

Let us here make a change of variables $`zw=f(t)z`$. Then matrix elements of $`U(z,t)`$ are written down by those of $`U(w)`$.

(d) Matrix Elements The matrix elements of $`U(z,t)`$ are

$$n|U(z,t)|m=\text{e}^{\left\{\frac{1}{2}+\frac{\overline{g(t)}}{|f(t)|^2}\right\}|w|^2+itm}n|U(w)|m.$$
(28)
(e) Trace For $`U(z,t)`$ we have an another decomposition :

$$U(z,t)=\text{e}^{\frac{i}{t}|z|^2}\text{e}^{\frac{i}{t}(za^{}+\overline{z}a)}\text{e}^{itN}\text{e}^{\frac{i}{t}(za^{}+\overline{z}a)},$$
(29)
so we have

$$\text{Tr}U(z,t)=\text{e}^{\frac{i}{t}|z|^2}\text{Tr}\text{e}^{itN}=\frac{\text{e}^{i|z|^2/t}}{1\text{e}^{it}}(\text{as Abel sum}),$$
(30)
see also (18).

We here note $`\text{Tr}U(z,0)=\text{Tr}U(z)=\pi \delta ^{(2)}(z)`$ , (14). Then we have a natural question :

Question Is the following correct ?

$$\text{lim}_{t0}\text{Tr}U(z,t)=\text{Tr}U(z,0)\text{for}z𝐂.$$
(31)
At first sight the answer of this question seems to be no, but is yes. This has been solved by S. Sakoda.

Making use of (28) and (22) we have the Glauber formula.

(f) Glauber Formula Let $`A`$ be any observable. Then we have

$$A=_𝐂\frac{|f(t)|^2[d^2z]}{\pi }\text{Tr}[AU^{}(z,t)]U(z,t).$$
(32)
A comment is in order. When taking the limit $`t0`$ these properties reduce to corresponding ones in the preceding section because $`f(t)1\text{and}g(t)1/2`$ from (19) and (20).

## 4 Discussion

In this paper we listed the six basic properties of coherent operators and investigated these ones for extended coherent operators. We had a natural extension. The details of several caluculations performed in Section 3 will be published in .

Here let us give a perspective (or application) of our work. We can also extend the squeezed operator in  as follows.

$$V(z,t)=\text{e}^{zK_+\overline{z}K_{}+itK_3}\text{for}z𝐂\text{and}t𝐑$$
(33)
where

$$K_+\frac{1}{2}\left(a^{}\right)^2,K_{}\frac{1}{2}a^2,K_3\frac{1}{2}\left(a^{}a+\frac{1}{2}\right).$$
(34)
We have studied some basic properties of this operator in , so that we can consider the product of two unitary operators

$$U(z,t)V(w,s)\text{for}z,w𝐂\text{and}t,s𝐑.$$
(35)
It may be possible to search the same lines as in  and , . These will be published in an another paper, .

Acknowledgment.
The author wishes to thank S. Sakoda for his helpful comments and suggestions.


text:
# 1 Production of dileptons from a)𝐽/𝜓 and b)𝜂_𝑐, where the blob represents the bound state wave function.

Over the past decade many different experiments in Nucleus-Nucleus (A-A) collisions  have consistently shown excess dilepton production with dimuon invariant mass in the range 1.5-2.5 GeV for $`\mu ^+\mu ^{}`$ pairs. These data sets were interpreted by appropriately scaling proton-nucleus (p-A) data. Essentially in all these data sets the dilepton sources are either Drell-Yan pairs or decays of $`J/\psi `$ or $`D\overline{D}`$ . The bulk of the data appears to agree in general by consideration of these sources alone. However there is a significant departure between the observed dilepton pairs and the theoretical estimate based on the above sources, in the intermediate mass range (IMR) (i.e. $`\mu ^+\mu ^{}`$ invariant mass in the range 1.5-2.5 GeV).

In the literature there have been several attempts to explain this discrepancy. Some of these explanations are based on decrease in $`\rho `$ meson mass due to thermal effects in $`e^+e^{}`$ data , D-rescattering , enhanced $`D\overline{D}`$ production, in-flight $`\pi ^+\pi ^{}`$ decaying to $`e^+e^{}`$ , and so on. Many of these physical processes suggested as an explanation are interesting in their own right but in all of these the explanation for excess dileptons in the IMR is at best partial and generally tend to be relevant in a regime different from the IMR. Fireball hydrodynamics with adjustable parameters however seem to contribute in the IMR regime in Pb-Pb at 158 GeV/nucleon central collisions .

Present data for di-leptons from p-A and A-A collisions over the entire kinematic regime upto 5 GeV agrees with the conventional QCD explanation in terms of Drell-Yan process and vector mesons except for the IMR region mentioned above. The overall picture involving charm quarks  is that when there is sufficient energy exchange in a collision, protons have non-negligible charm content ($`c\overline{c}`$ pairs), and substantial high energy gluons which in turn can decay to $`c\overline{c}`$ pairs. These $`c\overline{c}`$ pairs occassionally form bound states such as $`J/\psi `$ by emitting a soft gluon to maintain color balance, or can further polarize $`u\overline{u}`$ or $`d\overline{d}`$ from the surrounding medium to form $`D\overline{D}`$ pairs. Although the present theoretical understanding cannot predict absolute numbers for these processes it is possible to check the consistency of this picture with various p-A and A-A data. By and large the data agrees with various quantitative checks.

This picture however, also suggests that other charm meson bound states such as $`\eta _c,\psi ^{}`$ and $`\chi ^{}`$’s are produced as well. The relative abundance of each of these mesons is expected to be constrained by their sizes. In particular the larger the size of the charm meson the less likely it is to get formed in the hadronic plasma . This expectation is clearly borne out by the data, viz. in any experimental setup $`J/\psi `$ which is a 1S orbital state is produced about 100 times more than $`\psi ^{}`$ which is a 2S orbital state. Both these resonance peaks are clearly visible in the dimuon data.

In this paper we are concerned about the production of the $`\eta _c`$ meson. This is a 1S orbital state and is expected to have the same size as a $`J/\psi `$; furthermore it has almost the same mass. They differ only in their spin and hence it is expected that in any collision where $`c\overline{c}`$ quarks are produced these can form $`\eta _c`$ with about $`1/3`$ probability as $`J/\psi `$. In fact, any suppression mechanisms  due to the hot hadronic medium will equally affect both the mesons. Consequently it is fair to expect that in any of the experimental setups the production cross sections of $`J/\psi `$ and $`\eta _c`$ are similar. The $`\eta _c`$’s once produced will typically decay into lighter hadrons, which makes the direct detection of $`\eta _c`$ virtually impossible. However there is a small (estimable) cross section for it to decay into $`\gamma \mu ^+\mu ^{}`$ or $`\gamma e^+e^{}`$. The relative decay probabilities of $`\eta _c`$ to $`\gamma \mu ^+\mu ^{}`$ and that of $`J/\psi `$ to $`\mu ^+\mu ^{}`$ is essentially determined by the electromagnetic interaction of the charm quark.

The reason for this is clear from Fig. 1 which shows $`J/\psi `$ decay in Fig 1a and that of $`\eta _c`$ in Fig 1b where $`\gamma ^{}`$ in turn decays into a $`\mu ^+\mu ^{}`$ pair. If we restrict ourselves to large invariant mass for the $`\gamma ^{}`$ ( $`>`$ 1 GeV in Fig1b), it is reasonable to expect the charm quark inside the loop to be almost free. In other words, if the $`J/\psi `$ spin dependent wave function is $`\gamma _\mu \varphi `$ and the $`\eta _c`$ wave function is $`\gamma _5\varphi `$ where $`\varphi `$ refers to the remaining spin and spatial part of the wave function, then (ref. Fig 1a. and Fig 1b.) when the mass of $`\gamma ^{}`$ is above 1 GeV, we can take $`\varphi `$ to be well appproximated by free quark and anti-quark propagators upto an overall constant. Hence in the ratio of these processes this constant is irrelevant. In Fig 1b. we integrate the kinematic space for physical $`\gamma `$ and find the $`\eta _c`$ contribution to $`\mu ^+\mu ^{}`$ pair is given by  (for $`M^2M_\eta ^2`$)

$$\frac{dN}{dM}=\frac{1}{3}N_\psi \frac{1}{60}\frac{\alpha }{\pi }\frac{m_c^2}{m_\psi ^2}\left[\frac{M_\eta ^2}{M^2}1\right]\frac{|_0^1\frac{dx}{x}\mathrm{ln}\frac{m_c^2x(1x)M_\eta ^2}{m_c^2x(1x)M^2}|^2}{|_0^1𝑑xx(1x)\mathrm{ln}\frac{m_c^2x(1x)M_\psi ^2}{\mathrm{\Lambda }^2}|^2}$$
(1)
where $`m_c`$ is the mass of the charm quark and $`M_\eta `$ and $`M_\psi `$ are the masses of the $`\eta _c`$ and $`J/\psi `$ mesons. $`\mathrm{\Lambda }`$ is a cut-off parameter which regulates the logarithmic divergence in Fig.1a. Here $`N_\psi `$ is the number of $`J/\psi `$ events per unit mass around the $`J/\psi `$ peak. In reality, in the p-A and A-A experiments, the $`J/\psi `$ peak is broadened due to detector resolution. Consequently, we interpret $`N_\psi `$ as the total number of $`J/\psi `$ events under the broadened peak.

Now we would like to remark on the limitations of this preliminary analysis. In any experimental set up the detection of $`\mu ^+\mu ^{}`$ is limited by various cuts - in particular, rapidity and Collins-Soper  angle cuts. It is with these constraints that the experiment determines $`N_\psi `$, the number of $`J/\psi `$ events seen. $`J/\psi `$ undergoes a two body decay to dimuons and hence will have different acceptance ratios to that of the three body decay of $`\eta _c`$. These differences can be taken into account in detail, by doing a Monte Carlo simulation of the experiment using PYTHIA, for instance. In this preliminary analysis we shall ignore these details. Consequently we are assuming that the acceptance ratios of $`J/\psi `$ and $`\eta _c`$ events are about the same. It is known, for example, that the experimental acceptance ratio is about 15% for $`J/\psi `$ events. For $`\eta _c`$ we expect that when most of the energy is carried away by the dimuons alone and in addition, we also do not detect the real photon, it is possible that the acceptance ratio for $`\eta _c`$ may not be very different. However, for sufficiently small invariant mass ($`<`$ 1 GeV) this assumption will fail. In spite of ignoring these details, we feel that the qualitative features of our analysis can be of importance in understanding the excess dimuon production in the IMR.

The data on S-U collisions at 200 GeV/nucleon  is one of the best studied experimentally in terms of statistics and comparisons between central vs peripheral collisions. Furthermore, most of the features of the earlier experiments with different nuclei are essentially contained, with better statistics, in this data. We therefore take the S-U data as the typical example for doing our analysis and take $`\mathrm{\Lambda }`$ = .85 GeV to reproduce the excess dimuon events in the central S-U collisions as shown in Fig.2a. Using the same value for $`\mathrm{\Lambda }`$ we then check with the peripheral collision data, as shown in Fig.2b. Both appear to be reasonably satisfactory. It shows, for example, that in the IMR nearly 10 - 14 % of events (depending on peripheral or central), appear to come from $`\eta _c`$. Taking $`m_c=1.6`$ GeV, it turns out that the logarithm in the denominator in Eq.(1) is very sensitive to the choice of $`\mathrm{\Lambda }`$. This is only to be expected since the $`J/\psi `$ mass is very close to $`c\overline{c}`$ threshold. We further notice that after choosing $`\mathrm{\Lambda }`$, the exact choice of $`m_c`$ in the range 1.1 to 1.8 GeV is insignificant to our estimates.

For the chosen values of $`\mathrm{\Lambda }`$ and $`m_c`$, we can check the p-A data against the corresponding $`N_\psi `$. We find that in these data sets the $`\eta _c`$ contribution, although substantial in IMR, is still less than the background as estimated by the experimental groups . In fact, both in p-A and A-A the $`\eta _c`$ contribution is less than the $`D\overline{D}`$ contribution (see Fig. 2). From  it is clear that in the p-A data the $`D\overline{D}`$ contribution is an order of magnitude lower than the combinatorial background. Hence the $`\eta _c`$ contribution is also much below the background.

Next, data for Pb-Pb at 158 GeV/nucleon data is fitted with the same parameters $`m_c`$ and $`\mathrm{\Lambda }`$ obtained earlier from the S-U data. $`N_\psi `$ here is taken from the Pb-Pb data (for central or peripheral collisions as the case may be). This is shown in Fig. 3a and 3b (central and peripheral). Peripheral data is reasonably well accounted for with our $`\eta _c`$ production mechanism. (This peripheral data can be sensitive to the experimental cuts and therefore this agreement should be re-examined in a more detailed analysis incorporating all the experimental cuts). However in the central collision data Fig.3a the discrepancy between our explanation and the data is still about 40% in the IMR. This discrepancy in the total number of events is rather less sensitive to the experimental cuts. Consequently it is fair to infer that our explanation in terms of $`\eta _c`$ production for central collisions of Pb-Pb at 158 GeV/nucleon is not complete.

For $`e^+e^{}`$ data we have not redone this analysis since the corresponding $`N_\psi `$ is not quoted by the experimental group. It is clear from the analysis however, that there can be substantial contribution here as well from $`\eta _c`$.

The net result of this analysis is that the $`\eta _c`$ contribution to IMR of dimuon events is substantial and can fit the data satisfactorily in all A-A data excepting central Pb-Pb at 158 GeV/nucleon. There is substantial other evidence which suggests that in this data the hadronic plasma may have undergone QGP transition; consequently the hydrodynamic evolution parameters can play an important role .

Finally we would like to speculate over what other candidates from QCD phenomenology may play a role in the dilepton data. It has long been suspected that QCD has glueballs, although experimentally there is as yet no strong evidence \- only some candidate events. In the context of A-A collisions, we can expect that a gluon-rich medium either hadronic or QGP in nature can produce such particles. These in turn can form glueballs of various spin, which again decay into standard mesons and baryons. This makes the detection of glueballs extremely tricky. Occassionally due to electromagnetic interaction through quark loops they do decay into photons and lepton pairs. Estimates of these processes is essentiallly hampered by our lack of understanding of glueball production. Noting that most of A-A collision data other than Pb-Pb are well explained by our $`\eta _c`$ scenario, it would not be out of place to surmise that perhaps some glueballs may be produced after QGP in Pb-Pb collision. Let us now expand on this scenario.

From Fig.3a, after accounting for $`\eta _c`$, the major discrepancy in the data occurs in the narrow region 1.7 - 2.3 GeV. Scalar or pseudoscalar glueballs, lighter than 2 GeV cannot contribute to the above region as their decay goes through a $`\gamma ^{}\gamma `$ process. Thus only heavier than 3 GeV glueballs will contribute over the IMR. Since our present theoretical prejudices coming from either lattice  or sum rule techniques  suggest that the lightest glueballs are around 1.5 - 2 GeV, these light scalar or pseudoscalars cannot contribute to the discrepancy in Fig.3a. On the other hand the range 1.7- 2.3 GeV is reasonably narrow, which suggests that perhaps this is due to the $`1^{}`$ vector glueball with a mass of about 2 GeV decaying to dimuons. The production rate of this vector glueball is not estimable. If we presume that all the discrepancy is to be accomodated by the $`1^{}`$ vector glueball then we can infer from Fig. 3a that about $`34\times 10^3`$ dimuon events are due to this glueball.

Finally a word of caution. It is clear that in central Pb-Pb at 158 GeV/nucleon some interesting different physics such as the hydrodynamic fireball expansion model as envisaged by  may be at work. In such a case our association of all the extra events to glueballs alone would be incorrect. It would be worthwhile therefore if there were other corroborative evidence in the data to support glueball production.

Keeping this in mind, we can speculate that in RHIC, where the QGP phase is produced and lasts longer, the effects of glueballs would be enhanced. The net effect of this in the dimuon spectrum in the IMR would be an enhancement of the bump in that region, thereby increasing the discrepancy between the data and present explanations by much more than 40%.

To summarize, the standard QCD explanation for the dimuon spectrum naturally involves the $`\eta _c`$ contribution as well. Interestingly, this contribution can be self-consistently estimated and the inclusion of this effect explains the various experimental data. However, central Pb-Pb data which is believed to undergo perhaps a QGP phase transition continues to show a discrepancy which perhaps needs explaining through a qualitatively new mechanism.

Acknowledgements: We are grateful to E. Scomparin for sending us the data sets for S-U and Pb-Pb scattering and for useful comments. We would like to thank the anonymous referee for clarifying certain matters and thereby helping us present the paper with more clarity.


text:
# 1 Introduction

## 1 Introduction

In recent years, reaction–diffusion systems have been studied by many people. As mean–field techniques, generally, do not give correct results for low–dimensional systems, people are motivated to study exactly–solvable stochastic models in low dimensions. Moreover, solving one–dimensional systems should in principle be easier. Exact results for some models on a one–dimensional lattice have been obtained, for example in \[1–12\]. Different methods have been used to study these models, including analytical and asymptotic methods, mean field methods, and large–scale numerical methods.

Some interesting problems in non–equilibrium systems are non–equilibrium phase transitions described by phenomenological rate equations, and the way the system relaxes to its steady state. Kinetic generalizations of the Ising model, for example the Glauber model or the Kawasaki model, are such phenomenological models and have been studied extensively \[13–18\]. Combination of the Glauber and the Kawasaki dynamics has been also considered \[19–21\].

In this article, we want to study an asymmetric generalization of the zero–temperature Glauber model on a lattice with boundaries. There are also sources (or sinks) of particles at the end points of the lattice. We study the dynamics of the particle density, and specially the large time behavior of the system. In the thermodynamic limit, the system shows two kinds of phase transition. One of these is a static phase transition, the other a dynamic one. The static phase transition is controlled by the reaction rates, and is a discontinuous change of the behavior of the derivative of the stationary particle density at the end points, with respect to the reaction rates. The dynamic phase transition is controlled by the injection- and extraction- rates of the particles at the end points, and is a discontinuous change of the relaxation time towards the stationary configuration.

## 2 Asymmetric Glauber model at zero temperature

In the ordinary Glauber model, the interaction is between three neighboring sites. Spin flip brings the system to equilibrium with a heat bath at temperature $`T`$. A spin is flipped with the rate $`\mu :=1\mathrm{tanh}\frac{J}{kT}`$ if the spin of both of its neighboring sites are the same as itself, and is flipped with the rate $`\lambda :=1+\mathrm{tanh}\frac{J}{kT}`$ if the spin of both of its neighboring sites are opposite to it. At domain boundaries, the spins are flipped with unit rate. So the interactions can be written as,

$`AAAA\mathrm{}A`$ and $`\mathrm{}\mathrm{}\mathrm{}\mathrm{}A\mathrm{}\mu `$
$`A\mathrm{}AAAA`$ and $`\mathrm{}A\mathrm{}\mathrm{}\mathrm{}\mathrm{}\lambda `$
$`AA\mathrm{}A\mathrm{}\mathrm{}`$ and $`\mathrm{}\mathrm{}A\mathrm{}AA1`$
where spin up and spin down are denoted by $`A`$ and $`\mathrm{}`$. One can interpret an up spin as a particle, and a down spin as a hole. At zero temperature, the Glauber dynamics is effectively a two–site interaction :

$$A\mathrm{}\{\begin{array}{cc}AA\hfill & \\ \mathrm{}\mathrm{}\hfill & \end{array}$$
(1)
$$\mathrm{}A\{\begin{array}{cc}AA\hfill & \\ \mathrm{}\mathrm{}\hfill & \end{array},$$
(2)
where all the above processes occur with the same rate.

One can consider the following interactions, as an asymmetric generalization of the zero–temperature Glauber model.

$$A\mathrm{}\{\begin{array}{cc}AA\hfill & u\hfill \\ \mathrm{}\mathrm{}\hfill & v\hfill \end{array}$$
(3)
$$\mathrm{}A\{\begin{array}{cc}AA\hfill & v\hfill \\ \mathrm{}\mathrm{}\hfill & u\hfill \end{array}.$$
(4)
If $`uv`$, the above system has left–right asymmetry. The above system on an infinite lattice has been investigated in , where its $`n`$–point functions, its equilibrium states, and its relaxation towards these states are studied. It can be easily shown that the time evolution equation for the average densities of the system with the above interactions are the same as that of a system with the following interactions, where diffusion is also present:

$$A\mathrm{}\{\begin{array}{cc}\mathrm{}A\hfill & \mathrm{\Lambda }\hfill \\ AA\hfill & u\mathrm{\Lambda }\hfill \\ \mathrm{}\mathrm{}\hfill & v\mathrm{\Lambda }\hfill \end{array}$$
(5)
$$\mathrm{}A\{\begin{array}{cc}A\mathrm{}\hfill & \mathrm{\Lambda }^{}\hfill \\ AA\hfill & v\mathrm{\Lambda }^{}\hfill \\ \mathrm{}\mathrm{}\hfill & u\mathrm{\Lambda }^{}\hfill \end{array}.$$
(6)
Consider a lattice with $`L`$ sites and an asymmetric zero–temperature Glauber dynamics as the interaction. The rates of injection and extraction of particles in the first site (final site) are $`a`$ and $`a^{}`$ ($`b`$ and $`b^{}`$), respectively. The time evolution equations for the average densities are then

$`\dot{n}_k`$ $`=`$ $`(u+v)n_k+un_{k1}+vn_{k+1}\mathrm{for}k1,L`$ (7)
$`\dot{n}_1`$ $`=`$ $`a(a+a^{}+v)n_1+vn_2`$ (8)
$`\dot{n}_L`$ $`=`$ $`b(b+b^{}+u)n_L+un_{L1}.`$ (9)
First, let us calculate the profile of average densities at large times. At large times, the system goes to its stationary state, and the time–derivatives of the left–hand sides vanish. One can see then that the solution to the above system is

$$n_k(\mathrm{})=\alpha +\beta \left(\frac{u}{v}\right)^k.$$
(10)
Putting this in (7), one can easily find $`\alpha `$ and $`\beta `$. In the thermodynamic limit ($`L\mathrm{}`$), the solution becomes

$`n_k(\mathrm{})`$ $`=`$ $`{\displaystyle \frac{a}{a+a^{}}}+{\displaystyle \frac{ba^{}ab^{}}{(a+a^{})(b+b^{}+uv)}}\left({\displaystyle \frac{u}{v}}\right)^{kL},u>v`$ (11)
$`n_k(\mathrm{})`$ $`=`$ $`{\displaystyle \frac{b}{b+b^{}}}+{\displaystyle \frac{ab^{}a^{}b}{(b+b^{})(a+a^{}+vu)}}\left({\displaystyle \frac{u}{v}}\right)^{k1},u<v.`$ (12)
It is seen that the density profile is flat at the left end ($`kL`$) for $`u>v`$ and its value is independent of the reaction rates. But as $`v`$ exceeds $`u`$, the density profile acquires a finite slope, proportional to $`\mathrm{ln}(u/v)`$. This is the static phase transition previously mentioned.

Now return to the dynamics of the system. Let us write the homogeneous part of (7) as

$$\dot{n}_k=h_k^ln_l,$$
(13)
and find the eigenvalues and eigenvectors of the operator $`h`$. One finds

$`Ex_k`$ $`=`$ $`(u+v)x_k+ux_{k1}+vx_{k+1},k1,L`$ (14)
$`Ex_1`$ $`=`$ $`(a+a^{}+v)x_1+vx_2,`$ (15)
$`Ex_L`$ $`=`$ $`(b+b^{}+u)x_L+ux_{L1},`$ (16)
where the eigenvalue and eigenvector have been denoted with $`E`$ and $`x`$, respectively. The solution to these equations is

$$x_k=\alpha z_1^k+\beta z_2^k,$$
(17)
where $`z_i`$’s satisfy

$$E=(u+v)+vz+\frac{u}{z},$$
(18)
and

$`v(\alpha z_1^2+\beta z_2^2)(E+a+a^{}+v)(\alpha z_1+\beta z_2)`$ $`=`$ $`0`$ (19)
$`u(\alpha z_1^{L1}+\beta z_2^{L1})(E+b+b^{}+u)(\alpha z_1^L+\beta z_2^L)`$ $`=`$ $`0.`$ (20)
To have nonzero solutions for $`x`$, these last two equations should be dependent, the criterion for which is

$$(u+z_1\delta a)(vz_2^{L+1}+z_2^L\delta b)(u+z_2\delta a)(vz_1^{L+1}+z_1^L\delta b)=0,$$
(21)
where (18) has been used to eliminate $`E`$, and $`\delta a:=a+a^{}u`$ and $`\delta b:=b+b^{}v`$. Defining

$`Z_i`$ $`:=`$ $`z_i\sqrt{{\displaystyle \frac{v}{u}}}`$ (22)
$`A`$ $`:=`$ $`{\displaystyle \frac{\delta a}{\sqrt{uv}}}`$ (23)
$`B`$ $`:=`$ $`{\displaystyle \frac{\delta b}{\sqrt{uv}}},`$ (24)
(21) is simplified to

$$Z^{(L+1)}(1+AZ)(1+BZ)Z^{L+1}(1+AZ^1)(1+BZ^1)=0.$$
(25)
The eigenvalue $`E`$ satisfies

$$E=(u+v)+\sqrt{uv}(Z+Z^1).$$
(26)
Two obvious solutions of the equation (25) are $`Z=\pm 1`$. But these generally don’t correspond to eigenvalues and eigenvectors. In fact for these solutions, $`Z`$ and $`Z^1`$ are the same, so that (17) should be modified to

$$x_k=(\alpha +\beta k)(\pm 1)^k,$$
(27)
and it is not difficult to see that these do not fulfill the boundary conditions unless $`\alpha =\beta =0`$. Equation (26) can be written as a polynomial equation of order $`2(L+1)`$, and hence has $`2L`$ more roots in addition to $`Z=\pm 1`$. For these $`2L`$ roots, if $`Z`$ is a root $`Z^1`$ is another root, and these two correspond to one eigenvalue and one eigenvector. So the $`L\times L`$ matrix $`h`$ does have $`L`$ eigenvalues and eigenvectors.

For $`A=B=0`$, (25) is very simple and its nontrivial solutions are

$$Z_{(s)}=e^{i\pi s/(L+1)},1sL,$$
(28)
and their inverses. All of these are phases and the real–part of the corresponding eigenvalues satisfy

$$\mathrm{Re}(E)(u+v)+2\sqrt{uv}\mathrm{cos}\left(\frac{\pi }{L+1}\right)<0.$$
(29)
The maximum of the real–part of the eigenvalues determines the relaxation time toward the stationary average–density profile. That is

$$\tau =\left[u+v2\sqrt{uv}\mathrm{cos}\left(\frac{\pi }{L+1}\right)\right]^1.$$
(30)
In the limit $`L\mathrm{}`$, this is simplified to

$$\tau =(u+v2\sqrt{uv})^1.$$
(31)
The general solution to (7) is seen to be

$$n_k(t)=\underset{s,m}{}\frac{2}{L+1}\mathrm{exp}[E_{(s)}t]n_k(0)\left(\frac{u}{v}\right)^{(km)/2}\mathrm{sin}\left(\frac{\pi sk}{L+1}\right)\mathrm{sin}\left(\frac{\pi sm}{L+1}\right).$$
(32)
Now consider the general case. The equation (25) can be written as

$$G(Z):=F(Z)F(Z^1)=0,$$
(33)
where

$$F(Z):=Z^{(L+1)}(1+AZ)(1+BZ).$$
(34)
If $`Z`$ is a phase, it satisfies (25) provided $`F(Z)`$ is real. Consider the phase of $`F(Z)`$ for unimodular $`Z`$. We have

$$\varphi [F(z)]=(L+1)\varphi (Z)+\varphi (1+AZ)+\varphi (1+BZ),$$
(35)
where $`\varphi `$ denotes the phase of its argument. As the phase of $`Z`$ is changed from $`0`$ to $`2\pi `$, the change of the phase of $`F(z)`$ is

$$\mathrm{\Delta }\varphi [F(Z)]=2\pi (L+1)+\mathrm{\Delta }\varphi (1+AZ)+\mathrm{\Delta }\varphi (1+BZ).$$
(36)
For $`|A|,|B|<1`$, the phase changes of $`(1+AZ)`$ and $`\varphi (1+BZ)`$ is zero, as $`Z`$ moves on the whole unit circle. So,

$$\mathrm{\Delta }\varphi [F(Z)]=2\pi (L+1),\text{for}|A|,|B|<1.$$
(37)
But this means that the phase of $`F(Z)`$ will be an integer multiple of $`\pi `$ for at least $`2(L+1)`$ points on the unit circle. So all of the solutions of (25) are still phases, although they may be not uniformly spaced on the unit circle. Two of these are $`Z=\pm 1`$. The remaining $`2L`$ points correspond to $`L`$ eigenvalues and eigenvectors for $`h`$. One concludes that for $`|A|,|B|<1`$,

$$\tau =(u+v2\sqrt{uv}\mathrm{cos}\theta )^1,$$
(38)
for some $`\theta `$. Specially, at $`L\mathrm{}`$, the relaxation time is the same as the relaxation time for $`A=B=0`$, that is the same as (31).

If, for example, $`|A|>1`$, then the total phase change of $`(1+AZ)`$ is no longer zero. It is $`2\pi `$. One may then lose two of the roots of the unit circle. Note that the mere fact that $`\mathrm{\Delta }\varphi [F(Z)]=2\pi L`$ does not mean that there are just $`2L`$ solutions of (25) on the unit circle, since the phase of $`F(Z)`$ needn’t be monotonic. To find values of $`A`$ and $`B`$, for them the number of the solutions of (25) on the unit circle is $`2L`$ or $`2(L1)`$, consider the function $`G(Z)`$ at the points $`Z=\pm 1`$. Increasing $`|A|`$ or $`|B|`$, two of the roots on the unit circle tend to $`1`$ or $`1`$, and then move out of the unit circle and on the real line. At the point that this occurs, either $`G^{}(1)`$ or $`G^{}(1)`$ become zero, as there will be multiple roots at $`\pm 1`$. So the criterion for each change (losing $`2`$ roots of the unit circle) is that either $`G^{}(1)`$ or $`G^{}(1)`$ become zero. The curves in the $`AB`$ plane, corresponding to these changes are

$$L(1+A)(1+B)+1AB=0,G^{}(1)=0,$$
(39)
and

$$L(1A)(1B)+1AB=0,G^{}(1)=0.$$
(40)
These curves divide the plane into six regions:

* all of the solutions are phases.
* $`2L`$ phase solutions, 2 real negative solutions.
* $`2(L1)`$ phase solutions, 4 real negative solutions.
* $`2L`$ phase solutions, 2 real positive solutions.
* $`2(L1)`$ phase solutions, 4 real positive solutions.
* $`2(L1)`$ phase solutions, 2 real negative solutions, 2 real positive solutions.

In the above, by real solutions it is meant real solutions besides the trivial solutions $`\pm 1`$. Note, however, that not all of this plane is physical. The physical region is that part of the plane, which corresponds to nonnegative values for the injection and extraction rates. Returning to the definitions of $`A`$ and $`B`$, equation (22), it is seen that

$`A_{\mathrm{min}}`$ $`=`$ $`\sqrt{{\displaystyle \frac{u}{v}}}`$ (41)
$`B_{\mathrm{min}}`$ $`=`$ $`\sqrt{{\displaystyle \frac{v}{u}}}.`$ (42)
The point $`(A_{\mathrm{min}},B_{\mathrm{min}})`$ itself is on the curve

$$A_{\mathrm{min}}B_{\mathrm{min}}=1.$$
(43)
One can see that, unless $`u=v`$, part of the physical region is in the region IV, where two of the solutions of (25) are real positive (and of course inverse of each other). This makes the maximum of the real part of $`E`$ larger than $`(u+v)+2\sqrt{uv}`$, and correspondingly the relaxation time larger than (31).

In the thermodynamic limit $`L\mathrm{}`$, the regions of the $`AB`$ plane are greatly simplified. In fact the curves corresponding to $`G^{}(\pm 1)=0`$ become

$$A=1,\text{or}B=1,G^{}(1)=0,$$
(44)
and

$$A=1,\text{or}B=1,G^{}(1)=0.$$
(45)
In the case either $`|A|`$ or $`|B|`$ are greater than 1, the real roots of (25) are simply

$$Z=A,A^1,|A|>1$$
(46)
and

$$Z=B,B^1,|B|>1$$
(47)
Then, if for example $`A`$ is negative and less than $`1`$ and $`A<B`$, the maximum real part of $`E`$ is $`(u+v)\sqrt{uv}(A+A^1)`$, and the relaxation time of the system is

$$\tau =[u+v+\sqrt{uv}(A+A^1)]^1,$$
(48)
which is greater than (31). This is a phase transition which occurs at $`A=1`$. For $`A>1`$, the relaxation time is constant, (31). For $`A<1`$, it is $`A`$–dependent. The minimum of $`A`$ is $`\sqrt{u/v}`$, for which one of the eigenvalues of $`h`$ become zero, and the relaxation time becomes infinite. The same effect is seen for $`B<1`$. As mentioned before, for $`u=v`$ (the ordinary zero–temperature Glauber model) no part of the region IV is in the physical region, and this transition does not occur. The phase transition discussed here, is the dynamical phase transition mentioned before.

Acknowledgement
The authors would like to thank Institute for Studies in Theoretical Physics and Mathematics for partial support.