format: text

text:
# 1.65 𝜇m (H-band) surface photometry of galaxies. IX: photometric and structural properties of galaxies.

## 1 Introduction

The most commonly used galaxy classification system is the one originally developed by Hubble (realm (1936)), and brought to completion by Sandage (atlas (1961)). Among the main reasons for its success are the fact that it is relatively easy to use, and the fact that it has proven successful in ordering the “realm of nebulae” according to important physical characteristics. The correlations between quantities related to a galaxy stellar population (like galaxy color) and Hubble type are well known (see Roberts & Haynes RH94 (1994) for a recent review). However, the Hubble classification system suffers from a number of important drawbacks that are becoming more and more evident as the samples of galaxies studied in detail become larger both at low and high redshift. Among these drawbacks, its qualitative and quantized nature, that prevents fitting relations as a function of type on a solid statistical basis, and the fact that a large (and ever-growing, with the completion of new surveys) fraction of galaxies cannot fit well inside the scheme, are the most limiting features.

These and other limitations of the Hubble system have been recognized for a long time, and a number of extensions to the scheme, or alternative classification schemes, have been proposed through the years. Best-known among such efforts are the extension to the Hubble system devised by de Vaucouleurs (deVa59 (1959)), and the alternative scheme proposed by Morgan (Mor58 (1958), Mor59 (1959)). Still, neither of these schemes is widely used by the astronomical community. More recently, Whitmore (Whit84 (1984)) proposed a sound, quantitative scheme for the classification of spiral galaxies, that is based on two parameters: “scale” and “form” (following Whitmore’s terminology). The scale parameter is a combination of the galaxy absolute blue luminosity and isophotal radius, while the form is a combination of the galaxy color and bulge-to-total light ratio. This idea has been further developed by Bershady et al. (BJC (2000); hereafter BJC), who have presented a classification scheme encompassing all morphological types. In doing so, they were forced to increase the dimensionality of the classification volume from two to three dimensions, that they define to be a spectral index (galaxy color), a scale (galaxy surface brightness), and a form (galaxy light concentration, or asymmetry) parameter. This classification scheme is fully quantitative, can be used to classify galaxies over a broad range of redshifts (provided that images with enough spatial resolution can be obtained), and is quite successful in separating normal galaxies into classes that are largely coincident with those one would obtain using Hubble types.

One important feature for a classification scheme should be its ability to describe not only the “current” status of a galaxy (by “current” we mean the status at the epoch of the observation, the latter being of course dependent on the galaxy redshift), but also its evolutionary history. This goal however is complicated by the degeneracy that many galaxy parameters show with respect to different possible formation and evolution histories. A well-known example is the age-metallicity degeneracy for stellar populations in early-type galaxies (Worthey Wor94 (1994)), but we observe a rather similar color-star formation history degeneracy for late-type, star-forming galaxies (Gavazzi et al. 2001b ). To achieve the goal of obtaining a classification scheme that might provide information on a galaxy’s past evolutionary history we are building a large dataset of multifrequency observations of nearby galaxies, that covers objects of all morphological types over an extended range of luminosities and environmental conditions, from the cores of dense clusters to relatively under-dense regions in the field. We began our effort by focusing on late-type star-forming galaxies, and the results of the analysis of these systems are summarized in Gavazzi et al. (1996c ), and Gavazzi & Scodeggio (GS96 (1996)). Since then we have extended our sample to also cover early-type galaxies, and we have completed a photometric near-infrared ($`H`$-band) imaging survey that covers a large fraction of our total sample, to measure in a completely homogeneous fashion the photometric and structural properties of nearby normal galaxies. Previous papers in this series describe the photometric survey (Gavazzi et al. 1996a  – Paper I; Gavazzi et al. 1996b  – Paper II; Boselli et al. Bos97 (1997); Gavazzi et al. 2000a  – Paper III; Boselli et al. Bos00 (2000) – Paper IV; Gavazzi et al. 2001a  – Paper VII) and provide details on the derivation of the photometric parameters used here (Gavazzi et al. 2000b  – Paper V). In this work we describe the observed correlations between photometric and structural properties for the galaxies in our sample. For the first time this kind of analysis is presented using $`H`$-band near-infrared imaging data. This offers two key advantages with respect to the more often used optical bands imaging: a much smaller influence from dust extinction on the observed parameters, and the capability of tracing the bulk of the luminous matter mass, without being heavily influenced by the small fraction of young stars that often dominate galaxy emission in the blue optical bands. Two other papers in this series discuss the correlation between their photometric, structural, and dynamical properties (Pierini et al. Dan01 (2002)), and the correlation between these properties and the galaxy star formation history (Boselli et al. Bos01 (2001)). A future paper will discuss the practical application of our classification scheme to a sample of more distant galaxies.

## 2 The data

The analysis of photometric and structural properties of normal galaxies we carry out in this paper is based on a nearly complete Near Infrared $`H`$-band survey of an optically selected sample of nearby galaxies. Full details about the sample selection, the observations, and the data reduction procedures can be found in the previous papers of this series (Papers I, II, III, IV, VII), and here only a brief summary is presented.

Our sample includes all galaxies from the CGCG catalogue (Zwicky et al. 1961- (1968)) (m$`{}_{\mathrm{pg}}{}^{}`$ 15.7) that belong to the Coma supercluster region, including the two clusters Coma (A1656) and A1367, and to the clusters A262 and Cancer, and galaxies from the VCC catalogue (Binggeli et al. Bing85 (1985)) (restricted to m$`{}_{\mathrm{pg}}{}^{}`$ 16.0) in the Virgo cluster. Observations were carried out from 1993 to 1997 with the 1.5 m TIRGO and with the Calar Alto 2.2 and 3.5 m telescopes equipped with the NICMOS3 $`256^2`$ pixel arrays cameras ARNICA and MAGIC, respectively, and more recently with the ESO NTT and TNG telescopes, equipped with the SOFI and NICS or ARNICA cameras. In total we have obtained observations for 1302 galaxies, that span all morphological types from E and dE to dIrr. However, because of the selection criterion of our master catalog, extremely low surface brightness galaxies are not included in our sample. Paper VII provides information on the completeness of our observations.

The surface photometry measurements used in this paper are discussed in detail in Paper V. Each galaxy two-dimensional light distribution was fitted with elliptical isophotes to derive a radial surface brightness profile. This profile was fitted with either a de Vaucouleurs $`r^{1/4}`$ or an exponential disk profile, or a combination of two profiles (either a de Vaucouleurs plus exponential disk, or two exponential profiles, the choice being based on a best fit criterion). Among the 1302 profiles, 192 are fitted with a pure de Vaucouleurs $`r^{1/4}`$ law and 369 with a pure exponential one. The remaining profiles require a Bulge+Disk (B+D) decomposition, with 54 objects showing signs of a truncated outer light distribution. We have also tried to fit the global surface brightness profile with a single Sérsic (Sersic (1968)) $`r^{1/n}`$ profile, but we could not reproduce at all the shape of the observed profile for a large fraction of the B+D galaxies, and therefore abandoned this experiment. A total magnitude $`H_\mathrm{T}`$ was derived for each galaxy by extrapolating the suitable profile to infinity, and adding the extrapolated light contribution to the one measured within the last fitted isophote. Empirical effective radius $`r_\mathrm{e}`$ (the radius within which is contained half of the galaxy light) and effective surface brightness $`\mu _\mathrm{e}`$ (the mean surface brightness within the effective radius) were measured on the observed profile, at the half light point, on the basis of the total magnitude determination. Bulge and disk parameters were derived from the fitted profiles, together with the bulge to total (B/T) light flux (the ratio of the bulge flux to the total flux from the galaxy). We also derived for each galaxy an isophotal radius $`r_\mathrm{H}(20.5)`$ determined in the elliptical azimuthally–integrated profiles as the radius at which the surface brightness reaches 20.5 H–mag arcsec-2, and a concentration index (C31), defined as in de Vaucouleurs (deVa77 (1977)) to be the model–independent ratio between the radii that enclose 75% and 25% of the total light $`H_\mathrm{T}`$. The median uncertainty in the determination of the total magnitude is 0.15 mag, while those on the determination of log $`r_\mathrm{e}`$ and $`\mu _\mathrm{e}`$ are 0.05 and 0.16 mag, respectively.

Although we do not use Hubble types as a basis for the classification of galaxies, we have tried to leave the possibility of a comparison with the traditional classification scheme open. For this reason in the following analysis we exclude from our sample a number of galaxies that have very uncertain morphological classification (90), and peculiar galaxies (69). The total sample used here is therefore composed of 1143 galaxies. Of these, 831 have color measurements available, and are used in the definition of our classification scheme (see Sect. 4).

## 3 The concentration index

A concentration index has been used before as a classification tool, most noticeably in the Morgan classification scheme (Morgan Mor58 (1958), Mor59 (1959)). Other classification schemes involving a concentration index have been proposed, usually trying to take advantage of the correlation existing between image concentration and mean surface brightness in galaxies (e.g. Okamura et al. Oka84 (1984); Kent Kent85 (1985); Doi et al. Doi93 (1993); Abraham et al. Abr94 (1994)). Most recently BJC have presented a global classification scheme for galaxies that includes a concentration index as one of the fundamental classification parameters, along with color, surface brightness, and galaxy asymmetry. All their parameters were derived from $`B`$-band imaging data. They find correlations to be present between concentration and either color, surface brightness, or asymmetry, the strongest one being that between concentration and color.

Here we adopt a definition of concentration index which is somewhat different from the one used by BJC. We define, following de Vaucouleurs (deVa77 (1977)), the concentration index C31 as the ratio between the radii that enclose 75% and 25% of the total galaxy light $`H_T`$, whereas BJC take the logarithm of the ratio between two radii (specifically, those that enclose 80% and 20% of the total light). As we discuss in section 4, we consider our “linear ratio” definition better suited to emphasize the changes of fundamental galaxy properties as a function of galaxy luminosity, as discussed already in Gavazzi et al. (1996c ). With our definition of concentration index, a pure de Vaucouleurs profile should have C31 = 7.02, while a pure exponential disk profile should have C31 = 2.80.

It is quite natural to expect a correlation between the value of C31 and other parameters traditionally used for galaxy classification, like Hubble type, or bulge-to-disk / bulge-to-total light ratios, as the relevance of the bulge inside a galaxy is one of the main parameters in the Hubble scheme, and the only reliable one for edge-on spiral galaxies. We show in the left panel of Fig. 1 the distribution of C31 values as a function of Hubble type. Two characteristics are clearly evident: there is a general correlation between C31 and type, as evidenced by the global decline in the median C31 value going from ellipticals to late type spirals (with dE galaxies representing an exception with respect to this trend; see the Figure inset), but there is also a large scatter within a given type, at least up to type Sb. Part of this scatter might result from misclassifications, due to the greater difficulty in classifying galaxies at the distance of the Coma supercluster with respect to those in the Virgo cluster, but a component of this scatter is due to a real dependence of C31 on galaxy luminosity, as discussed in section 4. It is therefore clear that concentration index and Hubble type are not completely equivalent classification criteria.

The same situation is present when considering the B/T light ratio. As shown in the right panel of Fig. 1 there exists also for B/T a global correlation with Hubble type, with a decline in the median value of B/T when going from ellipticals to late type spirals, and a large scatter within each Hubble type. Also in this case there are classification problems and a real dependence of the B/T ratio on galaxy luminosity that contributes to the scatter. Therefore also the B/T ratio is not a classification criterion equivalent to Hubble type, and it is not equivalent either to the concentration index criterion, as shown in Fig. 2, where we see the presence of a rather large scatter in the value of B/T for any given value of C31.

Keeping in mind these differences, and the fact that all three form criteria (Hubble type, concentration index, and B/T ratio) share some common basis in terms of their ability to differentiate among galaxies with different appearance, structure, and stellar populations, we choose to base our “form” classification on the concentration index C31. With respect to Hubble type it has the advantage of being a quantitative index, that can be easily measured for all galaxies independently from the observer’s viewing angle, and is less affected by distance related resolution problems. With respect to the B/T ratio, it has the advantage of being model independent, not requiring the decomposition of a galaxy into a bulge and a disk component, a process which is always affected by considerable uncertainties. Another advantage offered by the concentration index is its continuous distribution of values, as compared to the quantized distribution of Hubble types, and the peculiar distribution of B/T values (composite B+D photometric profiles show a continuous distribution of values, but pure de Vaucouleurs profiles have B/T equal to 1, and pure exponential disk profiles have B/T equal to 0 by definition).

Finally we remark that, at least in the $`H`$-band, where internal dust extinction effects are quite negligible, galaxy viewing angle has little influence on the derivation of C31 values. The distribution of C31 covers the full range of observed values at all galaxy inclinations (as derived from the observed axial ratio), but we observe a moderately significant trend toward higher C31 values at high inclination. The median C31 value for galaxies with inclination below 50 degrees is 2.99, while for galaxies with higher inclination is 3.56 (in both cases the distribution has a standard deviation of approximately 1.5).

## 4 Scale, form, and color

Whitmore (Whit84 (1984)) introduced the first quantitative galaxy classification scheme that incorporates both structural and spectral parameters, introducing the concepts of galaxy “scale” (in his analysis a combination of absolute blue luminosity and isophotal radius) and “form” (in his analysis a combination of $`BH`$ color and B/T light ratio). He concluded that these two parameters would provide an optimal classification scheme for spiral galaxies, to which his work was limited. His choice for the specific combinations that define the form and scale parameters was based on the observation that those quantities show the strongest correlations between themselves among the 30 parameters taken into consideration in his work. More recently BJC have presented a quantitative classification scheme along the same lines, extended to galaxies of all morphological types. However they need three parameters to obtain a satisfactory classification of all galaxies: a spectral index ($`BV`$ color), a scale (surface brightness), and a form (concentration index or image asymmetry) parameter.

### 4.1 The luminosity, color, C31 cube

Independently from BJC, we have defined a very similar framework to describe the main structural and photometric properties of normal galaxies. Our sample in fact does not include a significant fraction of starbursting, interacting, or extremely low surface brightness objects. In agreement with BJC, we conclude that at least three parameters are needed to provide a unique classification for all galaxies. Figs. 3 and 6 show two alternatives definitions for such a classification cube. The first one, presented in Fig. 3 is based on $`H`$-band luminosity, color (either $`BH`$ or $`BV`$, this second not shown in the figure), and concentration index C31. In the second one, presented in Fig. 6, the effective surface brightness $`\mu _\mathrm{e}`$ replaces the $`H`$-band luminosity. To keep the connection with the Hubble types classification open, in both figures different symbols identify galaxies within different ranges of Hubble types.

It is well known that both early- and late-type galaxies obey a color-luminosity relation (see, for example, Visvanathan & Sandage VS77 (1977); Bower et al. BLE92 (1992) for early-type galaxies, and Tully et al. TMA82 (1982); Gavazzi et al. 1996c  for spiral galaxies). Since the two relations are quite different, their origins are also believed to be different. The small color changes of early-type galaxies are most likely produced by metallicity variations within a coeval galaxy population (Arimoto & Yoshii AY87 (1987); see however Worthey et al. Wor96 (1996) for a discussion on the effects of age variations within a metal homogeneous population). The much broader color range covered by late-type galaxies is mainly the result of a broad range of star formation histories (Searle et al. SSB73 (1973); Kennicutt Kenn83 (1983); Sandage San86 (1986)), coupled with an equally broad range of dust extinction values. This color-luminosity correlation is well visible in the upper panel of Fig. 3.

BJC have also already described a correlation between color and $`B`$-band concentration index, although they did not discuss what we consider to be the main feature visible in the lower panels of Figs. 3 and 6, which is the strong non-linearity in the correlation between C31 and either color, $`H`$-band luminosity, or effective surface brightness. It is clearly evident from these figures that there are two main regimes of galaxy properties, that produce two almost orthogonal distributions in the parameter space we are considering. At low luminosity ($`L_\mathrm{H}<10^{10}L_{}`$) all galaxies consistently have very small C31, independently from their color (and from their Hubble type, for that matter), while only at high luminosity ($`L_\mathrm{H}>10^{10}L_{}`$) we observe galaxies with large values of C31. At the same time blue galaxies consistently have low values of C31, and it is only among the red ones that large values of C31 are present. Alternatively, we can say that galaxies with small C31 can have all colors and all luminosities, while those with a large C31 are only found among red luminous objects. We stress here the fact that this qualitative picture is not sensitive to the choice of photometric bands used to derive the color parameter, at least as far as optical or optical-near infrared colors are involved. The finding of a luminosity threshold for the appearance of high C31 values is in agreement with the finding discussed in Papers V and VII, that pure de Vaucouleurs profiles are absent from the classification of photometric profiles below the luminosity $`L_\mathrm{H}10^{10}L_{}`$.

As discussed in Sect. 2 our sample is composed mainly of galaxies from the Coma Supercluster and the Virgo Cluster. The different distance of these two structures translates into different completeness limits when luminosity is taken into consideration. To provide a more reliable indication of the relative density of galaxies in the various parts of our classification cube we reproduce it again in Fig. 4, plotting in this case only Virgo Cluster galaxies. Notice that even in this case, however, for luminosities $`L_\mathrm{H}<10^9L_{}`$ sample incompleteness becomes non-negligible, and the density of galaxies in the plot does not represent any more the true fraction of objects that lay in this part of the galaxies’ parameter space. The most striking difference between Figs. 3 and 4 concerns spiral galaxies. Within the global sample the color-magnitude relation for these objects appears much steeper than in the more complete Virgo sample, purely as a result of sample incompleteness. In fact our sample is an optically selected one, complete down to a certain optical magnitude (as discussed in Sect. 2). Therefore when one considers $`H`$-band data, the magnitude limit for the sample changes with galaxy color, shifting to brighter luminosities as the color changes from blue to redd (a diagonal line in Fig. 3).

To better illustrate the relation between the classification cube presented here and the classic Hubble classification scheme, we present in Fig. 5 the same plots as in Fig. 3, separating the galaxies in different panels according to their Hubble type. It is evident from these figures that the scatter within the C31-Hubble type relation can be partly explained by the luminosity (color) dependence of the concentration index, that is quite evident even when a restricted range of Hubble types is considered.

### 4.2 The $`\mu _\mathrm{e}`$, color, C31 cube

An almost equivalent picture to that presented in Fig. 3 is obtained substituting the effective surface brightness $`\mu _\mathrm{e}`$ to the $`H`$-band luminosity (see Fig. 6). There are nonetheless two noticeable differences between these two data cubes. The first one is the smaller scatter shown by the $`\mu _\mathrm{e}`$-color relation with respect to the luminosity-color one. This results mainly from a lack of low $`\mu _\mathrm{e}`$ ellipticals, while late-type galaxies show similar behaviour in the two diagrams. This effect could however be due to an observational bias, as red objects with $`\mu _\mathrm{e}`$ fainter than 20 H–mag arcsec-2 are very difficult to observe and detect in the near infrared. As a result our sample is strongly incomplete with respect to dwarf-elliptical galaxies. Should the smaller scatter be confirmed even after the completion of our survey, it could probably be explained by the fact that $`\mu _\mathrm{e}`$ is not purely a scale parameter like luminosity, but a hybrid one, mixing both scale and form in its definition. The second difference is the total absence of a correlation between $`\mu _\mathrm{e}`$ and C31 for high C31 objects, compared to a small but significant trend between luminosity and C31 (the mean $`logL_H`$ for objects with $`4<`$C$`{}_{31}{}^{}<5`$ is 10.6, while for objects with C$`{}_{31}{}^{}>9`$ it is 11.0). Another interesting result is the fact that the observed strong correlation between $`\mu _\mathrm{e}`$ and color, which confirms earlier results presented by Gavazzi et al. (1996c ), is in disagreement with the results presented by BJC (see the upper panels of their Fig. 3), that confirm instead the so-called Freeman’s law (Freeman Free70 (1970)). This is yet another indication that the validity of Freeman’s law is limited to one specific photometric band.

It could be argued that this second classification cube is better suited than the one discussed in the previous section to represent normal galaxy properties, both because of the smaller scatter observed in the $`\mu _\mathrm{e}`$-color relation with respect to the $`L_H`$-color one (should this feature be confirmed when observations of our sample will be completed), and because surface brightness is a distance-independent quantity (cosmological dimming correction aside), whereas luminosity is not. However what could be better in terms of a classification scheme might not be equally relevant in building a physical model for galaxy formation and evolution. In that case total near infrared luminosity, which is a good tracer of galaxy mass, should be certainly considered as a more fundamental parameter than surface brightness.

## 5 The Photometric Fundamental Plane

The possibility of replacing $`H`$-band luminosity with effective surface brightness in the classification cube, as discussed above, is due to the presence of a significant correlation between those two quantities. In fact a number of correlations are known to exist, at least for early-type galaxies, that involve luminosity, effective surface brightness, and effective radius. Luminosity correlates with effective radius (Fish Fish64 (1964)), stellar velocity dispersion (Faber & Jackson FJ76 (1976)), and effective surface brightness (Binggeli et al. Bing84 (1984)), while effective radius and effective surface brightness correlate among themselves (Kormendy Kor77 (1977); see Guzmán et al. Guz93 (1993) for a global discussion about these correlations). These relations exhibit a larger scatter than can be accounted for by measurement errors alone, and this lead to the idea that early-type galaxies populate a plane in the 3–parameter space of effective radius, effective surface brightness and stellar velocity dispersion, independently introduced by Djorgovski & Davis (DD87 (1987)) and by Dressler et al. (7S (1987)). Modified versions of this Fundamental Plane relation have since then been proposed, that use color (de Carvalho & Djorgovski deCar89 (1989)) and luminosity (Scodeggio et al. SGH97 (1997)) as a substitute for the stellar velocity dispersion.

A Photometric Fundamental Plane correlating galaxy luminosity, effective radius $`r_\mathrm{e}`$ and effective surface brightness $`\mu _\mathrm{e}`$ still exists111although one must keep in mind that only two of the three parameters involved are statistically independent, as a result of our definition of effective radius and effective surface brightness: given the total flux $`f_\mathrm{t}`$ and the effective radius $`r_\mathrm{e}=r(f_\mathrm{t}/2)`$, we have $`\mu _\mathrm{e}=2.5\mathrm{log}(f_\mathrm{t}/2\pi r_\mathrm{e}^2)+Z_\mathrm{p}`$ even when considering galaxies of all morphological types, as shown in Fig. 7. However, as it could be expected, the picture one obtains when considering galaxies of all types is more complex than the one provided by early-type objects alone. While the shape and scatter of the luminosity-$`\mu _\mathrm{e}`$ relation are not influenced by the sample morphological mix, both the $`r_\mathrm{e}`$-luminosity and the $`r_\mathrm{e}`$-$`\mu _\mathrm{e}`$ relation show a significant segregation of galaxies as a function of their concentration index (or, alternatively, their morphological type). In the case of the $`r_\mathrm{e}`$-luminosity relation the difference in galaxy properties between low- and high-C31 galaxies is not extreme, and both classes of objects follow a well defined correlation, although with significantly different slopes and scatter (see the left panels in Fig. 8). The effect is more evident in the case of the $`r_e`$-$`\mu _\mathrm{e}`$ correlation, as shown in the right panels of Fig. 8. It is clear from the figure that only for high C31 galaxies, irrespective of their morphological type, $`r_\mathrm{e}`$ and $`\mu _\mathrm{e}`$ are correlated to produce the well known Kormendy relation of early-type galaxies. For low C31 galaxies, instead, no correlation is observed between the two parameters, except for a lack of objects with large $`r_\mathrm{e}`$ and very bright $`\mu _\mathrm{e}`$ (the bottom left corner of the diagram). This result is in qualitative agreement with the finding of Graham & de Block (GdeB01 (2001)) that very late-type spirals (thus galaxies with low C31) have smaller and fainter disks than earlier-type ones. The zone of avoidance at large $`r_\mathrm{e}`$ and very bright $`\mu _\mathrm{e}`$, which is common to galaxies of all C31 values, has been already discussed (for the case of disk galaxies only) by Dalcanton et al. (Dal97 (1997); see their Fig. 4), and has its origin in the existence of an exponential cutoff at the high-mass end of the mass function of galaxies (closely mirrored by the cutoff in the galaxy luminosity function).

## 6 Bulge and disk parameters

The main motivation for using a concentration index in our classification scheme instead of something like a bulge to disk light ratio was the fact that measurements of C31 do not require model-dependent decompositions of a galaxy light distribution, which carry with them systematic uncertainties very difficult to estimate. Still those decompositions are routinely carried out, and are used to analyze in fine detail models of galaxy formation. Courteau et al. (Court96 (1996)) were the first to report the existence of a correlation between disk and bulge scale lengths for two samples of late-type spirals, that they interpreted as evidence for the presence of secular evolution in driving bulge formation. Their result has been recently confirmed also by Khosroshahi et al. (Khos00 (2000)). However these authors, as well as Graham & Prieto (GP99 (1999)), find systematic changes in the ratio of disk to bulge scale length as a function of morphological type, and they interpret this result as evidence against the formation of bulges via secular evolution of the disk of spiral galaxies.

We have analyzed the relation between bulge and disk parameters for the 667 galaxies in our sample that required a bulge+disk decomposition of their surface brightness profile. We find that both the disk and bulge effective radius and disk and bulge effective surface brightness are correlated, although with a large scatter. These correlations are presented in Fig. 9, where different symbols identify galaxies of different morphological type (symbols are the same as in Figs. 3 and 6). There is also some marginal evidence for systematic changes with morphological type in the ratio between the two radii, or the difference between the two surface brightnesses, as shown in Fig. 10. In contrast with the finding of Graham & Prieto (GP99 (1999)), we find that Sa galaxies have a larger ratio between disk and bulge effective radius than later-type spirals, and they also have a much larger difference in effective surface brightness than later types. This difference cannot be considered really significant, as Graham & Prieto (GP99 (1999)) themselves have shown that a change in the function used to fit the surface brightness profile of the bulge (an exponential profile instead of a Sérsic one) leads to a complete reversal of their results, bringing them in agreement with ours. This same caveat applies to the comparison with the results obtained by Courteau et al. (Court96 (1996)): they used a fixed profile shape for fitting the bulge photometric profile, while we have chosen either a de Vaucouleurs $`r^{1/4}`$ or an exponential disk profile according to an objective $`\chi ^2`$ criterion. The large scatter which is superimposed on the correlations shown in Fig. 9 seems to argue against the presence of a single, homogeneous mechanism that could regulate the bulge formation and its interaction with the surrounding disk.

## 7 Discussion and conclusions

A classification scheme for any kind of object, either real or totally abstract ones, is really useful only if it can be used to sort properties of those objects other than those explicitly used for the classification itself. For example, the Harvard stellar spectra classification scheme, created purely on the basis of spectral shape and spectral line patterns, is still in use because it was later recognized that it sorts stars in terms of surface temperature, and ultimately in terms of mass. As far as galaxies are concerned, Hubble’s tuning fork scheme is by far the most commonly use classification tool. The origin of its success, besides the fact that it is relatively easy to use, is mainly due to the fact that it helps sort important physical galaxy properties, like current stellar population and star formation history, or gas and dust content.

However this sorting is only a partial one, and while clear trends exist from many galaxy properties as a function of Hubble’s type, it is also true that these properties show a large scatter within each type, so that there is a large overlap in their distribution between different types (see Roberts & Haynes RH94 (1994) for a recent review). We have shown two examples of this situation in Fig. 1. These large scatters naturally lead to the idea of including some new parameters in the classification scheme, to increase its discriminating power. However, other limitations within Hubble’s scheme, namely its qualitative and quantized nature, and its failure to classify a significant fraction of known galaxies, severely limit the usefulness of any such extension. We have therefore decided to explore the possibility of replacing Hubble’s scheme altogether. This is not a new idea, as complete replacements were already proposed in the past, most noticeably by Morgan(Mor58 (1958), Mor59 (1959)), Whitmore (Whit84 (1984)), and most recently by BJC. Building on Whitmore’s ideas, we have independently obtained a classification scheme that turns out to be qualitatively very similar to the one proposed by BJC, and is illustrated by Fig. 3.

Contrary to BJC, we do not make here any attempt at reproducing the sub-division among Hubble’s types within our scheme, as we are only interested in building a scheme that can give a complete description of galaxy properties based on a small number of quantitative and easy-to-measure parameters. In agreement with BJC we find that normal galaxies (we do not have peculiar, strongly interactive, or extremely low surface brightness galaxies in our sample) can be described by a three parameter model. Using Whitmore terminology, we can call these three parameters galaxy scale (for us the $`H`$-band luminosity), form (for us the concentration index C31), and color (we use the $`BH`$ color, but any optical or optical-near infrared color can be used for this purpose). The need for three parameters is clearly shown by objects that belong to the two different property regimes existing along the strongly non-linear relations between C31 and either color or luminosity. At low luminosity (see the bottom panels of Fig. 3) we need color information to distinguish between a dE and a dIrr galaxy with the same C31, while at high luminosity we need C31 information to distinguish between an S0 and an Sc galaxy with the same color.

The three parameters are relatively easy to derive from photometric observation data, and do not require any modeling of the galaxy light distribution. However the measurement of C31 imposes constraints on the resolution of the imaging data that could be used to classify galaxies, and it is clear that even using our model it would be extremely difficult to rely on ground-based observations to classify high-redshift galaxies.

An alternative representation of the classification cube can be obtained substituting galaxy luminosity with effective surface brightness, taking advantage of the correlation existing between the two quantities (see Figs. 6 and 7). This is perhaps the most unexpected result of our analysis, as it is in contrast with previous results based on observations carried out at shorter wavelengths. We find that a single, well-defined relation exists between luminosity and $`\mu _\mathrm{e}`$ in the $`H`$-band, as opposed to the two separate regimes observed by Binggeli et al. (Bing84 (1984)) in the $`B`$-band, and between $`BH`$ or $`BV`$ color and $`\mu _\mathrm{e}`$, as opposed to the much shallower relation between $`BV`$ and $`\mu _\mathrm{e}`$ in the $`B`$-band presented by BJC. Actually, for $`BV<0.8`$ the BJC relation is completely flat, as demonstrated by the vertical boundary between late- and intermediate-type galaxies in the upper panel of their Fig. 3, and as should be expected on the basis of Freeman’s law (Freeman Free70 (1970)). On the contrary, we find that surface brightness in the near-infrared is not a scale-free parameter, but it depends on the galaxy luminosity, confirming earlier results by Gavazzi et al. (1996c ) and de Jong (deJ96 (1996)).

These differences highlight the importance of selecting the “right” photometric band for mapping the structural properties of galaxies. In this respect near-infrared bands have two key advantages over optical ones: reduced dust extinction effects, and reduced effects of current star formation activity on both the galaxy total luminosity and surface brightness distribution. For these reasons, near-infrared luminosity traces quite accurately galaxy mass (see for example Gavazzi et al. 1996c , and Scodeggio et al. tilt (1998)). Using this luminosity-mass equivalence we can obtain an accurate picture of the role that galaxy mass has in determining the global structure of a galaxy, via the influence it has in determining the galaxy luminous matter mean density and global distribution.

If the observed correlation between mass and surface brightness is already an indication in favor of a scale-dependent galaxy collapse mechanism, the strong non-linearity of the relation between mass and concentration index C31 should pose very strong constraints on such a scenario. Strong light concentrations (high C31 values), that are somewhat correlated to the presence of a dominant spheroidal, bulge-like component (but not in a simple way, as the scatter in the C31-B/T ratio seen in Fig. 2 demonstrates), appear only for luminous red galaxies, at $`L_H>10^{10}L_{}`$. Below that limit galaxies do not have enough mass to drive the concentration of mass and luminous matter towards their center. Above the limit this concentration is possible, but it does not take place automatically for all galaxies. If it happens, however, it must happen quite early during the evolution of a galaxy, or at least involve only old stellar populations without the onset of new star formation activity. In fact high values of C31 are present only among very red galaxies that must have formed a large fraction of their stars during the first few (one to three) billion years of their history.

This important role of galaxy mass in determining the galaxy structure matches a similar role that mass has in determining the star formation history of galaxies, at least for spiral ones (see Gavazzi & Scodeggio GS96 (1996) and also Boissier & Prantzos BP00 (2000)). Unfortunately the link between galaxy mass and near-infrared luminosity is the only well established one along the path from a phenomenological classification scheme to a self consistent physical model, and we don’t have other simple physical properties that can be unequivocally associated with star formation (i.e. color) and matter distribution within a galaxy, although it is obvious that angular momentum must be playing an important role in determining galaxy structure (see, for example, Dalcanton et al. Dal97 (1997); Boissier & Prantzos BP00 (2000); Pierini et al. Dan01 (2002)). It is therefore not possible to build, based on the present data, a simple physical model that could be used, for example, to explore the origin of the scatter that is still present in all relations involving galaxy structural parameters.


text:
# Quantum solution to the Newcomb’s paradox

## 1 Introduction

There is a common belief that the characteristic size of the brain’s integral parts is too big to allow for quantum effects being important . But recent experiments show that separated objects of the size of a golf ball can form quantum entangled states even in a room temperature . Physicists successfully apply quantum mechanics to describe a lot of complex system that may have in principle arbitrary size, including black holes or even the whole Universe. Are there any reasons for quantum modeling of phenomena related to brain activity, consciousness or social behaviour? One can give an answer to this question only after construction and thorough verification of respective models . Below we consider a problem easily susceptible of modeling as a quantum game that should shed some light on the solutions that quantum theory may offer.

In 1960 William Newcomb, a physicist, intrigued the philosopher Robert Nozick with a claim that in an elementary game characterized by the matrix $`M`$

$$M:=\left(\begin{array}{cc}\$1000& \$\mathrm{1\hspace{0.17em}001\hspace{0.17em}000}\\ 0& \$\mathrm{1\hspace{0.17em}000\hspace{0.17em}000}\end{array}\right)$$
(1)
giving the pay-off of the player 1 in all possible situations, the player 1 is not able to chose his strategy without having any measure of occurring a posteriori of any of the four possible events. Rows correspond the player 1’s strategies: feminine $`|\mathrm{𝟢}_1`$ and masculine $`|\mathrm{𝟣}_1`$111The use of the adjectives feminine and masculine to underline the character of the strategies will be explained later, see also the Gardner book and columns to opponent’s strategies $`|\mathrm{𝟢}_2,|\mathrm{𝟣}_2`$. It so happens even despite the fact that the feminine strategy dominates the masculine one (that is the pay-off is greater regardless of the opponents strategy). The choice of the masculine strategy $`|\mathrm{𝟣}_1`$ is more profitable when the event corresponding to the off-diagonal elements of $`M`$ do not occur and the rest have almost equal probabilities. This might happen if the opponent is able to foresee the player 1 moves. Due to this paradoxical property the above game with indefinite (hidden) set of occurrences became for philosophers, economists and theologians a graceful theme of speculations about free will and its consequences . The disputes, often referred to as newcombmania , deserve a thorough analyzis from the quantum game theory point of view -. The development of the probability theory provide us with many intriguing examples where ambiguous specification of the appropriate probability measures resulted in contradiction (Bertrand  and Banach-Tarski  paradoxes are the most famous ones). One can still find people who regardless of this facts continue philosophical disputes while ignoring the necessity of precise definition of the probabilistic measures in their models. We would like to show that quantum theory may be of help in settling the ambiguities.

## 2 Quantum description of the game

Quantum game theory exploits the formalism of quantum mechanics in order to offer the players new classes of strategies. Interesting generalization of well known classical games have been put forward . There are arguments that quantum strategies may offer extraordinary tools for biologists -. Economics being the theater of various games and conflicts should not despise these new ideas . We will describe player’s strategies as vectors (often referred to as states) in Hilbert spaces $`_i`$ where the subscripts $`i=1,2`$ distinguish between the player 1 and 2. It is convenient to define the strategy density operator $`𝒲`$

$$𝒲=\underset{r,s=1}{\overset{2}{}}W_{rs}|𝗋\mathrm{𝟣}_1|𝗌\mathrm{𝟣}_2_1𝗋\mathrm{𝟣}|_2𝗌\mathrm{𝟣}|$$

where $`(W_{rs})`$ is a matrix with nonnegative entries such that $`_{r,s}W_{rs}=1`$ and $`|𝗋_1|𝗌_2_1𝗋|_2𝗌|`$, $`r,s\{\mathrm{𝟢},\mathrm{𝟣}\}`$ are projective operators on the states of the game, $`|𝗋_1|𝗌_2_1_2`$. For our aims it will be sufficient to use two dimensional Hilbert spaces for the players’ strategies. The states of the classical setting (mixed strategies) are represented by a diagonal matrix $`(W_{rs})`$. Non-diagonal elements of $`(W_{rs})`$ describe situations (strategies) that are out of the reach for classical players. Following the classical terminology we will call the pay-off observable $``$ a Hermitian operator corresponding to the matrix $`(\text{1})`$:

$$:=\underset{r,s=1}{\overset{2}{}}M_{rs}|𝗋\mathrm{𝟣}_1|𝗌\mathrm{𝟣}_2_1𝗋\mathrm{𝟣}|_2𝗌\mathrm{𝟣}|.$$

Therefore, according to the classical interpretation of the game, the player 1’s expected pay-off is equal to the sum of diagonal elements (trace) of the product of $`M`$ and the transpose of $`W`$:

$$E():=\text{Tr}𝒲=\underset{rs}{}M_{rs}W_{rs}=\text{Tr}MW^T.$$

## 3 Newcomb’s paradox

M. Gardner proposed the following fabulous description of a game with pay-off given by the matrix $`(\text{1})`$ . An alien Omega (or Alf?) being a omniscent representative of alien civilization (player 2) offers a human (player 1) a choice between two boxes. The player 1 can take the content of both boxes or only the content of the second one. The first one is transparent and contains $1000. Omega declares to have put into the second box that is not transparent $1000000 (strategy $`|\mathrm{𝟣}_2`$) but only if he foresaw that the player 1 decided to take only the content of that box ($`|\mathrm{𝟣}_1`$). A male player 1 thinks: If Omega knows what I am going to do then I have the choice between $1000 and $1000000. Therefore I take the $1000000 (strategy $`|\mathrm{𝟣}_1`$). A female player 1 thinks: Its obvious that I want to take the only the content of the second box therefore Omega foresaw it and put the $1000000 into the box. So the one million dollar is in the second box. Why should I not take more – I take the content of both boxes (strategy $`|\mathrm{𝟢}_1`$). The question is whose strategy, male’s or female’s, is better? One cannot give unambiguous answer to this question without precise definition of the measures of the events relevant for the pay-off.

## 4 Human’s and Omega’s strategies

Omega as representative of an advanced alien civilization is certainly aware of quantum properties of the Universe that are still obscure or mysterious to humans. The boxes containing pay-offs are probably coupled. One can suspect this because the human cannot take content of the transparent box only ($1000). The female player is sceptical about the possibility of realization of the Omega’s scenario for the game. She thinks that the choice of the male strategy results in Omega putting the one million dollar in the second box, and after this being done no one can prevent from her taking the content of the both boxes in question (ie $1001000). But Meyer proposed a quantum tactics  that, if adopted by Omega, allows Omega to accomplish his scenario. Let us note that Omega may not be able to foresee the future . For it aims it is sufficient that it is able to discern human intentions regardless of their will or feelings on the matter. The obstacles to this implied by the no-cloning theorem can be overcome by means of teleportation : Omega has must be able to intercept and then return human’s strategies. The presented below manipulations leading to thwarting humans are feasible with contemporary technologies. The course of the game may look as follows. At the starting-point, the density operator $`𝒲`$ acting on $`_1_2`$ describes the human’s intended strategy and the Omega’s strategy based on its prediction of human’s intentions. The actual game must be carried on according to quantum rules that is players are allowed to change the state of the game by unitary action on $`𝒲`$ . The human player can only act on her/his $`q`$-bit Hilbert space $`_1`$. Omega’s tactics must not depend on the actual move performed by the human player (it may not be aware of the human strategy): its moves are performed by automatic device that couples the boxes. The Meyer’s recipe leads to:

1. Just before the human’s move, Omega set the automatic devise according to its knowledge of human’s intention. The device executes the tactics $``$, where $``$ is the identity transform (Omega cannot change its decision) and $``$ is the well known Hadamard transform frequently used in quantum algorithms: $`F:=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)`$.
2. The human player with the probability $`w`$ uses the female tactics $`𝒩`$ , where $`𝒩`$ is the negation operator222$`𝒩|\mathrm{𝟢}=|\mathrm{𝟣}`$, $`𝒩|\mathrm{𝟣}=|\mathrm{𝟢}`$ and with the probability $`1w`$ the male tactics $``$.
3. At the final step the boxes are being opened and the built-in coupling mechanism performs once more the transform $``$ and the game is settled.

## 5 The course of the game and its result

Let us analyze the evolution of the density operator $`𝒲`$. The players’ tactics, by definition, could have resulted in changes in the (sub-)space $`_1`$ only therefore it suffices to analyze the human’s strategies. In a general case the human can use a mixed strategy: the female one with the probability $`v`$ and the male one with the probability $`1v`$. Let us begin with the extreme values of $`v`$ (pure strategies). If the human decided to use the female strategy ($`v=1`$) or the male one ($`v=0`$) then the matrices $`𝒲_i`$, $`i=0,1`$ corresponding to the density operators

$$𝒲_0=\underset{r,s=1}{\overset{2}{}}W_{0}^{}{}_{rs}{}^{}|𝗋\mathrm{𝟣}_1|\mathrm{𝟢}_2_1𝗌\mathrm{𝟣}|_2\mathrm{𝟢}|$$

and

$$𝒲_1=\underset{r,s=1}{\overset{2}{}}W_{1}^{}{}_{rs}{}^{}|𝗋\mathrm{𝟣}_1|\mathrm{𝟣}_2_1𝗌\mathrm{𝟣}|_2\mathrm{𝟣}|$$

are calculated as follows:

$$\begin{array}{cc}& \left(\begin{array}{cc}v& 0\\ 0& 1v\end{array}\right)\frac{1}{2}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)\left(\begin{array}{cc}v& 0\\ 0& 1v\end{array}\right)\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)=\frac{1}{2}\left(\begin{array}{cc}1& 2v1\\ 2v1& 1\end{array}\right)\hfill \\ & \frac{w}{2}\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)\left(\begin{array}{cc}1& 2v1\\ 2v1& 1\end{array}\right)\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)+\frac{1w}{2}\left(\begin{array}{cc}1& 2v1\\ 2v1& 1\end{array}\right)=\frac{1}{2}\left(\begin{array}{cc}1& 2v1\\ 2v1& 1\end{array}\right)\hfill \\ & \frac{1}{4}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)\left(\begin{array}{cc}1& 2v1\\ 2v1& 1\end{array}\right)\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)=\left(\begin{array}{cc}v& 0\\ 0& 1v\end{array}\right).\hfill \end{array}$$

It is obvious that independently of the used tactics, human’s strategy takes the starting form. For the mixed strategy the course of the game is described by the density operator

$$𝒲=v𝒲_0+(1v)𝒲_1$$

which also has the same diagonal form at the beginning and at the end of the game:

$$𝒲=v|\mathrm{𝟢}_1|\mathrm{𝟢}_2_1\mathrm{𝟢}|_2\mathrm{𝟢}|+(1v)|\mathrm{𝟣}_1|\mathrm{𝟣}_2_1\mathrm{𝟣}|_2\mathrm{𝟣}|.$$

Therefore the change of mind resulting from the female strategy cannot lead to any additional profits. If the human using the female tactics (that is changes his/her mind) begins the game with the female strategy then at the end the untransparent box will be empty and he/she will not get the content of the transparent box: the pay-off will be minimal (0). If the human acts just the opposite the transparent box must not be opened but nevertheless the pay-off will be maximal ($100000). Only if the human begins with the female strategy and then applies the male tactics the content of the transparent box is accessible. If restricted to the classical game theory Omega would have to prevent humans from changing their minds. In the quantum domain the pay-off $`M_{21}`$ (female strategy and tactics) is possible (the phrase la donna mobile gets a quantum context): humans regain their free will but they have to remember that Omega has (quantum) means to prevent humans from profiting from altering their decisions. In that way quantum approach allows to remove the paradox from the rationally defined dilemma. One can also consider games with more alternatives for the human player. The respective larger pay-off matrices would offer even more sophisticated versions of the Newcomb’s observation. But even then there is a quantum protocol that guarantees that Omega keeps its promises (threats) .

## 6 Market interpretation of the game

It is obvious that the above scenario cannot be realized if the actual conditions would differ from Omega’s promises. For example, Omega may not be able to predict humans intentions or its understanding of the rules of the game differs from that implied by their expression in human language (cultural differences). There may be much dispute over the question what Omega really has in mind? We would like to consider one of the variant that may be interesting in the context of quantum market games . This may result from pondering over the meaning of the term Omega adopts the same strategy.

Players in a quantum market game sometimes buy and sometimes sell. A demand representation of the player’s strategy is a Fourier transform of his strategy used while supplying the goods . In a simplified model where player’s strategies span a finite dimensional Hilbert space we should apply discrete Fourier transform which transforms the demand representation of the strategy, being $`m`$-tuple of complex numbers $`d|\psi `$ to the supply representation given by

$$s|\psi =\frac{1}{\sqrt{m}}\underset{d=0}{\overset{m1}{}}\text{e}^{\frac{2\pi \text{i}}{m}sd}d|\psi .$$

If $`m=2`$ then the discrete Fourier transform reduces to the Hadamard transform $``$ which we have already met. In our case the Hadamrd transform switches maximally localized strategies with the the maximally indefinite strategies and vice versa, eg. $`d|\psi =[d=0]\stackrel{}{}s|\psi =\frac{1}{\sqrt{2}}`$ (the Iverson notation  is used: $`[expression]`$ denotes the logical value (1 or 0) of the sentence $`expression`$).

We introduce the nonhomogeneous complex coordinate $`z\overline{}`$ to parameterize player’s strategies

$$|\psi _z:=|\mathrm{𝟢}+z|\mathrm{𝟣}.$$

If the "buying human" decides to use the strategy $`|\psi _z_1`$ and the other side of the bargain (Omega) want to play in the same way and therefore uses the supply representation of human’s strategy setting its $`q`$-bit to $`(|\mathrm{𝟢}_2+z|\mathrm{𝟣}_2)=|\mathrm{𝟢}_2+\frac{1z}{1+z}|\mathrm{𝟣}_2`$ then the quantum state of the game takes the form333Cases when the player may in fact play also ”against himself” often happen in market description: demand or supply result from self-consistent strategy of all players. For example in a stock exchange one big bid or transaction can influence the whole market

$$𝒲_z=\frac{(1+z)(1+\overline{z})}{2(1+z\overline{z})^2}(|\mathrm{𝟢}_1+z|\mathrm{𝟣}_1\left)(_1\mathrm{𝟢}|+\overline{z}_1\mathrm{𝟣}|)(|\mathrm{𝟢}_2+\frac{1z}{1+z}|\mathrm{𝟣}_2\right)(_2\mathrm{𝟢}|+\frac{1\overline{z}}{1+\overline{z}}_2\mathrm{𝟣}|).$$

Therefore, as in the previous discussion, the female strategy gives not higher a pay-off. The expectation value of the human’s pay-off, $`\text{Tr}𝒲_z`$, is maximal for a superposition of male and female strategies with phase shifted by $`\pi `$ (i.e. for $`z=1`$). In this case the human is better off than in the previous case but she or he must be cautious because the phase shift by $`\pi `$ ($`z=1`$) does not change the respective probabilities but result in the lowest expectation value of the pay-off ($500). Classical human’s strategies correspond to $`z=0`$ (female) and $`z=\pm \mathrm{}`$ (male). The expectation values of the human pay-off with respect to the adopted strategy are presented in Figure 1.

Enthusiasts for newcombmania will certainly find a lot of new quantum solutions to the Newcomb game.


text:
# Nonlinear anomalous diffusion equation and fractal dimension: Exact generalized gaussian solution

## I Introduction

Anomalous diffusion has been the subject of numerous investigations in the last years. In particular, anomalous diffusion has played a fundamental role in the analysis of a wide class of systems such as diffusion in plasmas berryman , turbulent flows shlesinger , transport of fluid in porous media spohn , and diffusion on fractals stephenson . In general, anomalous diffusion may be classified employing the mean square displacement, $`(\mathrm{\Delta }x)^2`$. When $`(\mathrm{\Delta }x)^2t^\eta `$, the value $`\eta >1`$ characterizes a superdiffusive process, $`\eta <1`$ a subdiffusive one and $`\eta =1`$ a ‘normal’ diffusion.

We consider here a quite large class of anomalous diffusion, namely those associated with the $`N`$-dimensional nonlinear equationfootnote

$$\frac{\rho (𝐫,t)}{t}=\{K[\rho (𝐫,t)]^\nu \}[\mu 𝐅(𝐫,t)\rho (𝐫,t)]\alpha (t)\rho (𝐫,t),$$
(1)
where $`K=Dr^\theta `$, $`\nu `$, $`\theta `$, and $`D`$ are real parameters, $`\mu `$ is the drift coefficient, $`𝐅(𝐫,t)`$ is the external force, and $`\alpha (t)`$ is a time-dependent source. From the Gauss theorem one verifies that this general equation preserves the norm $`\rho (𝐫,t)d^Nr`$ if $`\alpha (t)=0`$ and $`K[\rho (𝐫,t)]^\nu \mu 𝐅(𝐫,t)\rho (𝐫,t)`$ goes sufficiently rapid to zero when $`r=|𝐱|\mathrm{}`$. The difference between the Fokker-Planck-like equation (1) and the usual one is that the diffusion term depends on a power of $`\rho (𝐫,t)`$ and the diffusion coefficient is spatially dependent. Notice that this equation for $`\theta =0`$ and $`\alpha (t)=0`$ has the form:

$$\frac{\rho (𝐫,t)}{t}=D^2[\rho (𝐫,t)]^\nu [\mu 𝐅(𝐫,t)\rho (𝐫,t)],$$
(2)
and it is usually called porous media equation. The above equation admits important physical applications such as percolation of gases through porous media muskat , thin liquid films spreading under gravity buckmaster , some self-organizing phenomena carlson , and surface grow spohn , among others (see Ref. tsallis1  and references therein). Furthermore, this equation has been related to nonextensive Tsallis statistics. This connection was first pointed by Plastino and Plastino plastino  and revisited by others authors tsallis1 ; stariolo ; compte ; compte2 ; martinez ; lisa2 ; lisa1 .

Another important kind of anomalous diffusion emerges from Eq. (1) when $`\nu =1`$, $`\alpha (t)=0`$ and $`𝐅(𝐫,t)=0`$ . In this case, the equation

$$\frac{\rho (𝐫,t)}{t}=[K\rho (𝐫,t)]$$
(3)
was related to turbulent diffusion in the atmosphere richardson  (see also Ref. procaccia2 ). Here $`Kr^{4/3}`$. In general, in a $`d`$-dimensional space, $`K`$ is proportional to $`\frac{1}{r^{d1}}\frac{}{r}r^{d1\theta }\frac{}{r}+\frac{A}{r^{2+\theta }}`$ ($`A`$ is an operator, which depends on the angular variables). Thus,

$$\stackrel{~}{\mathrm{\Delta }}=\frac{1}{r^{d1}}\frac{}{r}r^{d1\theta }\frac{}{r}$$
(4)
is the radial part of the operator $`K`$ to be considered in the study of $`d`$-dimensional spherical symmetrical solutions of Eq. (3). If $`d`$ is interpreted as a positive real number, it plays a role of a fractal dimension embedding in a $`N`$-dimensional space. Thus, Eq. (3) becomes, in a fractal framework, the O’Shaughnessy-Procaccia equationprocaccia ,

$$\frac{\rho (r,t)}{t}=D\stackrel{~}{\mathrm{\Delta }}\rho (r,t).$$
(5)
We are ready now to unify Eqs. (2) and (5). In fact, if we consider only diffusive terms and the radial part of these equations we have mala

$$\frac{\rho (r,t)}{t}=D\stackrel{~}{\mathrm{\Delta }}[\rho (r,t)]^\nu .$$
(6)
Of course, if $`\theta =0`$ and $`\nu =1`$, the usual diffusion equation is recovered, i.e.,

$$\frac{\rho (𝐫,t)}{t}=D^2\rho (𝐫,t).$$
(7)
In the previous scenario, Eq. (1) can be generalized in order to be also employed in a fractal context, when spherical symmetric solutions are considered and $`K`$ is replaced by $`\stackrel{~}{\mathrm{\Delta }}`$. This paper is dedicated to investigate a class of these solutions in the presence of a source term and with a drift term, as well as the stationary solutions when a general drift term is present.

The present work is organized as follows. In the next section we review an exact solution (generalized gaussian) for Eq. (6) and its limit cases. This generalized gaussian essentially gives the basis for the remaining solutions investigated here. In Sec. III, a formal solution for the nonlinear Fokker-Planck equation Eq. (6) with a source term is obtained. Sec. IV deals with the application of our formal solution for a time power law source. In Sec. V, we present a generalized gaussian solution for the nonlinear Fokker-Planck with a linear spherical force term. Furthermore, in Sec. VI, we consider the stationary solutions for nonlinear Fokker-Planck equation with a general drift term. A summary of our results is presented in the last section.

## II Unified power-law and stretched exponential solution

Motivated by the Gaussian solution of the usual diffusion, Eq. (7), it is natural to employ a generalization of such solution when we are investigating Eq. (6). In this direction, we follow Ref. mala  employing the ansatz:

$$\rho (r,t)=\left[1(1q)\beta (t)r^\lambda \right]^{1/(1q)}/Z(t)$$
(8)
if $`1(1q)\beta (t)r^\lambda 0`$, and $`\rho (r,t)=0`$ if $`1(1q)\beta (t)r^\lambda <0`$ (cutoff condition). By replacing this generalized gaussian function in Eq.(6), we verify that $`\lambda =2+\theta `$, $`q=2\nu `$,

$$\beta (t)=\beta _0[1+At]^{\lambda /[\lambda +d(1q)]},$$
(9)
and

$$Z(t)=Z_0[1+At]^{d/[\lambda +d(1q)]},$$
(10)
where $`A=D\lambda (2q)[\lambda +d(1q)]\beta _0Z_0^{q1}`$, $`\beta _0=\beta (0)`$ and $`Z_0=Z(0)`$. It is also important to note that the product $`\beta Z^{\lambda /d}`$ is time independent.

From the solution (8) the mean value for a generic function $`f(r)`$ can be calculated. In particular, the mean value for $`r^\sigma `$ is given by

$$r^\sigma =\frac{_0^{\mathrm{}}r^\sigma \rho (r,t)r^{d1}𝑑r}{_0^{\mathrm{}}\rho (r,t)r^{d1}𝑑r}=C_\sigma \beta ^{\sigma /\lambda },$$
(11)
with

$$C_\sigma =\frac{\mathrm{\Gamma }(\frac{\sigma +d}{\lambda })}{\mathrm{\Gamma }(\frac{d}{\lambda })}\{\begin{array}{cc}\frac{\mathrm{\Gamma }(\frac{1}{1q}+\frac{d}{\lambda }+1)}{(1q)^{\sigma /\lambda }\mathrm{\Gamma }(\frac{\sigma +d}{\lambda }+\frac{1}{1q}+1)}\hfill & \text{ }q<1\hfill \\ \frac{\mathrm{\Gamma }(\frac{1}{q1}\frac{d+\sigma }{\lambda })}{(q1)^{\sigma /\lambda }\mathrm{\Gamma }(\frac{1}{q1}\frac{d}{\lambda })}\hfill & q>1.\hfill \end{array}$$
(12)
Observe that the existence of $`r^\sigma `$ imposes restrictions over the parameters: $`\lambda +d(1q)>\sigma (q1)`$ and $`\lambda >0`$. In addition, the restriction $`q<2`$ ($`\nu >0`$) is necessary for $`\rho (r,t)`$ to be real. In the following, we assume that the parameters obey the above restrictions.

An important particular case of $`r^\sigma `$ is the mean square displacement, $`r^2`$, since such mean value is commonly employed to classify the diffusion process. Thus, by taking $`\sigma =2`$ in Eqs. (11) and (12) we have

$$r^2=C_2\beta _0^{2/\lambda }[1+At]^{2/[\lambda +d(1q)]}.$$
(13)
Consequently, $`r^2t^{2/[2+\theta +d(\nu 1)]}`$ for a sufficient large $`t`$. Now, it is easy to verify that $`\theta =d(1\nu )`$ describes a ‘normal’ diffusion, even when $`\theta 0`$ and $`\nu 1`$. In fact, there is a competition between $`\theta `$ and $`\nu `$ in such way that the anomalous diffusive regime induced by $`\theta 0`$ is compensated by a convenient one with $`\nu 1`$. Furthermore, this competition leads to a subdiffusive process when $`\theta >d(1\nu )`$, and to a superdiffusive one when $`\theta <d(1\nu )`$.

Note, also, that an alternative point of view to verify how much a solution deviates from the usual one is to analyze the decay of $`\rho (0,t)`$. In other words, the way the solution decreases in this central point can inform about its deviation from the normal diffusion. In our case, Eq. (8) gives $`\rho (0,t)t^{d/[2+\theta +d(\nu 1)]}`$ in contrast with $`\rho (0,t)t^{d/2}`$, obtained from the usual diffusion.

We remark that important limiting cases can be obtained from Eq. (8). If $`\theta =0`$, the solution for the porous media equation without the drift term is recovered tsallis1 ; plastino ,

$$\rho (r,t)=\left[1(1q)\beta _1(t)r^2\right]^{\frac{1}{1q}}/Z_1(t),$$
(14)
with $`\beta _1(t)t^{2/[2+d(1q)]}`$ and $`Z_1(t)t^{d/[2+d(1q)]}`$. On the other hand, for $`\nu =1`$, the point-source solution for the O’Shaugnessy-Procaccia equationprocaccia  emerges,

$$\rho (r,t)=e^{\beta _2(t)r^{\theta +2}}/Z_2(t),$$
(15)
where $`\beta _2t^1`$ and $`Z_2(t)t^{d/(2+\theta )}`$. Moreover, Eq. (8) recovers the time-translated Gaussian solution for the Eq. (7),

$$\rho (r,t)=\frac{\rho _0}{[4\pi D(t_0+t)]^{d/2}}e^{\frac{r^2}{4D(t_0+t)}}.$$
(16)
In fact, Eqs. (9) and (10) can be respectively written as

$$\beta (t)=\stackrel{~}{A}(t_0+t)^{\lambda /[\lambda +d(1q)]}$$
(17)
and

$$Z(t)=\stackrel{~}{B}(t_0+t)^{d/[\lambda +d(1q)]},$$
(18)
with $`\stackrel{~}{A}=\beta _0A^{\lambda /[\lambda +d(1q)]}`$, $`\stackrel{~}{B}=Z_0A^{d/[\lambda +d(1q)]}`$, and $`t_0=1/A`$. Eq. (16) is obtained from Eq. (8) with $`\beta (t)`$ and $`Z(t)`$ given above in the limit $`\nu =1`$ and $`\lambda =2`$, and by using the normalization $`\mathrm{\Omega }_d\rho (r,t)r^{d1}𝑑r=\rho _0`$, where $`\mathrm{\Omega }_d=2\pi ^{d/2}/\mathrm{\Gamma }(d/2)`$ is the solid angle in a $`d`$-dimensional space.

To conclude this review, we notice that our point-like solution can be connected with the Tsallis entropymala .

## III Nonlinear Fokker-Planck equation with time-dependent source term

Nonlinear diffusion with absorption arises in many areas of science such as engineering, biophysics, solid state physics and reactor, among othersscott ; murray ; aris ; meyrs . This section is dedicated to the analysis of the generalized Fokker-Planck equation presented in the introduction with a source term. This kind of equation is expected to provide a basis for phenomenological applications, where a large class of nonlinear anomalous diffusion with source plays an important role.

Let us consider the nonlinear Fokker-Planck equation (1). The source term in this equation can be removed by an appropriate change in the solution:

$$\rho (𝐫,t)=\mathrm{exp}\left(_0^t\alpha (t^{})𝑑t^{}\right)\widehat{\rho }(𝐫,t).$$
(19)
In this way, $`\widehat{\rho }(𝐫,t)`$ obeys the equation

$$\frac{\widehat{\rho }(𝐫,t)}{t}=\{\widehat{K}[\widehat{\rho }(𝐫,t)]^\nu \}[\mu 𝐅(𝐫,t)\widehat{\rho }(𝐫,t)],$$
(20)
where $`\widehat{K}(t)=K\mathrm{exp}\left[(1\nu )_0^t\alpha (t^{})𝑑t^{}\right]`$. Thus, Eq. (20) has the same form of the corresponding one without the source term, but with a time-dependent diffusion coefficient $`\widehat{K}`$. Note that this time-dependence is induced by the nonlinear term $`\rho ^\nu `$, disappearing when $`\nu =1`$.

As emphasized in the introduction, below Eq. (7), the replacement of the operator $`(\widehat{K})`$ by $`\stackrel{~}{\mathrm{\Delta }}`$ unable us also to analyze situations with fractal dimension in the sense of O’Shaughnessy-Procaccia. Thus, to investigate a radial nonlinear anomalous diffusion equation with a source term and non-integer dimension $`d`$, we consider

$$\frac{\rho (r,t)}{t}=D\stackrel{~}{\mathrm{\Delta }}[\rho (r,t)]^\nu \alpha (t)\rho (r,t)$$
(21)
instead of Eq. (1). In the interests of brevity we omitted the drift term. The case with an external force is discussed in Sec. V. To eliminate the source term, we proceed analogously as in the previous paragraph. Therefore, we obtain

$$\frac{\widehat{\rho }(r,t)}{t}=\widehat{D}\stackrel{~}{\mathrm{\Delta }}[\widehat{\rho }(r,t)]^\nu ,$$
(22)
with $`\widehat{D}(t)=D\mathrm{exp}\left[(1\nu )_0^t\alpha (t^{})𝑑t^{}\right]`$. Now, we redefine the time variable introducing an effective time as

$$\tau (t)=_0^t\widehat{D}(t^{})𝑑t^{},$$
(23)
so Eq. (22) can be written as

$$\frac{\widehat{\rho }(r,\tau )}{\tau }=\stackrel{~}{\mathrm{\Delta }}[\widehat{\rho }(r,\tau )]^\nu .$$
(24)
The solution of this equation was presented and discussed in Sec. II. Thus, it is given by Eqs. (8), (9) and (10) with $`D=1`$ and $`t`$ replaced by $`\tau `$. Another possible generalized Fokker-Planck is

$$\frac{[\rho (𝐫,t)]^\eta }{t}=\{K[\rho (𝐫,t)]^\nu \}\{\mu 𝐅(𝐫,t)[\rho (𝐫,t)]^\eta \}\alpha (t)[\rho (𝐫,t)]^\eta .$$
(25)
This equation, with $`K`$ and $`\alpha `$ being constants, was discussed in Ref. tsallis2 . Note also that Eq. (25) reduces to Eq. (1) when we replace $`[\rho (r,t)]^\eta `$ by $`\stackrel{~}{\rho }(r,t)`$. In the same direction, this kind of generalization can be incorporated in Eq. (21). In this paper, we will not present details about these possibilities because their solutions can be directly obtained from Eqs. (1) or (21) by using the identification $`\stackrel{~}{\rho }=\rho ^\eta `$.

## IV Some kinds of sources

In order to have an enhanced understanding of the relevance of the source terms in the anomalous diffusion equation, we apply the previous general procedure to study some specific cases. In this context, the simplest situation to be investigated is when $`\alpha (t)=\alpha _0`$, where $`\alpha _0`$ is a constant. We also investigate a more generic source term, $`\alpha (t)=\alpha _nt^n`$.

### IV.1 Constant source

Now, taking into account Eq.(21), we analyze the case $`\alpha (t)=\alpha _0`$. Thus, the effective time, Eq. (23), becomes

$$\tau (t)=D\frac{\mathrm{exp}[(1\nu )\alpha _0t]1}{(1\nu )\alpha _0}.$$
(26)
This effective time exhibits three different behaviors. For $`(1\nu )\alpha _0<0`$, the effective time goes exponentially to $`\tau _{\mathrm{}}=D/[(\nu 1)\alpha _0]`$, i.e., a stabilization occurs in the diffusive process. For $`(1\nu )\alpha _0>0`$, there is an exponential increase leading to $`\tau (t)D\mathrm{exp}[(1\nu )\alpha _0t]/[(1\nu )\alpha _0]`$, when $`t`$ is large. It is interesting to remind that this kind of behavior arises due to the nonlinearity of the diffusion equation. In this direction, we have $`\tau =Dt`$ when $`\nu =1`$, which corresponds to $`q=1`$ as long as $`q=2\nu `$.

By using Eqs. (9) and (10) with $`t`$ replaced by $`\tau `$, we obtain

$$\beta (\tau (t))=\beta _0\left(1+A\frac{\mathrm{exp}[(1\nu )\alpha _0t]1}{(1\nu )\alpha _0}\right)^{\lambda /[\lambda +d(1q)]}$$
(27)
and

$$Z(\tau (t))=Z_0\left(1+A\frac{\mathrm{exp}[(1\nu )\alpha _0t]1}{(1\nu )\alpha _0}\right)^{d/[\lambda +d(1q)]}.$$
(28)
Thus, from Eqs.(8) and (19) we verify that

$$\rho (r,t)=e^{\alpha _0t}\left[1(1q)\beta (\tau (t))r^\lambda \right]^{1/(1q)}/Z(\tau (t))$$
(29)
if $`1(1q)\beta (\tau (t))r^\lambda 0`$, and $`\rho (r,t)=0`$ if $`1(1q)\beta (\tau (t))r^\lambda <0`$.

In this solution, when $`(1\nu )\alpha _0<0`$, we have, for large time,

$$\rho (r,t)e^{\alpha _0t}\left[1(1q)\beta _{\mathrm{}}r^\lambda \right]^{1/(1q)}/Z_{\mathrm{}},$$
(30)
where

$$\beta _{\mathrm{}}=\beta _0\left\{1+A/[(\nu 1)\alpha _0]\right\}^{\lambda /[\lambda +d(1q)]}$$
(31)
and

$$Z_{\mathrm{}}=Z_0\left\{1+A/[(\nu 1)\alpha _0]\right\}^{d/[\lambda +d(1q)]}.$$
(32)
This means that $`\rho (r,t)`$ decays (grows) exponentially in the time when $`\alpha _0>0(\alpha _0<0)`$, and is spatially limited for $`q<1`$ as long as $`\rho (r,t)=0`$ if $`r^\lambda 1/[(1q)\beta _{\mathrm{}}]`$. In contrast, the spatial asymptotic power law behavior,

$$\rho (r,t)r^{\lambda /(1q)},$$
(33)
arises for sufficiently large $`t`$ and $`(1\nu )\alpha _0>0`$. On the other hand, when $`\nu =1`$, we have

$$\rho (r,t)t^{d/\lambda }\mathrm{exp}[\alpha _0tr^\lambda /(D\lambda ^2t)]$$
(34)
for large $`t`$.

Let us focus attention now on the mean square displacement. By using the effective time in Eq. (13), we obtain

$$r^2=C_2\beta _0^{2/\lambda }\left[1+A\frac{e^{(1\nu )\alpha _0t}1}{(1\nu )\alpha _0}\right]^{2/[\lambda +d(1q)]}.$$
(35)
The long time behavior for $`r^2`$ leads to an exponential growth when $`(1\nu )\alpha _0>0`$. For $`(1\nu )\alpha _0<0`$, this expression converges to an asymptotic value. For $`\nu =1`$, we have $`r^2t^{2/\lambda }`$, i.e., when $`\nu =1`$, the mean square displacement is not affected by the constant source term. These results show that the presence of a source term plays a relevant role on the diffusive process when $`\nu 1`$.

### IV.2 Source with a general power term

Now we analyze the more general case corresponding to the source $`\alpha _n(t)=\alpha _nt^n`$, where $`\alpha _n`$ and $`n`$ are constants. The related effective time, $`\tau _n(t)`$, is given by

$$\tau _n(t)=_{t_0}^tD\mathrm{exp}\left[\frac{(1\nu )\alpha _n}{n+1}\left(t^{n+1}t_0^{n+1}\right)\right]𝑑t^{},$$
(36)
where $`t_0`$ is a cutoff only used to avoid divergence when $`n1`$ and $`(1\nu )\alpha _n/(n+1)>0`$. When $`n>1`$, we get

$$\tau _n(t)=D\frac{1}{n+1}\left(\frac{\alpha _n(1\nu )}{1+n}\right)^{1/(1+n)}\left[\mathrm{\Gamma }\left(\frac{1}{1+n}\right)\mathrm{\Gamma }(\frac{1}{1+n},\frac{\alpha _n(\nu 1)t^{1+n}}{1+n})\right].$$
(37)
Here $`\mathrm{\Gamma }(c)`$ and $`\mathrm{\Gamma }(c,x)`$ are the gamma and the incomplete gamma functions, respectively. By using the asymptotic behavior for $`\mathrm{\Gamma }(c,x)`$, we obtain $`\tau _nt^n\mathrm{exp}[\alpha _n(\nu 1)t^{1+n}/(1+n)]`$ for $`n>1`$ and large $`t`$. For $`n=1`$, we have $`\tau _1(t)t^{\alpha _1(1\nu )+1}`$, and for $`n<1`$, we verify that $`\tau _n(t)t`$. Taking into account the asymptotic behavior of the mean square displacement, $`r^2\tau _n^{2/[\lambda +d(1q)]}`$, we identify two different behaviors when $`n>1`$: if $`\alpha _n(\nu 1)>0`$, the mean square displacement reaches a limit value, and it grows exponentially if $`\alpha _n(\nu 1)<0`$. For $`n=1`$, its value goes with $`t^{2[1(1q)\alpha _1]/[\lambda +d(1q)]}`$. Finally, for $`n<1`$, we obtain $`r^2t^{2/[\lambda +d(1q)]}`$.

## V External forces

Now we draw our attention in order to analyze the generalized radial Fokker-Planck equation with a drift term. We take Eq. (6) with a radial external force, i.e.,

$$\frac{\rho (r,t)}{t}=D\stackrel{~}{\mathrm{\Delta }}[\rho (r,t)]^\nu \frac{1}{r^{d1}}\frac{}{r}\left[r^{d1}F(r)\rho (r,t)\right],$$
(38)
where we set $`\mu `$ equal to unity. In this section we are going to consider a specific case of a particle moving in a harmonic spherical potential $`V(r)=kr^2/2`$. Thus, $`r=0`$ is a stable equilibrium position for positive $`k`$, and $`r=0`$ is an unstable equilibrium position for negative $`k`$. Eq. (38), with this potential, corresponds to the Uhlenbeck-Ornstein processUhlenbeck  in the particular case of $`\theta =0`$ and $`q=1`$ ($`\nu =2q`$).

By using, again, the ansatz (8) into Eq. (38), we verify that $`Z(t)`$ and $`\beta (t)`$ obey the equations:

$$\frac{1}{Z}\frac{dZ}{dt}=Dd\nu \lambda \beta Z^{q1}kd$$
(39)
and

$$\frac{1}{\beta }\frac{d\beta }{dt}=D\nu \lambda ^2\beta Z^{q1}+k\lambda .$$
(40)
The solution for these two nonlinear coupled equations is

$$Z(t)=Z_0\left[\frac{}{}1+\frac{1}{k}\left(D\nu \lambda \beta _0Z_0^{q1}k\right)\left(1e^{k[\lambda +d(1q)]t}\right)\right]^{d/[\lambda +d(1q)]}$$
(41)
and

$$\beta (t)=\beta _0\left[1+\frac{1}{k}\left(D\nu \lambda \beta _0Z_0^{q1}k\right)\left(1e^{k[\lambda +d(1q)]t}\right)\right]^{\lambda /[\lambda +d(1q)]},$$
(42)
whose initial conditions are $`Z(0)=Z_0`$ and $`\beta (0)=\beta _0`$.

Since $`\lambda +d(1q)>0`$, for positive $`k`$, we have a stationary solution and $`r^2`$ goes to a constant value. When $`k<0`$, there is no stationary state and the $`r^2`$ increases exponentially for large $`t`$. Our solution with $`\theta =0`$ corresponds to one obtained by Tsallis and Bukman tsallis1  for the one-dimensional case when $`x=0`$. We also remark that O’Shaugnessy-Procaccia equation, Eq. (5), when an external force is incorporated, is a special case of our results if we take $`\nu =1`$.

Before concluding this section, we note that in the presence of the source term, $`\alpha (t)\rho `$, a similar procedure as presented in Sec. III can be employed, leading to Eq. (38) with $`D`$ replaced by $`\widehat{D}(t)`$. Thus, Eqs. (39) and (40) must be solved with $`\widehat{D}(t)`$ instead $`D`$ to obtain the generalized gaussian solution.

## VI Stationary solutions for the nonlinear Fokker-Planck equation

We can extend our study about Fokker-Planck-like equation by investigating stationary solutions. The usual Fokker-Planck equation has a stationary solution that corresponds to the canonical Boltzmann-Gibbs distribution, $`Pe^{\beta V}`$. On the other hand, the corresponding stationary solution of the nonlinear equation (2) islisa1

$$\rho (𝐫)=\rho _0\left[1(1q)\beta V(𝐫)\right]^{1/(1q)}$$
(43)
if $`1(1q)\beta V(𝐫)0`$, and $`\rho (𝐫)=0`$ if $`1(1q)\beta V(𝐫)<0`$. Here $`\beta =\mu \rho _0^{q1}/[(2q)D]`$, $`\rho _0`$ is a positive integration constant, and $`𝐅(𝐫)=V(𝐫)`$. Notice that $`\beta >0`$ when $`V(𝐫)`$ is a confining potential.

When we consider the unified nonlinear equation

$$\frac{\rho (𝐫,t)}{t}=\{K(𝐫)[\rho (𝐫,t)]^\nu \}[\mu (𝐫)𝐅(𝐫)\rho (𝐫,t)],$$
(44)
its stationary solution that generalizes Eq. (43) is

$$\rho (𝐫)=\rho _0\left[1\left(\frac{1q}{2q}\right)\rho _0^{q1}V_{eff}(𝐫)\right]^{1/(1q)}$$
(45)
if $`1\left(1q\right)\rho _0^{q1}V_{eff}/(2q)0`$, and $`\rho (r)=0`$ if $`1\left(1q\right)\rho _0^{q1}V_{eff}/(2q)<0`$, where we introduced the effective potencial $`V_{eff}(𝐫)=\frac{\mu (𝐫)}{K(𝐫)}V(𝐫)𝑑𝐫`$. Note that $`[\mu (𝐫)/K(𝐫)]V(𝐫)`$ must be written as a gradient so it is possible to obtain an effective potential in the $`N`$-dimensional case. In particular, for a $`3`$-dimensional example, this implies that $`[\mu (𝐫)/K(𝐫)]\times V(𝐫)=0`$. This condition can be easily accomplished if we take $`\mu =\mu (r)`$, $`K=K(r)`$, and $`V=V(r)`$. Thus, Eq.(45) can be written as

$$\rho (r)=\rho _0[1\beta ^{}(1q)V_{eff}(r)]^{1/(1q)},$$
(46)
with $`\beta ^{}=\rho _0^{q1}/(2q)`$ and $`V_{eff}(r)=\frac{\mu (r)}{K(r)}\frac{dV(r)}{dr}𝑑r`$. Following the development presented in the previous sections, we stress that the stationary solution (46) remains true for an arbitrary non-integer dimension. For the one-dimensional case, the solution (45) was obtained in Ref. lisa1 .

## VII Summary

In this paper we have obtained exact spherically symmetric solutions for a nonlinear Fokker-Planck equation with radial-dependent diffusion coefficient, time-dependent source, and spherical external force. We remark that the time-dependent solutions are generalizations of the gaussian ones.

Without both external force and source term, there is a competition between the fractal aspect and the nonlinear behavior (dictated by the parameter $`\nu `$), leading to superdiffusive, ‘normal’ and subdiffusive processes. On the other hand, in the presence of a source term, the mean square displacement is only affected in the nonlinear case, i.e., the source term modifies the superdiffusive, subdiffusive or ‘normal’ characters of the process when $`\nu 1`$. In particular, this behavior is notorious when a power law time-dependent source is considered.

The diffusive process with a external force $`F(r)=kr`$ leads to different behaviors to the time dependence of distribution function. For positive $`k`$, we have a stationary solution and the mean square displacement goes to a constant value. When $`k<0`$ there is no stationary state and the mean square displacement increases asymptotically in an exponential way.

In the case of the stationary solutions for our Fokker-Planck-like equation, we obtained a generalization of the canonical Boltzmann-Gibbs distribution. This generalized solution only exists if we are able to obtain an effective potential that includes the true potential, the drift and the diffusion coefficients.


text:
# Derivation of continuum stochastic equations for discrete growth models

## Abstract

We present a formalism to derive the stochastic differential equations (SDEs) for several solid-on-solid growth models. Our formalism begins with a mapping of the microscopic dynamics of growth models onto the particle systems with reactions and diffusion. We then write the master equations for these corresponding particle systems and find the SDEs for the particle densities. Finally, by connecting the particle densities with the growth heights, we derive the SDEs for the height variables. Applying this formalism to discrete growth models, we find the Edwards-Wilkinson equation for the symmetric body-centered solid-on-solid (BCSOS) model, the Kardar-Parisi-Zhang equation for the asymmetric BCSOS model and the generalized restricted solid-on-solid (RSOS) model, and the Villain–Lai–Das Sarma equation for the conserved RSOS model. In addition to the consistent forms of equations for growth models, we also obtain the coefficients associated with the SDEs.

In recent years, the study of nonequilibrium surface growth has attracted considerable interest in both analytical and computational physics. A number of discrete growth models and continuum stochastic equations have been proposed to describe the kinetic roughening properties of surface growth. By studying these models and equations, we classify them into universality classes according to their scaling behavior and associate the continuum stochastic equations with the given discrete growth models.

In general, two methods have been widely used to establish the correspondence between a continuum growth equation and a discrete growth model. One method is to use Monte Carlo simulations to obtain the scaling exponents from the discrete model and compare them with those of the corresponding continuum equation. The other is to derive the continuum equation analytically from a given discrete model. Computational methodologies have contributed significantly to our understanding of epitaxial growth over the past few years and continue to do sounabated. Analytic derivations include methods using the principle of symmetry or reparametrization invariance, and approaches starting from the master equation. In particular, a systematic method proposed by Vvedensky et al. has been successfully applied to the derivation of the continuum growth equations directly from the growth rules of the discrete model for several solid-on-solid discrete models. The derivation procedure of Vvedensky et al. consists of two steps. First, the discrete stochastic equation is derived for the discrete growth model beginning with the master-equation description of the microscopic dynamics of the discrete model. Second, the discrete equation is transformed into a continuous stochastic equation via regularization by expanding the nonanalytic quantities and replacing them with analytic quantities. In this regularization procedure, the step function is approximated by an analytic shifted hyperbolic tangent function, which is expanded in a Taylor series. As pointed out by Předota and Kotrla, the choice of regularization scheme for the step function is ambiguous. Thus, the coefficients in the derived continuum stochastic equation cannot be determined uniquely.

In this Rapid Communication, we present a method for deriving the continuum stochastic equations from the discrete growth models. Our method can be applied to most of the models that are accessible via the method of Vvedensky et al. In addition to the derivation of stochastic equations consistent with the numerical solutions for the discrete models, our method predicts the coefficients in the stochastic equations. Our method begins with mapping of the discrete models onto reaction-diffusion systems with hard-core particles and sets up the master equation of the microscopic dynamics in the form of the Schrödinger equation. Next, borrowing the method introduced in Ref., we derive the corresponding Fokker-Planck equation, and then the stochastic differential equation is obtained. We apply our method to three discrete growth models: the body-centered solid-on-solid (BCSOS) model, the generalized restricted solid-on-solid (RSOS) model, and the conserved RSOS (CRSOS) model.

The BCSOS model is one of the simplest microscopic growth models. Consider a surface built from square bricks rotated by $`\pi /4`$, without defects such as overhangs and vacancies. Each bond $`i`$ contains a step $`S_i=\pm 1`$. The growth dynamics is as follows: Choose one of the columns at random. If this column is at the bottom of a local valley ($`S_{i1}=1`$ and $`S_i=+1`$), a particle adsorbs with probability $`p`$ (and nothing happens with probability $`1p`$). If it is at the top of a local hill ($`S_{i1}=+1`$ and $`S_i=1`$), a particle desorbs with probability $`q`$ (and nothing happens with probability $`1q`$). If it is part of a local slope ($`S_{i1}=S_i`$), nothing happens. This model can be mapped onto the problem of the asymmetric exclusion process (ASEP). The process describes particles that hop independently with hard-core exclusion along a one-dimensional lattice with a bias that mimics an external driving force. By denoting the descending slope as the particle ($`A`$) and the ascending slope as the vacancy($`\mathrm{}`$), the microscopic growth dynamics can be mapped to a stochastic dynamical rule of the ASEP: the only transitions allowed for the site with neighboring bonds ($`i1`$,$`i`$) are

$`A\mathrm{}\mathrm{}A`$ $`\mathrm{with}\mathrm{probability}p,`$ (1)
$`\mathrm{}AA\mathrm{}`$ $`\mathrm{with}\mathrm{probability}q.`$ (2)
This model has been studied extensively in the literature by Monte Carlo simulations and more recently it was realized that it can be solved exactly. Furthermore, Derrida and Mallick calculated the diffusion constant associated with fluctuations of the current in the limit of $`pq`$ and derived the corresponding continuum stochastic equation with the coefficients correct up to the lowest order. Krug et al. found the coefficients for the totally asymmetric, that is, $`q=0`$, case. Here, we apply our method to the ASEP to derive the corresponding stochastic equation. For simplicity, we set $`p=1`$ and $`q=x`$ ($`0x1`$). For $`x=1`$ the system is symmetric, whereas for $`x=0`$ it reduces to the totally asymmetric case.

Introducing the annihilation and creation operators ($`\widehat{a}_i`$’s and $`\widehat{a}_i^{}`$’s, respectively) satisfying the mixed commutation relations explained in Ref. and defining the state vector $`|\psi ;t_CP(C;t)|C`$, the master equation can be written as a Schrödinger-like equation:

$$\frac{}{t}|\psi ;t=\widehat{H}|\psi ;t,$$
(3)
where $`P(C;t)`$ is the probability for the system to be in a given microscopic configuration $`C`$ at time $`t`$, and $`\widehat{H}`$, called the Hamiltonian, is an evolution operator expressed in terms of $`\widehat{a}_i`$’s and $`\widehat{a}_i^{}`$’s. From now on, any operator will have a caret (e.g., $`\widehat{a}`$,$`\widehat{b}`$) and any symbol without a caret should not be confused with an operator. Occasionally, the same symbol is used to represent an operator and a density (e.g., $`\widehat{a}_i`$ and $`a_i`$). The Hamiltonian generating the time evolution of the BCSOS model is found to be

$$\widehat{H}=\underset{i}{}\left(\widehat{a}_i\widehat{a}_{i+1}^{}+x\widehat{a}_i^{}\widehat{a}_{i+1}\right).$$
(4)
Since the diagonal terms that make $`\widehat{H}`$ stochastic have no role in our formalism, for simplicity, we omit these terms here and throughout this Rapid Communication. By involving the commutation relations between the Hamiltonian and some relevant operators such as $`\widehat{a}_i^{}\widehat{a}_i`$ and $`\widehat{a}_i^{}\widehat{a}_i\widehat{a}_j^{}\widehat{a}_j`$, and using the property of the projection state $`|(\widehat{a}_i^{}+\widehat{a}_i)=|`$, where $`|_CC|`$, we find the Kramers-Moyal coefficients corresponding to the above Hamiltonian:

$`C_i`$ $`=`$ $`a_{i+1}+a_{i1}2a_i(1x)\left[a_{i+1}(1a_i)a_i(1a_{i1})\right],`$ (5)
$`C_{ij}`$ $`=`$ $`\left[a_i+xa_{i+1}(1x)a_ia_{i+1}\right](\delta _{i+1,j}\delta _{ij})`$ (7)
$`+\left[a_{i1}+xa_i(1x)a_ia_{i1}\right](\delta _{i,j}\delta _{i1j}).`$
Notice the absence of the caret on the $`a`$’s. From the coefficients $`C_i`$ and $`C_{ij}`$, we write down the discrete stochastic equation

$$\frac{a_i}{t}=C_i+\xi _i,$$
(8)
where $`a_i`$ is the local density of the particle at $`i`$ and $`\xi _i(t)\xi _j(t^{})=C_{ij}\delta (tt^{})`$. Replacing $`a_i`$ by $`(1h)/2`$ and taking the continuum limit, we obtain the continuum stochastic equation for the height variable $`h`$ for the BCSOS model (and equivalently for the ASEP):

$$\frac{h}{t}=\frac{1x}{2}+\frac{1+x}{2}^2h\frac{1x}{2}(h)^2+\xi (r;t),$$
(9)
with $`\xi (r;t)\xi (r^{};t^{})=4(1+x)\rho (1\rho )\delta (rr^{})\delta (tt^{})`$ ($`\rho `$ is the stationary-state density, which is the same as the initial density). For the symmetric process ($`x=1`$), the corresponding equation is the EW equation, whereas for the asymmetric process ($`x1`$) it is the KPZ equation. (Compare the coefficients with those in Refs..)

Next we apply our method to the generalized RSOS growth model that was introduced by Neergaard and den Nijs. They studied this model using an elegant mean-field type approach and derived the deterministic part of the KPZ equation. Our method is able to produce the stochastic (noise) part of the KPZ equation as well as the same deterministic part as the method of Neergaard and den Nijs. The RSOS growth model describes the growth of simple cubic surfaces in which only monatomic steps are allowed. The heights at the nearest-neighbor columns can differ by only $`\mathrm{\Delta }h=0,\pm 1`$ (the RSOS constraint). After choosing one of the columns at random, one particle can be deposited at the site $`i+1/2`$ with a probability 1 according to the height differences at $`i`$ and $`i+1`$. Employing the same parameter as in Ref., we map this model onto the two-species hard-core particle system with the following processes:

$`\mathrm{}\mathrm{}AB:p_h,`$ $`\mathrm{}\mathrm{}BA:q_v,`$ (10)
$`AB\mathrm{}\mathrm{}:q_h,`$ $`BA\mathrm{}\mathrm{}:p_v,`$ (11)
$`\mathrm{}AA\mathrm{}:p_s,`$ $`\mathrm{}BB\mathrm{}:q_s,`$ (12)
$`A\mathrm{}\mathrm{}A:q_s,`$ $`B\mathrm{}\mathrm{}B:p_s,`$ (13)
where a particle $`A`$ ($`B`$) stands for the ascending (descending) bond and a vacuum $`\mathrm{}`$ represents the flat bond. The corresponding Hamiltonian is found to be

$`\widehat{H}={\displaystyle \underset{i}{}}[p_h\widehat{a}_i^{}\widehat{b}_{i+1}^{}+q_v\widehat{b}_i^{}\widehat{a}_{i+1}^{}+q_h\widehat{a}_i\widehat{b}_{i+1}+p_v\widehat{b}_i\widehat{a}_{i+1}`$ (14)
$`+p_s\widehat{a}_i^{}\widehat{a}_{i+1}+q_s\widehat{a}_i\widehat{a}_{i+1}^{}+q_s\widehat{b}_i^{}\widehat{b}_{i+1}+p_s\widehat{b}_i\widehat{b}_{i+1}^{}].`$ (15)
Following the same steps as above, we obtain the Kramers-Moyal coefficients $`C_i^\alpha `$ and $`C_{ij}^{\alpha \beta }`$ ($`\alpha ,\beta `$ indicate either $`A`$ or $`B`$) as follows :

$`C_i^A`$ $`=`$ $`{\displaystyle \frac{p_h+q_v}{2}}v_i(v_{i+1}+v_{i1})+{\displaystyle \frac{p_hq_v}{2}}v_i(v_{i+1}v_{i1})q_ha_ib_{i+1}p_vb_{i1}a_i`$ (18)
$`+{\displaystyle \frac{p_s+q_s}{2}}\left[a_{i+1}+a_{i1}2a_i+a_i(b_{i+1}+b_{i1})b_i(a_{i+1}+a_{i1})\right]`$
$`+{\displaystyle \frac{p_sq_s}{2}}\left[(a_{i+1}a_{i1})(12a_i)(a_{i+1}a_{i1})b_ia_i(b_{i+1}b_{i1})\right],`$
$`C_{ij}^{AA}`$ $`=`$ $`\{{\displaystyle \frac{p_h+q_v}{2}}v_i(v_{i+1}+v_{i1})+{\displaystyle \frac{p_hq_v}{2}}v_i(v_{i+1}v_{i1})+q_ha_ib_{i+1}+p_vb_{i1}a_i`$ (23)
$`+{\displaystyle \frac{p_s+q_s}{2}}\left[v_i(a_{i+1}+a_{i1})+a_i(v_{i+1}+v_{i1})\right]`$
$`+{\displaystyle \frac{p_sq_s}{2}}[v_i(a_{i+1}a_{i1})+a_i(v_{i1}v_{i+1})]\}\delta _{ij}`$
$`\left\{p_s(1a_ib_i)a_{i+1}+q_sa_i(1a_{i+1}b_{i+1})\right\}\delta _{i+1,j}`$
$`\left\{q_s(1a_ib_i)a_{i1}+p_sa_i(1a_{i1}b_{i1})\right\}\delta _{i1,j},`$
$`C_{ij}^{AB}`$ $`=`$ $`C_{ji}^{BA}=\delta _{i+1,j}(p_hv_iv_{i+1}+q_ha_ib_{i+1})+\delta _{i1,j}(q_vv_iv_{i1}+p_va_ib_{i1}),`$ (24)
where $`v_i=1a_ib_i`$ and $`C_i^B`$ ($`C_{ij}^{BB}`$) is obtained by the exchange $`ab`$ followed by $`i+kik`$ in $`C_i^A`$ ($`C_{ij}^{AA}`$). By introducing the local slope $`Dab`$ and the step density $`Sa+b`$, we derive the deterministic parts of the stochastic equations for these two parameters:

$`{\displaystyle \frac{D}{t}}`$ $`=`$ $`\left[f_q(1S)^2+s_d(1S)S+{\displaystyle \frac{1}{4}}h_g(S^2D^2)\right]`$ (25)
$`+`$ $`{\displaystyle \frac{1}{2}}s_m^2D+{\displaystyle \frac{1}{4}}(2s_ma_s)\left[DSSD\right],`$ (26)
$`{\displaystyle \frac{S}{t}}`$ $`=`$ $`2c_s(1S)^2{\displaystyle \frac{1}{2}}a_s(S^2D^2)+s_d\left[(1S)D\right]`$ (27)
$`+`$ $`{\displaystyle \frac{1}{2}}h_g[DSSD]c_s(1S)^2S`$ (28)
$`+`$ $`{\displaystyle \frac{1}{2}}s_m^2S{\displaystyle \frac{1}{4}}a_s\left[S^2SD^2D\right],`$ (29)
where we have used the same notation for the parameters as in Ref.; $`c_s=p_h+q_v`$, $`a_s=p_v+q_h`$, $`s_m=p_s+q_s`$, $`f_g=p_hq_v`$, $`h_g=p_vq_h`$, and $`s_d=p_sq_s`$. The equation for the step density $`S`$ is the same as that of Neergaard and den Nijs. As pointed out in Ref., the equation for the step density contains a mass term and the step density reaches its stationary value $`S_0=\left[1+\sqrt{a_s/(4c_s)}\right]^1`$ after a characteristic time $`\tau _s=(2\sqrt{a_sc_s})^1`$. Thus, although there are two order parameters, only the local slope fluctuates at time scales larger than $`\tau _s`$. At larger time scales, the step density does not behave as an independent dynamic variable. It follows local fluctuations in the slope of the surface:

$$S=S_0+\frac{1}{2}\tau _sa_sD^2+\tau _s\left[s_d(1S_0)\frac{1}{2}h_gS_0\right]D+\mathrm{}.$$
(30)
Substituting Eq. (30) for $`S`$ and identifying $`D=h`$, the equation for the height variable becomes

$$\frac{h}{t}=v_{\mathrm{}}+\nu ^2h+\frac{\lambda }{2}(h)^2+\xi ,$$
(31)
where

$`v_{\mathrm{}}`$ $`=`$ $`f_g(1S_0)^2+s_d(1S_0)S_0+{\displaystyle \frac{1}{4}}h_gS_0^2,`$ (32)
$`\nu `$ $`=`$ $`{\displaystyle \frac{1}{4}}\tau _s[2s_d(1S_0)h_gS_0][h_gS_04f_g(1S_0)+2s_d(12S_0)]+{\displaystyle \frac{1}{2}}s_m(1S_0)+{\displaystyle \frac{1}{4}}a_sS_0,`$ (33)
$`\lambda `$ $`=`$ $`{\displaystyle \frac{1}{2}}h_g+{\displaystyle \frac{1}{2}}\tau _sa_s[h_gS_04f_g(1S_0)+2s_d(12S_0)],`$ (34)
and $`\xi (x,t)\xi (x^{},t^{})=D_{\xi \xi }\delta (xx^{})\delta (tt^{})`$ with $`D_{\xi \xi }=c_s(1S_0)^2+s_mS_0(1S_0)+a_sS_0^2/4`$. This equation is the KPZ equation corresponding to the general RSOS model. To compare these coefficients with numerical work, let us consider the simple RSOS model introduced by Kim and Kosterlitz (KK model). The KK model corresponds to $`p_h=p_v=p_s=1`$ and $`q_h=q_v=q_s=0`$. We obtain the corresponding coefficients $`v_{\mathrm{}}=4/9`$, $`AD_{\xi \xi }/(2\nu )=2/3`$, $`\lambda =5/6`$. These values are consistent with the estimated values from a numerical study by Krug et al..

Recently, a different growth model with a RSOS condition has been proposed and studied by Kim et al.. Instead of rejecting the particle when the RSOS condition is not satisfied, this model allows the deposited particle to hop to the nearest site where the RSOS condition is satisfied. Thus, this model has the constraint of a conserved growth condition and is called the conserved RSOS model. The detailed derivation of Villain–Lai–Das Sarma (VLD) equation from the CRSOS model will be published elsewhere. Here we only sketch the procedure and report the result. The procedure is similar to the one used to get the KPZ equation from the RSOS model except for some complicated calculation. After some algebra, we find the VLD equation

$$\frac{h}{t}=\stackrel{~}{\nu }^4h+\stackrel{~}{\lambda }^2(h)^2+\eta ,$$
(35)
where $`\stackrel{~}{\nu }=(2112\sqrt{2})/2`$, $`\stackrel{~}{\lambda }=(103\sqrt{2})/2`$, and $`\eta (x,t)\eta (x^{},t^{})=D_{\eta \eta }\delta (xx^{})\delta (tt^{})`$ with $`D_{\eta \eta }=(2\sqrt{2}1)/2`$. In deriving the above equation, we kept only the most relevant terms, and found that neither the EW nor the KPZ term exists in the growth equation. It is known, however, that higher order terms of the form $`(h)^{2n+1}`$ $`(n1)`$ generate the EW term by the dynamic renormalization group. Hence we should investigate the possibility of occurrence of these terms. Indeed, we found that the dangerous term of the form $`(h)^{2n+1}`$ does not arise in the derivation of the VLD equation. Consequently, we concluded that the continuum equation of the CRSOS model is the VLD equation.

Although the VLD equation was derived by Huang and Gu using the master-equation description with the regularization procedure proposed by Vvedensky et al., there were some ambiguities in choosing the regularization of the step function and thus the coefficients could not be predicted, whereas, in our derivation, we are able to predict the coefficients for the VLD equation corresponding to the CRSOS model. Unfortunately, however, there is no numerical method, to our knowledge, to find $`\stackrel{~}{\nu }`$ and $`\stackrel{~}{\lambda }`$ for a microscopic model. Numerical studies up to now can argue only that $`\stackrel{~}{\lambda }`$ may be positive, which is consistent with our derivation.

In deriving the VLD equation, we recognized the interesting aspect of the RSOS model. When we allow hopping processes only up to distance $`l_0`$, the stochastic equation for the height variable eventually becomes the KPZ one. This is contradictory to previously reported simulation results. Kim and Yook studied the RSOS model with finite-distance hopping by Monte Carlo simulation and concluded that there is a phase transition at finite $`l_0`$ from the KPZ class ($`l_0=0`$) to the VLD class ($`l_0=\mathrm{}`$). However, this conclusion seems to be a finite-size effect. In fact, we have confirmed our argument by carrying out a Monte Carlo simulation with sufficiently large system size for several $`l_0`$.

In summary, we have presented a formalism for deriving the continuum stochastic differential equations corresponding to discrete growth models. Applying the formalism to the BCSOS model, we derived the EW equation for the symmetric process and the KPZ equation with exact coefficients for the asymmetric process. The RSOS model was also studied with the general probabilities for possible processes. We derived the KPZ equation with a fluctuating noise part as well as the deterministic part, which is the same as the result of Neergaard and den Nijs. For the special case with $`p`$ $`=`$ 1 and $`q`$ $`=`$ 0 (the KK model), our coefficients are consistent with the numerical results of Krug et al.. Finally, we applied our formalism to the conserved RSOS model. For the CRSOS model, we found that the VLD equation is the corresponding continuum SDE. However, if we allow only finite hopping ($`l_0<\mathrm{}`$), the system belongs to the KPZ class eventually. We also predict the coefficients of the VLD equation for the CRSOS model.

This work was supported by Grant No. 2000-2-11200-002-3 from the BRP program of the KOSEF, and by the Brain Korea 21 Project at Seoul National University.