text:
# The Evolution of Stellar Exponential Discs

## 1 Introduction

One of the most notable properties of disc galaxies is that the surface brightness profiles of their stellar discs are remarkably close to exponential distributions over a large dynamic range in radii (eg. de Vaucouleurs 1959; Freeman 1970; de Jong 1996a). This has traditionally been held to imply an exponential surface density profile for gas infalling to the galactic plane, which in turn places particular constraints on the detailed angular momentum distribution of baryonic material in the protogalaxy. Gunn (1982), for example, showed that, under the assumption of angular momentum conservation for each gas element under infall, a uniformly rotating uniform sphere collapses to an approximately exponential profile, however the resulting function is only exponential over $`3`$ scale lengths compared to the more impressive 4-6 scale lengths exhibited by real galaxies (eg. de Jong 1996a). Calculations with other initial halo density profiles and angular momentum distributions produce similar results (eg. van der Kruit 1987; Dalcanton, Spergel & Summers 1997; Bullock et al 2000).

The exponential profile is so ubiquitous, and its appearance when plotted as magnitude versus radius so deceptively simple, that one may be tempted to regard it as a natural or obvious outcome of gravitational collapse. In fact, gravitational collapse calculations in other astronomical contexts never give rise to exponential disc profiles (Cassen & Moosman 1981; Terebey, Shu & Cassen 1984), with the more usual outcome being profiles of approximately power law form. In the galactic context, the difficulty of producing exponential profiles over a large number of scale lengths through angular momentum conserving collapse from reasonable initial conditions has recently been highlighted by Efstathiou (2000): the outer regions of galactic discs contain so little mass, that for power law halo density profiles, the outer disc material must originate from a small range of radii in the halo. On the other hand, these regions contain a large dynamic range of angular momentum, so that these conditions can only be met simultaneously if the angular momentum profile at the halo suffers an abrupt, and physically implausible, upturn at some radius (see Figure 6(b) of Efstathiou 2000).

Exponential profiles have the unique feature that the scale length on which quantities change is the same at all radii. In the case of a centrifugally-supported disc, the exponential scale length is determined (in a fixed potential) purely by the average angular momentum per unit mass for the entire disc. Such a system thus has the remarkable property that material over a significant range in radii (and with a correspondingly significant range in specific angular momenta) is somehow ‘informed’ about the average specific angular momentum of the disc as a whole. The key question is whether this property is a result of the formation process of the disc, or somehow acquired through subsequent evolution.

Disc formation models have moved on from the picture of smooth dissipational collapse within a dark halo potential, and recent numerical calculations have focused on the more complex formation histories resulting from hierarchical growth (eg. Kauffmann, White & Guiderdoni 1993; Cole et al 1994; Navarro, Frenk & White 1995). A notable outcome of such calculations is the strong redistribution of angular momentum during collapse. Attention has so far mainly focused on the net transfer of angular momentum from baryonic to dark matter and the consequent implications for the overall sizes of galactic discs (eg. Navarro & Benz 1991; Navarro & Steinmetz 1997; Weil, Eke & Efstathiou 1998). A related issue is, however, the redistribution of angular momentum within the baryonic material and its effect on the resulting disc profile. The large dynamic range in size scales that have to be captured in such simulations, and the difficulties in controlling the effective viscosity in SPH calculations that are close to the limit of adequate resolution, mean that it is not yet clear whether discs that are exponential over a number of scale lengths are indeed the expected outcome.

Observational evidence against a picture of exponential infall followed by in situ star formation is provided by the radial profiles of gas in disc galaxies. The gas surface density profiles are generally observed to be much flatter then the stellar profiles. In the inner regions of disc galaxies, the gas is the minority component and represents the residue following star formation: it is not difficult, by appropriate adjustment of the star formation law, to generate gas profiles within the optical disc that are compatible with observations using an initially exponential gas profile (eg. Boissier & Prantzos 1999). In the outer regions, however, the gas is the dominant component and its profile should presumably reflect its infall distribution. It seems highly suspicious, in the context of in situ models, that the gas profiles at large radii are so different from the exponential profiles that characterise the stars at smaller radii: if the stellar exponential profile derives from exponential infall of gas, surely one would expect to see that profile continued in the gas not yet incorporated in stars?

All of these difficulties are circumvented by the simple model first proposed by Lin and Pringle (1987a). These authors envisaged that gas in the disc not only turns into stars (on a timescale $`t_{}`$), but is also subject to some form of effective viscosity that gives rise to angular momentum redistribution and radial flows. The dominant physical process responsible for the viscosity is poorly constrained, but could plausibly result from cloud-cloud collisions and/or gravitational instabilities in the disc. If viscous redistribution occurs on a timescale $`t_\nu `$, then the resultant stellar profile is exponential for a wide range of initially less centrally concentrated gas profiles, provided only that $`t_\nu t_{}`$. This result can be understood, in the context of the above discussion, inasmuch as viscous processes in the disc (provided $`t_{}`$ is not much less than $`t_\nu `$) allow material at different specific angular momenta to communicate with each other, and thus to ‘inform’ themselves of the average specific angular momentum of disc material, which in turn sets the disc scale length. On the other hand, $`t_{}`$ cannot be much greater than $`t_\nu `$ because in that case, viscous evolution would run away too far towards its ultimate endpoint (i.e. ‘all the mass at the origin, all the angular momentum at infinity’, Lynden-Bell and Pringle (1974)) to be compatible with the observed profiles of disc galaxies. Fortunately it turns out that the generation of exponential profiles over a number of scale lengths does not require particularly fine tuning between $`t_{}`$ and $`t_\nu `$, and that satisfactory fits are obtained even if the ratio $`t_{}/t_\nu `$ deviates from unity by half an order of magnitude, or is a mild function of radius (Clarke 1989).

The viscous model thus relieves cosmological models of the burden of generating the ‘right’ angular momentum profile (though, by the same token, it also erases information about the profile on infall, and thus weakens the role of gaseous and stellar disc profiles as cosmological probes). It also predicts that the radial gas surface density distribution beyond the optical radius should be rather flat, in accordance with observations. The viscous model was able to satisfy the various observational constraints available a decade ago (mainly surface density profiles and chemical constraints imposed by the Milky Way), although it was not manifestly superior to conventional in situ models in this respect (Clarke 1989, Yoshii & Sommer-Larsen 1989). Then, as now, the main drawback to viscous models is the lack of a detailed theory linking $`t_{}`$ to $`t_\nu `$, although it is not difficult to invoke physical processes (e.g. self-gravitational instability of the disc) that one might expect to simultaneously promote star formation and angular momentum redistribution in the disc (Lin & Pringle 1987b; Olivier, Primack & Blumenthal 1991).

A wealth of new observational data pertaining to the evolution of disc galaxies has emerged in recent years from both high redshift galaxy surveys – which permit a direct study of how galactic properties evolve with epoch – as well as detailed studies of the fossil stellar populations in our own Milky Way and other local systems. Given that the general arguments in favour of viscous models appear at least as strong now as they did a decade ago, it would seem timely to revisit this model and investigate its predictions in the light of new and forthcoming observations. In this paper, we explore the structural evolution and radial growth of stellar discs subject to simultaneous star formation and viscous evolution by considering models in which the initial gas disc is pre-assembled as well as models in which the gas disc is assembled over an extended period due to the continued infall of material from the protogalactic halo. In both cases, the gaseous and stellar discs are naturally built up from the inside-out, though the mechanism is different in each case: in pre-assembled models, gas is transferred to larger radii only as a result of viscous torques in the disc, whereas in the case of continued infall, the dominant cause of the inside-out growth is the increasing specific angular momentum of infalling material at later times (whether or not viscous flows are also important).

We do not address in this paper how well viscous models can be made to reproduce in detail the present-day properties of the Milky Way and other nearby spirals; for a recent discussion of this aspect see Thon & Meusinger (1998). We also bypass the issue of how one sets up the conditions in the collapsing gas that will reproduce the correct absolute values of the disc scale lengths. The resultant scale length in viscous evolution models is set purely by the average specific angular specific angular momentum of the gas infallen at a given time hence the problem of forming disc systems with large scale lengths (as discussed by Mo et al 1998) is one of cosmological initial conditions and not one that viscous evolution can alleviate or affect in any way. In the present paper, we simply assume that the forming disc has access to material of the requisite high angular momentum. We describe in Section 2 the model and free parameters and present our results in Section 3. Section 4 discusses model constraints from existing and future observational data on both local and distant spirals.

## 2 Model Parameters

We envisage a general scenario for disc galaxy formation in which infalling material gradually builds up the disc over a time interval $`t=0`$ to $`t=t_{infall}`$ and (simultaneous) star formation and viscous redistribution occur over the interval $`t=t_{sf}`$ ($`t_{infall}`$) to $`t=t_{now}`$ ($`t_{infall}`$). The case $`t_{sf}=t_{infall}`$ corresponds to the case of star formation in a pre-assembled gas disc whereas $`t_{now}t_{infall}>>t_{sf}`$ corresponds to the extreme infall case where gas is accreted up until the present epoch.

### 2.1 Infall Calculation

We calculate the surface density profile with which material arrives in the disc plane by assuming an idealised (spherically symmetric) density distribution and angular velocity distribution for baryonic material in the protogalaxy. We then map this material down into the disc plane by invoking detailed angular momentum conservation, neglecting both the self-gravity of the baryonic material and the response of the halo to the infalling baryons. We justify this rigid halo approximation on the grounds that the initial distribution of baryons in the halo is poorly known (and likely to be considerably more complex than the smooth spherically symmetric distribution invoked here) and also because in any case, the behaviour of the viscous evolution is only weakly dependent on the shape of the initial distribution (as opposed to its total mass and angular momentum). We calculate the gas infall rate (for $`tt_{infall}`$) by requiring that all material infalling from within a spherical radius $`r_{inf}(t)`$ has reached the disc plane by a time $`t`$. This infall radius is fixed by the requirement that the free fall time within $`r_{inf}`$ is much less than the instantaneous cosmic shear time ($`H^1`$) - or, equivalently, that the mean enclosed density is much larger ($`200`$) than the mean density of the Universe at that epoch. As successive shells of matter are ‘released’ and allowed to fall in, we map them onto the disc plane (without applying the (small) correction for the finite infall time of each shell). We do not explicitly consider the gas cooling time in our model; numerical simulations suggest that the cooling radius tracks the virial radius for Milky Way-type galaxies hence that all gas within $`r_{inf}`$ is indeed able to infall on a free fall time (Somerville & Primack 1999). Our default halo model assumes that the baryons initially follow an $`r^2`$ density distribution (ie. isothermal sphere), so that for an Einstein-de Sitter cosmology, the total infall rate onto the galaxy is constant in time for $`t<t_{infall}`$. At any given radius however, the infall history is characterised by an initial lag of variable length (corresponding to the time required for a shell of high enough angular momentum to fall in), a peak and then a roughly exponential decline. This prescription thus produces a crude on/off switch on the infall, so that, if $`t_{infall}`$ is comparable to $`t_{now}`$, it implies considerably more infall at recent times and at large radii than is usually assumed in galactic infall models, where the infall is arbitrarily assumed to decline exponentially in both time and radius (eg. Boissier & Prantzos 1999; Chiappini et al 1997) in order to mimic inside-out growth. We motivate the prescription employed here by the wish to explore the response of viscous discs to extreme infall conditions.

Finally, in order to prescribe the surface density distribution generated by each shell as it falls in, we must assign an initial angular momentum distribution for the baryons in the halo. The angular momentum of the dark matter and baryons is presumably acquired through tidal torques from neighbouring structures (Peebles 1969) and thus the average angular momentum of infalling material (whether the infall be smooth or discrete) is expected to increase with time. We assume the halo baryons initially rotate rigidly on spheres with angular velocity $`r^1`$, which implies that the mean value of the spin parameter $`\lambda `$ for each shell is constant. For a power law angular velocity law, the surface density distribution generated by each shell is self-similar. We scale the overall dimensions and normalisation of the system to match the observed profile of the Milky Way. As stated before, our calculation thus sidesteps the difficulty of how, given the usual assumptions about the value of $`\lambda `$ and the size scale for collapse as a function of redshift, it is possible to produce large, Milky Way-scale, discs at redshifts greater than $`1`$ (eg. Mo et al 1998).

Figure 1 shows the final gas surface density profile at the completion of the infall phase (hence the initial starting point in the case of a pre-assembled gas disc). Beyond the central region, the gas surface density slowly declines and then turns over abruptly at the maximum infall radius, $`r_{init}`$. The value of $`r_{init}`$/RD depends primarily on the initial halo profile and rotation law, and to a lesser extent on the star formation and viscosity prescriptions. If the disc surface density decreases outwards throughout the disc, $`r_{init}>3R_D`$; the absolute minimum value of $`r_{init}`$, which corresponds to an annulus of matter at this location, is $`2R_D`$. This latter configuration does not viscously evolve to a smooth exponential profile however. For the flat rotation curve and halo profile assumed here (and neglecting disc self-gravity), $`r_{init}5R_D`$. An exponential fit (dashed line) provides a good description of the radial behaviour of the initial gas profile over $`2.5`$ scale lengths.

### 2.2 Disc Evolution

The evolution of gas surface density, $`\mathrm{\Sigma }_{gas}(r,t)`$ in a disc that is subject to simultaneous viscous redistribution and star formation is described by:

$$\frac{\mathrm{\Sigma }_{gas}}{t}=\frac{1}{r}\frac{}{r}\left\{\frac{/r[\nu \mathrm{\Sigma }_{gas}r^3(d\mathrm{\Omega }/dr)]}{d/dr(r^2\mathrm{\Omega })}\right\}\frac{\mathrm{\Sigma }_{gas}}{t_{}}$$
(1)
where $`r,\mathrm{\Omega },\nu `$ and $`t_{}`$ represent the galactocentric radius, angular velocity, kinematic viscosity and local star formation timescale. The evolution of the stellar surface density $`\mathrm{\Sigma }_{}(r,t)`$ is then given by:

$$\frac{\mathrm{\Sigma }_{}}{t}=\frac{\mathrm{\Sigma }_{gas}}{t_{}}$$
(2)
The gas surface density in these models is assumed to be the sum of the neutral and molecular components. We calculate the evolution of the gaseous and stellar components of the disc, for times $`t>t_{sf}`$, assuming that the gas evolves in the fixed potential of a halo with flat rotation law (i.e. $`\mathrm{\Omega }r^1`$). This simplification will not be correct for the very inner regions of the disc, and so we do not attempt any detailed comparison of our results with structures in the central parts of disc galaxies.

In order to obtain a resulting exponential stellar profile, we follow Lin & Pringle (1987a) and fix the relationship between the star formation and viscosity prescriptions such that the characteristic timescale for material to move a fractional radius of order unity is comparable with the timescale on which it is converted into stars (i.e. $`t_{}t_\nu `$). With this constraint, and given an assumed rotation law, the only free parameter is the star formation prescription. Our default star formation law relates the star formation rate per unit area with the gas column density, $`\mathrm{\Sigma }_{gas}`$, as $`\mathrm{\Sigma }_{gas}R^1`$, which corresponds to an $`n=1`$ Schmidt law, modulated by the local angular frequency (eg. Wyse & Silk 1989). A prescription of this sort (with power-law $`n=1.5`$) has recently been shown to be in good agreement with the observed radial star formation rate in the Milky Way (Boissier & Prantzos 1999). We also consider a model with a dependence of $`\mathrm{\Sigma }_{SFR}\mathrm{\Sigma }_{gas}^{1.5}`$ as well as a model in which the star formation law steepens consderably in the outer disc, $`\mathrm{\Sigma }_{SFR}(R>3R_D)\mathrm{\Sigma }_{gas}/R^2`$. While the former has been shown to provide a reasonable description of the star formation rates observed in the inner regions of galactic discs, the two component law may provide a better description of the steep declines in star formation rate observed at large radii, perhaps resulting from the flaring of the gas disc in these parts (eg. Ferguson et al 1998). We stress that in this paper we are seeking generic properties of viscous models, and their inter-relation with continued infall, that do not depend strongly on the star formation prescription employed. While varying the star formation prescription also requires changing the form of the viscosity law so that the approximate equality in the respective timescales is preserved, the dominant viscosity-producing mechanism in discs is so poorly constrained that this does not create any obvious physical problem111For a given star formation law of the form $`\mathrm{\Sigma }_{SFR}\mathrm{\Sigma }_{gas}^aR^b`$, the corresponding viscosity prescription is $`\nu \mathrm{\Sigma }_{gas}^{a1}R^{2b}`$.

We compute the evolution of a viscous star forming disc subject to infall using a standard first order explicit scheme with $`200`$ gridpoints, equispaced in $`R`$ with $`R_{out}/R_{in}=150`$. Zero torque boundary conditions are applied at the inner and outer edge of the grid. In practice, we prescribe $`t_{sf}`$ and $`t_{infall}`$ and then pursue the calculation until a time $`t_{now}`$ at which the system is deemed to resemble the Milky Way - specifically, we require that the gas fraction at a radius of $`2.7`$ exponential scale lengths222We adopt R$`{}_{}{}^{}=8.0`$ kpc (Reid 1993) and RD=3.0 kpc (Sackett 1997). is $`0.4`$.

## 3 Results

### 3.1 The Case of a Pre-Assembled Gas Disk

We first consider the results of models in which star formation (hence viscous redistribution) switches on after the initial gas disc is already in place – ie. pre-assembled. This corresponds to the scenario where all or most of the infall phase is complete before the onset of star formation ($`t_{sf}t_{infall}`$). Figure 2 shows snapshots of the gaseous and stellar surface density distributions at various fractional times since the onset of star formation for three different star formation prescriptions ($`\mathrm{\Sigma }_{SFR}\mathrm{\Sigma }_{gas}/R`$, $`\mathrm{\Sigma }_{SFR}\mathrm{\Sigma }_{gas}^{1.5}`$ and a $`\mathrm{`}`$two component’ law described by $`\mathrm{\Sigma }_{SFR}\mathrm{\Sigma }_{gas}/R^n`$, where $`n=1`$ out to 3 scale lengths and $`n=2`$ beyond that). The profiles are normalised by matching the final (ie. present day) gas and stellar surface density in the solar neighbourhood to their measured values (we assume $`\mathrm{\Sigma }_{gas}=14.5`$M/pc<sup>2</sup> (Olling & Merrifield 1998), $`\mathrm{\Sigma }_{stars}=40`$M/pc2 (Gould et al 1996)). Figure 2 illustrates several notable features. Firstly, different star formation prescriptions lead to very similar exponential behaviours (over at least 5 scale lengths) for the stellar surface density profiles; models really only deviate from each other at early times and at large (unobservably faint) radii. Interestingly, the two component model produces rather smooth profiles despite the abrupt change in the form of the star formation and viscosity laws at roughly three scalelengths. These particular profiles deviate somewhat from exponential behaviour beyond 5–6 scale lengths, but the effect is subtle. Gas profiles are more sensitive to the different star formation laws, but generally exhibit exponential declines in their outer regions with e-folding scales significantly larger than that of the stars. At late times, the scale length of the gas at large radii, where the gas is still the majority component, exceeds the stellar scale length by a factor of approximately $`2`$. Secondly, we find the stellar component exhibits self-similar growth out to $`5`$ scale lengths over the entire star forming period, implying an invariance in the stellar scale length with time, as first noted by Olivier et al (1991). This is further quantified in Figure 3 which shows that the scale length lies well within $`20\%`$ of its present-day value throughout the disc’s history. Thus, while viscous discs are characterised by inside-out growth, which has led to the speculation that the resulting scale length should increase with time (eg. Zhang & Wyse 2000), we here show that this growth is simply manifest as an expansion in the radial extent of the exponential stellar disc due to the viscous transport of gas from smaller radii. This effect is unlikely to be verifiable by studies of high redshift galaxies, however, since the strongest evolution occurs in regions of the disc so faint that they are inaccesible to direct observation.

Figure 4 shows the star formation rate (dashed lines) normalised to the peak rate as a function of time at $`2.7`$ (‘solar neighbourhood’), $`5`$ and $`10`$ exponential scale lengths for the particular star formation law $`\mathrm{\Sigma }_{SFR}\mathrm{\Sigma }_{gas}/R`$; other star formation prescriptions exhibit similar behaviours. For a Freeman exponential disc ($`\mathrm{\Sigma }_B=21.65`$ mag/$`\mathrm{}^{^{\prime \prime }}`$), 5 scale lengths corresponds to an expected surface brightness of $`1\%`$ of sky, or about the current limit for quantitative studies of surface brightness in external galaxies. The top axis shows age of the disc in Gyr, where we have assumed the present-day stellar disc age of the Milky Way to be 11 Gyr as determined from the ages of the oldest stars in the Hipparcos sample (Binney, Dehnen & Bertelli 2000). The disc at 2.7 and 5 scale lengths exhibits slowly declining star formation rates in the pre-assembled model, decreasing by factors of roughly $`5`$ and $`2`$ over the disc’s lifetime respectively. On the other hand, the disc at 10 exponential scalelengths is characterised by a slowly increasing star formation rate (bottom panel Figure 4) reflecting the finite time required to transport gas outwards to these parts by viscous torques; the arrival of material in the extreme outermost disc is therefore delayed long after the onset of star formation at smaller radii. This delay is not a major effect at $`5`$ scale lengths however, even for the most compact plausible initial gas distributions.

Also shown in Figure 4 is the fraction of the final stellar mass in place as a function of time (solid lines). These curves are contrasted with the analytic predictions of a simple in situ model (dashed-dotted line) obeying the same star formation law (note that in the in situ model the total - gas $`+`$ stars - radial profile is exponential by construction at all times) calculated as:

$$\mathrm{\Sigma }_{SFR}=\frac{\mathrm{\Sigma }_{gas}}{t}=\frac{k\mathrm{\Sigma }_{gas}}{R}$$
(3)
where the adopted star formation efficiency, $`k`$, is chosen to be consistent with the viscous model. Integration yields

$$\frac{\mathrm{\Sigma }_{stars}(t)}{\mathrm{\Sigma }_{stars}(t_{now})}=\frac{1e^{kt/R}}{1e^{kt_{now}/R}}$$
(4)
where we have normalised the mass in stars present at a given time to that of the present day. While the mean age of the stellar population present at $`5`$ scale lengths is somewhat smaller than at $`2.7`$ scale lengths, the effect is rather small and the bulk of the stars at this radius have formed fairly early on. Specifically, while 50% of the final stellar surface density is in place after $`0.3`$ of the star formation duration (or $``$ 3.5 Gyr for our assumed disc age) at $`2.7`$ scale lengths, this same amount is in place after $`0.4`$ of the star formation period (or $``$ 5 Gyr) at $`5`$ scale lengths. Interestingly, there does not appear to be any great difference between the growth rate of stellar mass at $`5`$ scale lengths in the viscous and in situ models, implying that viscous flows have little impact on the star formation history at these radii; observations of stellar age distributions alone in these parts could not distinguish between viscous and in situ scenarios. Beyond the radius of the initial gas disc, the stellar population becomes increasingly young with radius.

### 3.2 The Case of Extreme Infall

We now contrast the models above with viscous models in which $`t_{infall}t_{now}`$. We refer to this as the $`\mathrm{`}`$extreme infall’ case, since the total mass infall rate onto the galaxy has been constant up to nearly the present day, and star formation in the disc started early on in the infall period ($`t_{star}\frac{1}{10}t_{infall}`$). Figure 5 shows the stellar and gaseous surface density profiles for different star formation prescriptions at various times since the onset of star formation. For all models, the final stellar and gaseous profiles are almost identical for the extreme infall viscous model and the pre-assembled viscous model (contrast the highest stellar profiles and lowest central gas profiles in Figures 2 and 5). The temporal evolution of these profiles is very different however. While the disc is able to continuously adjust to the changing profile of the infallen material so as to generate an approximately exponential stellar profile at all times, the scale length of the exponential evolves systematically with time, simply reflecting the mean specific angular momentum of the material which has fallen in by that epoch. The stellar profiles are always exponential over at least 4–5 scale lengths. For the density and angular momentum distributions in the halo assumed here ($`\rho r^2,\mathrm{\Omega }r^1`$), the total angular momentum scales with the square of the infallen mass, so that the mean specific angular momentum is simply proportional to the infallen mass. Our default infall prescription, which makes specific assumptions about the halo properties and the cosmological model (see Sec 2.1), implies a constant infall rate for $`tt_{infall}`$ and hence the mean specific angular momentum of the infallen material rises linearly with time over this period. Figure 6 demonstrates that the stellar scale length indeed displays this behaviour.

A fundamental assumption in our infall model is that the mean specific angular momentum of the accreted material systematically increases with time. We feel this is a valid assumption, regardless of the exact nature of the infall, and should be preserved. For example, if the dominant mode of infall is high velocity cloud accretion (Blitz et al 1999), one would still expect those clouds falling in at late times to have originated further from the Milky Way and to have suffered on average more tidal torquing than clouds which fell in at early times. While individual clouds will carry a range of angular momenta (and indeed in our models a given shell maps onto a range of radii), the mean trend should be for the average angular momentum to increase with time. It is worthwhile investigating how much of the scale length evolution is due to other aspects of our adopted infall model however. Figure 7 shows the growth of the stellar exponential scale length for an isothermal halo model in which the total infall rate is allowed to decrease exponentially with time. In this case, the amount of recent infall (ie. t $``$ tnow) is reduced relative to that in the constant infall case by a factor $``$ 11 Gyr/$`\tau _{inf}`$, where 11 Gyr is the present-age of the disc and $`\tau _{inf}`$ the e-folding time of the infall rate. In the limit of $`\tau _{inf}>>11`$ Gyr, this infall prescription mimics the constant infall case whereas for $`\tau _{inf}<<11`$ Gyr it mimics the pre-assembled disc case, with essentially all the infall occuring early on; Figure 7 shows an intermediate case where $`\tau _{inf}=`$5 Gyr. We also examine an infall model in which the baryons in the halo are described by a Navarro-Frenk-White (NFW) profile in which $`\rho (r/r_s)^1(1+r/r_s)^2`$ with the parameter rs defining the scale radius. Navarro, Frenk & White (1997) argue on the basis of high resolution N-body simulations that this form better characterizes the equilibrium density profiles of dark matter halos than an isothermal sphere. In terms of infall rates, it leads to a somewhat higher infall rate than the isothermal sphere at early times and somewhat lower at late times. Figure 7 indicates that the growth of the scale length in these alternate infall models is very similar to that in the constant infall case and thus that our results on the scale length evolution are robust. Contrasting Figures 6 and 7 with Figure 3, one concludes that if clear evidence were to emerge for an increase in stellar scale length of galaxies over cosmological time, it would, within the framework of our model, be an unambiguous indicator of the late infall of high angular momentum material and could not merely be attributed to the effects of viscous evolution.

A consistent picture emerges from analysis of the star formation history at different locations in the disc. Figure 8 shows the star formation rate (dashed line) normalised to the peak rate as a function of time for our default star formation law. Very few stars form initially at the solar neighbourhood as a direct result of the low gas density (only gas viscously transported from smaller radii is present). As the peak infall rate at this radius is reached (indicated by the arrow in the top panel of Figure 8), the star formation rate rapidly increases and then levels off to an almost constant value until the present epoch. The fraction of the final stellar mass in place as a function of time for singular isothermal sphere and NFW halo baryon distributions is also shown (solid and dashed-dotted lines respectively). As expected, the disc at $`2.7`$ scale lengths is predominantly young in the extreme infall model, with only a small dependence on the specific form of the infall prescription. Specifically, 50% of the final stellar mass in place after 0.7 of the star forming duration, or $`8`$ Gyr if a disc age of 11 Gyr is assumed (note that in this case very few 11 Gyr stars are expected to exist outside the innermost regions of the disc, where the first infall is received). At $`5`$ and 10 scale lengths, the stellar population is conspicuously young in the extreme infall model. The disc at 5 scale lengths receives material from direct infall of halo gas only at the very end of the infall phase, so that prior to this, gas only reaches this location as a result of outward viscous transport from an initially much more compact gas distribution. This is especially true of the disc at 10 scale lengths, which never receives direct infall in our model. Consequently, the dominant star-forming phase in these regions only gets under way at late times. We draw attention to the fact that this combined infall/viscous model may produce an unusual chemical signature in the very outer regions, comprising an older metal rich population (formed from pre-enriched gas transported outwards from small radii at early times) and a young metal poor population containing stars formed from recent unenriched infall. We plan to investigate detailed predictions for the chemical signatures of viscous evolution and infall in a future paper.

## 4 Discussion

We have explored models for the evolution of disc galaxies subject to star formation, viscous redistribution of angular momentum, and cosmologically-motivated gaseous infall. We have focused attention on the extreme cases of these models and the rather different predictions which arise for the temporal evolution of the exponential scale length, and the radial growth rate of galactic stellar discs. In this section, we discuss how these model predictions compare with relevant observational data for both local and distant disc galaxies.

### 4.1 Constraints from the Local Universe

Stellar populations in the Milky Way and other nearby spirals provide the $`\mathrm{`}`$fossil record’ of star formation and galaxy assembly over cosmic lifetimes. The wealth of high quality data available for stars in our local neighbourhood of the Milky Way disc makes it a logical starting point for model comparisons.

An important prediction of our models is the star formation history – which is reflected in the stellar age distribution – at various radii. The pre-assembled and extreme infall models make very different predictions for the star formation history of the disc at 2.7 exponential scale lengths. In the former case, the star formation rate declines smoothly by a factor of $`5`$ over the disc’s lifetime with the consequence that a significant fraction of disc stars form early on. In the latter case, however, very few stars form initially with the rate then increasing steeply to a roughly constant value. As a consequence, most stars form late in this model. Quantitatively, for an assumed disc age of 11 Gyr, more than 50% of the final stellar mass is in place at solar neighbourhood after $`3.5`$ Gyr in the pre-assembled model whereas less than 5% is in place by same time in extreme infall model.

It is well-established that stars with a considerable range in ages exist in our local neighbourhood, a non-negligible fraction of which appear to be $`\mathrm{`}`$old’, ie. ages $`>8`$ Gyr (eg. Edvardsson et al 1993, Binney et al 2000, Rocha-Pinto et al 2000). Rocha-Pinto et al (2000) have recently re-determined the star formation history of the local Galactic disk from a sample of 552 late-type dwarfs exhibiting chromospheric activity. Stellar ages are derived using a new metallicity-dependent chromospheric activity-age relation (see their Table 3). Assuming this sample provides an unbiased representation of the total thin disc population, we use their age distribution to infer the cumulative mass in stars present in the solar neighbourhood as a function of time as shown in Figure 9. This curve shows a good agreement with the predictions of our extreme infall viscous model. While the Rocha-Pinto data suggest a slightly higher star formation rate at early times (hence a small excess of very old stars), the overall similarity of the two star formation histories at later times is quite striking. We reiterate that our star formation history is largely determined by the infall history at this radius, and not by viscous flows. Although the possibility remains that some observational biases exist within the Rocha-Pinto et al sample – for example, the onset and duration of stellar chromospheric activity is poorly understood (eg. Balnius & Vaughan 1985) – it seems unlikely that these biases could alter the derived star formation history enough to make it more closely resemble the predictions of our pre-assembled model.

Knowledge of the age distributions of stars at large radii in discs (for example, at 5 scale lengths) would have even greater leverage for model comparisons, since it is precisely these parts where the predictions of the extreme infall and pre-assembled models deviate the most. Unfortunately, essentially no information is currently available on the range of stellar ages present in these parts. Due to our location within the Milky Way disc, this issue may be more easily addressed in the future by deep studies of the resolved stellar populations in nearby external galaxies than in our own.

The effects of angular momentum transport alone mean that viscous models naturally predict inside-out growth for galactic discs whether or not infall is also present. Associating the $`\mathrm{`}`$mean age’ of the stellar population with the age by which 50% of the final stellar mass is in place, our models suggest an age gradient of $`0.7`$ Gyr over the radial range spanning 2–5 scale lengths in discs (see Figures 4 and 8). Inside-out formation has long been suspected in disc galaxies on the basis of radial colour and metallicity gradients (eg. Tinsley & Larson 1978). For example, de Jong (1996b) finds that radial colour gradients (in the sense of bluer colours at larger radii) extending over several scale lengths are common in spirals and are largely consistent with trends in stellar age alone. Using stellar population synthesis models, Bell and de Jong (2000) have explicitly quantified the radial age variation implied by these colour gradients and find the typical large spiral to exhibit an age gradient of $`0.8`$ Gyr per K-band scalelength, which also agrees with the predictions of our models.

Finally, we note that there is still significant debate as to the importance of continued gaseous infall for disc galaxy evolution, with few existing direct observational constraints. Galactic high velocity clouds may be a signature of gas accretion onto the Milky Way, but the lack of information concerning the distances of these clouds hinders a clear understanding of their origin. Blitz et al (1999) have carried out a simple simulation of the formation and evolution of the Local Group in which 106 particles (each representing a 106 $`\mathrm{M}_{}`$ gas cloud) are subject to the gravitational forces of M31, the Milky Way and the tidal field of the other Local Group members. Any cloud falling towards the Milky Way or M31 which passes within 100 comoving kpc of their centres is assumed to be accreted by that galaxy. They find present-day neutral gas accretion rates of $`1`$$`\mathrm{M}_{}`$ yr-1 for both galaxies, with the rates being somewhat higher in the distant past. The overall evolution is consistent with an exponential decline with an e-folding time of about $`5\times 10^9`$ yr. On the other hand, our extreme infall model, when normalised to the Milky Way, requires a constant infall rate of $`4.5`$$`\mathrm{M}_{}`$ yr-1. Our infall models therefore require gas accretion rates that are slightly higher but of the same order of magnitude as those found in Blitz’s simulation.

### 4.2 Constraints from Distant Galaxies

An exciting recent development has been the availability of quantitative structural parameters for moderate redshift disc galaxies which can be used to directly map the size evolution of galaxies with cosmic epoch (eg. Roche et al 1998, Lilly et al 1998, Simard et al 1999, Giallongo et al 2000). Based on a sample of discs drawn from the CFRS and LDSS redshift surveys, Lilly et al (1998) find the size function of disc scale lengths in disc-dominated galaxies (ie. bulge-to-total ratios $`0.5`$) to be roughly constant to z$`1`$, at least for large discs with R$`{}_{D}{}^{}>3.2h_{50}^1`$ kpc. In a similar study, Simard et al (1999) analyse the size-magnitude relation defined by disc-dominated field galaxies in the DEEP survey. Carefully considering possible selection effects, these authors also conclude that there is no significant evidence for any evolution in the size-magnitude distribution of disc galaxies over the redshift range $`0.1<z<1.1`$ and that a significant number of discs lie close to the canonical Freeman relation at all redshifts probed.

Taken together, these studies suggest only a very modest change in stellar exponential scale length with epoch out to $`z`$ 1 (however see Mao et al (1998) for an alternative interpretation of the same dataset). For an Einstein de Sitter Universe, the interval between a redshift of unity and the present epoch can be associated with the period from $`0.31.0`$ times the star formation interval in our models, assuming that star formation in the disc begins early on. Inspection of Figures 6 and 7 indicates that the disc scale lengths are expected to increase by about a factor of 3 over this redshift range if gaseous infall is important. Thus, in contrast to the Milky Way constraints discussed in the previous section, it would appear that observations of moderate redshift galaxies lend support for a scenario in which continued infall – as we have modelled it – is not important for building galactic discs. We note that the global star formation rates of the discs which form in our pre-assembled and extreme infall models decrease and remain roughly constant with time respectively, neither of which are inconsistent with current observations of large disc galaxies.

It should be kept in mind, though, that quantitative characterization of distant galaxy profiles is still plagued by a number of difficulties. For the highest redshift galaxies in the aforementioned samples, exponential profiles are generally fit over 1–2 disc scale lengths and hence are particularly susceptible to errors in the bulge and disc separation. Furthermore, studies to date have measured structural parameters in a single passband over the redshift range of interest. If strong colour gradients are present in these systems, this could conspire with the $`k`$-correction to mimic a constant scale length in a fixed passband over a range in redshift even though the rest-frame scale lengths are evolving. Colour gradients significantly in excess of those measured in local galaxies would be required for this to be an important effect however.

## 5 Conclusions

We have explored models for the evolution of disc galaxies that invoke simultaneous star formation and viscous redistribution of angular momentum, as well as a cosmologically-motivated prescription for gaseous infall. The existence of galaxies which exhibit exponential disc profiles over significantly more than $`5`$ scale lengths (eg. Weiner et al 2000, Barton & Thompson 1997) lend strong support for viscous models as viscosity may be the only viable mechanism for transporting material to these extreme radii, while maintaining a smooth exponential profile in the stellar disc. Following Lin & Pringle (1987a), we fix the relationship between the star formation and viscosity prescriptions such that the characteristic timescale for material to move a fractional radius of order unity is comparable with the timescale on which it is converted into stars (i.e. $`t_{}t_\nu `$). While the inclusion of viscous flows is essential for ensuring that the stellar disc is always exponential, we find that such flows play essentially no role in determining the evolution of the disc scale length, or the star formation history over the inner $`5`$ scale lengths (corresponding to the maximum infall radius in our model). Indeed, the dominant process governing the size evolution and star formation history of the luminous stellar disc is that of gaseous infall.

We summarize our main results as follows:
1. In models in which the main infall phase precedes the onset of star formation (ie. a pre-assembled gas disc), we find the exponential scale length to be rather invariant with time over the entire history of the disc. This result holds true for a variety of popular star formation prescriptions. The value of the stellar disc scale length is primarily determined by mean specific angular momentum of the initial gas disc.
2. In the pre-assembled model, the star formation rate smoothly declines with time within the inner 5 scale lengths (approximately the maximum infall radius, $`r_{init}`$). As a result, there is a considerable $`\mathrm{`}`$old’ (ie. $`>8`$ Gyr) population present at these radii, with half of the final stellar mass in place more than 6 Gyr ago. The onset of star formation is delayed at radii greater than $`r_{init}`$ as the outer disc is formed entirely by material viscously transported from smaller radii; the stellar populations in these parts are thus increasingly young with radius.
3. If infall continues during the main star-forming phase of the galaxy’s evolution, the stellar disc remains approximately exponential at all times, but now the scale length increases towards the present epoch, in simple proportion to the average specific angular momentum of the material fallen in to the galactic disc to date. For example, for an isothermal halo density profile and constant spin parameter, the scale length is directly proportional to the infallen mass which is proportional to time. Regardless of whether viscosity is at work or not, late infall of high angular momentum material is the only way to increase the stellar scale length of discs within the framework of our models. 4. In contrast to the case of a pre-assembled disc, the disc stellar populations in the infall models are predominantly young, with more than half of the stars even at the location of solar neighbourhood predicted to have ages $`<5`$ Gyr. The disc at 5 scale lengths receives direct infall of halo gas only at the very end of the accretion phase, so that prior to this, gas only reaches this location as a result of outward viscous transport. As in the pre-assembled case, all material at radii beyond the maximum infall radius has been viscously transported there from smaller radii.
5. We compare our model predictions with recent observations of disc galaxies sizes and stellar populations. A particularly good agreement is found between the solar neighbourhood star formation history predicted by our extreme infall model and the recent observational determination of this quantity by Rocha-Pinto et al (2000). On the other hand, the relative invariance of disc scale length with cosmic epoch, as derived from studies of moderate redshift galaxies, appears to lend support for our pre-assembled model. Future observations that will be important for clarifying this situation include detailed studies of resolved stellar populations at large galactocentric radii ($`>`$ 5 scale lengths) in the Milky Way and other nearby spirals in order to constrain the age and star formation history of the outer disc, and larger, deeper surveys of the structural parameters of moderate redshift disc systems in various passbands.

## 6 Acknowledgments

We are grateful to Dave Schade, Nicole Vogt, Rosie Wyse and especially Rachel Somerville for useful discussions during the course of this work.


format: text

text:
# Asymptotic tails of massive scalar fields in Schwarzschild background

## I Introduction

The late-time evolution of various fields outside a collapsing star has important implications for several aspects of black-hole physics. For example, the no-hair theorem  introduced by Wheeler in the early 1970’s, states that the external field of a black hole relaxes to a Kerr-Newman field characterized solely by the black-hole’s mass, charge and angular-momentum. Thus, it is of interest to explore the dynamical mechanism responsible for the relaxation of perturbations fields outside a black hole and to determine the decay rate of the various perturbations. In addition, the dynamical mechanism of generating perturbation fields ingoing to a black hole is of interest in relation to the problem of stability of the Cauchy horizon .

It was first demonstrated by Price , regarding scalar, gravitational and electromagnetic perturbations of Schwarzschild black hole exterior, that the fields die off at late time with an inverse power-law tail. It has been shown that for a spherical-harmonic wave mode of multiple number $`l`$, a $`t^{2l2}`$ or $`t^{2l3}`$ decay tail ($`t`$ being the Schwarzschild time coordinate) dominates at late time, depending on the initial conditions. The same tail behavior occurs also for Reissner-Nordström black holes and fields of different spin. Recently the linear perturbation analyses and the nonlinear simulations have been done by Gundlach et al.  and , respectively. (See also .) In addition, the recent treatment of tails in the spacetime of a spinning black hole has been done by Barack and Ori .

Many previous works are mainly concerned with massless fields. However, the evolution of massive scalar fields will become important, for example, if one considers Kaluza-Klein theories, in which the Fourier modes can behave like massive fields. Further, the recent development of the Kaluza-Klein idea, such as Rendall-Sundrum model  in string theory, strongly motivates us to understand the evolutional features characteristic to the field mass in detail.

The physical mechanism by which late-time tails of massive scalar fields are generated may be qualitatively different from that of massless ones. In fact, it has been pointed out that the oscillatory inverse power-law tails

$$\psi t^{l{\scriptscriptstyle \frac{3}{2}}}\mathrm{sin}(mt),$$
(1)
dominates as the intermediate late-time, if the field mass $`m`$ is small  (see also ). It is clear from Eq. (1) that massive fields decay slower than massless ones, and waves with peculiar frequency $`\omega `$ quite close to $`m`$ can contribute to massive tail, while low frequency waves can contribute to massless tail. Though the oscillatory power-law form (1) has been numerically verified at intermediate late times, $`mMmt1/(mM)^2`$, it should be noted that the intermediate tails are not the final asymptotic behaviors; Another wave pattern can dominate at very late times, when it still remains very difficult to determine numerically the exact decay rate .

In the previous paper , hereafter referred to as paper I, we have analytically found that the transition from the intermediate behavior to the asymptotic one occurs in nearly extreme Reissner-Nordström background. The oscillatory inverse power-law behavior of the dominant asymptotic tail is approximately given by

$`\psi t^{5/6}\mathrm{sin}(mt),`$ (2)
independently of angular momentum wave number $`l`$, and the decay becomes slower than the intermediate ones. The origin of such an asymptotic tail is expected to be a resonance backscattering due to curvature-induced potential. This will be supported by the relationship between the field mass and the time scale when the $`t^{5/6}`$ tail dominates; For $`mM1`$, the smaller $`mM`$ is, the later the $`t^{5/6}`$ tail begins to dominate ($`mt1/(mM)^2`$, where $`M`$ is the black-hole mass). On the other hand, for large field mass $`mM1`$, the larger $`mM`$ is, the later the $`t^{5/6}`$ tail begins to dominate ($`mtmM`$). Therefore, the time scale when the $`t^{5/6}`$ tail dominates will become minimum at $`mMO(1)`$, which means that the most effective backscattering occurs for such massive scalar fields with the Compton wave length $`1/m`$ nearly equal to the horizon radius $`M`$.

In this paper, following paper I which was devoted to the analysis for nearly extreme Reissner-Nordström background, we investigate the asymptotic tails of a massive scalar field in Schwarzschild background. The purpose of this paper is to determine analytically the decay rate and to confirm that the asymptotic tail with the decay rate of $`t^{5/6}`$ is not peculiar to the nearly extreme black hole. Nevertheless, it is sure that the form of effective potential varies according to the ratio of the black-hole charge to the mass. Then, we also study the difference of the asymptotic tail behavior in Schwarzschild background from that in nearly extreme Reissner-Nordström. In Sec. II, to analyze time evolution of massive scalar fields, we introduce the black-hole Green’s function using the spectral decomposition method . Differently from the case of nearly extreme Reissner-Nordström in paper I, it is difficult to analyze the field equations for the arbitrary field mass $`m`$. In this paper, therefore, we consider the cases of small field mass $`mM1`$ in Sec. III and large one $`mM1`$ in Sec. IV, and we compare the results with those found in nearly extreme Reissner-Nordström background. Section V is devoted to a summary.

## II Green’s function analysis

### A Massive scalar fields in Schwarzschild geometry

We consider time evolution of a massive scalar field in Schwarzschild background with the black-hole mass $`M`$. The metric is

$$ds^2=\left(1\frac{2M}{r}\right)dt^2+\left(1\frac{2M}{r}\right)^1dr^2+r^2d\mathrm{\Omega }^2,$$
(3)
and the scalar field $`\varphi `$ with the mass $`m`$ satisfies the wave equation

$$(\mathrm{}m^2)\varphi =0.$$
(4)
Resolving the field into spherical harmonics

$$\varphi =\underset{l,m}{}\frac{\psi ^l(r)}{r}Y_{lm}(\theta ,\phi ),$$
(5)
hereafter we omit the index $`l`$ of $`\psi ^l`$ for simplicity, and we obtain a wave equation for each multiple moment:

$$\psi _{,tt}\psi _{,r_{}r_{}}+V\psi =0,$$
(6)
where $`r_{}`$ is the tortoise coordinate defined by

$$dr_{}=\frac{dr}{1\frac{2M}{r}},$$
(7)
and the effective potential $`V`$ is

$$V=\left(1\frac{2M}{r}\right)\left[\frac{l(l+1)}{r^2}+\frac{2M}{r^3}+m^2\right].$$
(8)
### B The black-hole Green’s function

The time evolution of a massive scalar field is given by

$$\psi (r_{},t)=\left[G(r_{},r_{}^{};t)\psi _t(r^{},0)+G_t(r_{},r_{}^{};t)\psi (r_{}^{},0)\right]𝑑r_{}^{}$$
(9)
for $`t>0`$, where the (retarded) Green’s function $`G(r,r^{};t)`$ is defined as

$$\left[\frac{^2}{t^2}\frac{^2}{r_{}^2}+V\right]G(r_{},r_{}^{};t)=\delta (t)\delta (r_{}r_{}^{}).$$
(10)
The causality condition gives the initial condition that $`G(r,r^{};t)=0`$ for $`t0`$. In order to find $`G(r,r^{};t)`$ we use the Fourier transform

$$\stackrel{~}{G}(r_{},r_{}^{};\omega )=_0^{}^+\mathrm{}G(r_{},r_{}^{};t)e^{i\omega t}𝑑t.$$
(11)
The Fourier transform is analytic in the upper half $`\omega `$ plane, and the corresponding inversion formula is

$$G(r_{},r_{}^{};t)=\frac{1}{2\pi }_{\mathrm{}+ic}^{\mathrm{}+ic}\stackrel{~}{G}(r_{},r_{}^{};\omega )e^{i\omega t}𝑑\omega $$
(12)
where $`c`$ is some positive constant. The usual procedure is to close the contour of integration into the lower half of the complex frequency plane shown in Fig. 1. Then, the late-time tail behaviors which are our main concern should be given by the integral along the branch cut placed along the interval $`m\omega m`$.

Now the Fourier component of the Green’s function $`\stackrel{~}{G}(r_{},r_{}^{};\omega )`$ can be expressed in terms of two linearly independent solutions for the homogeneous equation

$$\left(\frac{d^2}{dr_{}^2}+\omega ^2V\right)\stackrel{~}{\psi }_i=0i=1,2.$$
(13)
The boundary condition for the basic solution $`\stackrel{~}{\psi }_1`$ is that it describes purely ingoing waves crossing the event horizon, i.e.,

$$\stackrel{~}{\psi }_1e^{i\omega r_{}},$$
(14)
at $`r_{}\mathrm{}`$, while $`\stackrel{~}{\psi }_2`$ is required to damp exponentially at spatial infinity, i.e.,

$$\stackrel{~}{\psi }_2e^{\varpi r_{}},$$
(15)
at $`r_{}\mathrm{}`$, where $`\varpi \sqrt{m^2\omega ^2}`$. Because the complex conjugate $`\stackrel{~}{\psi }_1^{}`$ is also a solution for Eq. (13), $`\stackrel{~}{\psi }_2`$ can be written by the linear superposition

$$\stackrel{~}{\psi }_2=\alpha \stackrel{~}{\psi }_1+\beta \stackrel{~}{\psi }_1^{},$$
(16)
and the Wronskian is estimated to be

$$W(\omega )=\stackrel{~}{\psi }_1\stackrel{~}{\psi }_{2,r_{}}\stackrel{~}{\psi }_{1,r_{}}\stackrel{~}{\psi }_2=2i\omega \beta .$$
(17)
Using these two solutions, the Green’s function can be written by

$$\stackrel{~}{G}(r_{},r_{}^{};\omega )=\frac{1}{2i\omega \beta }\{\begin{array}{cc}\stackrel{~}{\psi }_1(r_{}^{},\omega )\stackrel{~}{\psi }_2(r_{},\omega ),\hfill & r_{}^{}>r_{},\hfill \\ \stackrel{~}{\psi }_1(r_{},\omega )\stackrel{~}{\psi }_2(r_{}^{},\omega ),\hfill & r_{}^{}<r_{}.\hfill \end{array}$$
(18)
The contribution $`G^C`$ from the branch cut to the Green’s function is reduced to

$`G^C(r_{},r_{}^{};t)`$ $`=`$ $`{\displaystyle \frac{1}{4\pi i}}{\displaystyle _{\mathrm{cut}}}{\displaystyle \frac{1}{\omega }}{\displaystyle \frac{\alpha }{\beta }}\stackrel{~}{\psi }_1(r_{}^{},\omega )\stackrel{~}{\psi }_1(r_{},\omega )e^{i\omega t}𝑑\omega .`$ (19)
Then the main task to evaluate $`G^C`$ is to derive the coefficients $`\alpha `$ and $`\beta `$.

### C Effective contributions to tail behaviors

Here we point out the at very late times $`mt1`$ the rapidly oscillating term $`e^{i\omega t}`$ leads to a mutual cancellation between the positive and the negative parts of the integrand in Eq.(19). Then, the effective contribution to the tail behaviors arises from the narrow range

$`ϵ2M\varpi 1,`$ (20)
which was also claimed in .

## III small-mass case

Introducing the non-dimensional variable defined as

$`x`$ $``$ $`{\displaystyle \frac{r}{2M}},`$ (21)
Eq.(13) is rewritten by

$`{\displaystyle \frac{d^2\stackrel{~}{\psi }}{dx^2}}+{\displaystyle \frac{1}{x(x1)}}{\displaystyle \frac{d\stackrel{~}{\psi }}{dx}}+\left[{\displaystyle \frac{4M^2\omega ^2x^2}{(x1)^2}}{\displaystyle \frac{4M^2m^2x}{x1}}{\displaystyle \frac{l(l+1)}{x(x1)}}{\displaystyle \frac{1}{x^2(x1)}}\right]\stackrel{~}{\psi }=0.`$ (22)
To study analytically the mode solutions $`\psi `$, the asymptotic matching between the inner and outer solutions was successfully used in paper I for the potential $`V`$ in nearly extreme Reissner-Nordström background. However, to apply the same method to Eq.(22), it is necessary to assume the field mass to be very small or very large. In this section we consider the small-mass case $`mM1`$.

### A Mode solutions

#### 1 The inner region $`1x1/m^2M^2`$

For the small-mass case such as $`mM1`$, truncating the terms of the order $`m^2M^2`$, we can approximate Eq.(22) by

$`{\displaystyle \frac{d^2\stackrel{~}{\psi }}{dx^2}}+\left({\displaystyle \frac{1}{x1}}{\displaystyle \frac{1}{x}}\right){\displaystyle \frac{d\stackrel{~}{\psi }}{dx}}+\left[{\displaystyle \frac{4M^2m^2}{(x1)^2}}+{\displaystyle \frac{4M^2m^2l^2l1}{x(x1)}}+{\displaystyle \frac{1}{x^2}}\right]\stackrel{~}{\psi }`$ $`=`$ $`0.`$ (23)
Now the solution $`\psi _1`$ satisfying the boundary condition (14) can be written using the hyper-geometric function $`F`$ as follows,

$`\psi _1`$ $`=`$ $`x^{K_{++}}(x1)^{2iMm}F(K_{},K_{};14iMm;1{\displaystyle \frac{1}{x}})`$ (24)
$`=`$ $`x^{K_{++}}(x1)^{2iMm}{\displaystyle \frac{\mathrm{\Gamma }(14iMm)\mathrm{\Gamma }(2\mu )}{\mathrm{\Gamma }(K_+)^2}}F(K_{},K_{};12\mu ;{\displaystyle \frac{1}{x}})`$ (26)
$`+x^{K_+}(x1)^{2iMm}{\displaystyle \frac{\mathrm{\Gamma }(14iMm)\mathrm{\Gamma }(2\mu )}{\mathrm{\Gamma }(K_{})^2}}F(K_+,K_+;1+2\mu ;{\displaystyle \frac{1}{x}}),`$
where $`\mu `$ and $`K_{\pm \pm }`$ are

$`\mu `$ $`=`$ $`\sqrt{\left(l+{\displaystyle \frac{1}{2}}\right)^28M^2m^2}`$ (27)
and

$`K_{\pm \pm }`$ $`=`$ $`{\displaystyle \frac{1}{2}}\pm \mu \pm 2iMm,`$ (28)
respectively. and we used the linear transformation formulas (15.3.6) of  in the second equality of Eq. (24). If estimated in the region $`x1`$, we obtain

$`\stackrel{~}{\psi }_1`$ $``$ $`{\displaystyle \frac{\mathrm{\Gamma }(14iMm)\mathrm{\Gamma }(2\mu )}{\mathrm{\Gamma }(K_+)^2}}x^{\mu +\frac{1}{2}}+{\displaystyle \frac{\mathrm{\Gamma }(14iMm)\mathrm{\Gamma }(2\mu )}{\mathrm{\Gamma }(K_{})^2}}x^{\mu +\frac{1}{2}},`$ (29)
which is necessary for asymptotic matching in the overlap region $`1x1/(mM)^2`$ with the outer solutions valid in the region $`x1`$.

#### 2 The outer region $`x1`$

For large $`x`$, Eq.(22) is approximated by

$`{\displaystyle \frac{d^2\stackrel{~}{\psi }}{dx^2}}+\left[{\displaystyle \frac{4M^2m^2}{x}}ϵ^2+{\displaystyle \frac{8M^2m^2l(l+1)}{x^2}}\right]\stackrel{~}{\psi }`$ $`=`$ $`0.`$ (30)
Then we can write the solutions using Whittaker’s functions. The mode solution $`\psi _2`$ satisfying the boundary condition (15) should be

$`\psi _2`$ $`=`$ $`W_{\sigma ,\mu }(2ϵx),`$ (31)
where

$`\sigma `$ $`=`$ $`{\displaystyle \frac{2M^2m^2}{ϵ}}ϵ.`$ (32)
If estimated in the region $`x1/ϵ`$, we obtain

$`\stackrel{~}{\psi }_2`$ $``$ $`{\displaystyle \frac{\mathrm{\Gamma }(2\mu )}{\mathrm{\Gamma }(\frac{1}{2}\mu \sigma )}}(2ϵ)^{\mu +\frac{1}{2}}x^{\mu +\frac{1}{2}}+{\displaystyle \frac{\mathrm{\Gamma }(2\mu )}{\mathrm{\Gamma }(\frac{1}{2}+\mu \sigma )}}(2ϵ)^{\mu +\frac{1}{2}}x^{\mu +\frac{1}{2}}.`$ (33)
#### 3 Matching

We can match the solution (29) valid at $`1x1/m^2M^2`$ with the solution (33) valid at $`1x1/ϵ`$ in the overlap region and determine the coefficients $`\alpha `$ and $`\beta `$ in Eq.(19). Then the ratio is given as follows,

$`{\displaystyle \frac{\alpha }{\beta }}(|\omega |,ϵ)|_{|\omega |m}`$ $`=`$ $`[(2ϵ)^{\mu +\frac{1}{2}}{\displaystyle \frac{\mathrm{\Gamma }(2\mu )^2\mathrm{\Gamma }(1+4iMm)}{\mathrm{\Gamma }(\frac{1}{2}\mu \sigma )\mathrm{\Gamma }(K_+)^2}}`$ (35)
$`+(2ϵ)^{\mu +\frac{1}{2}}{\displaystyle \frac{\mathrm{\Gamma }(2\mu )^2\mathrm{\Gamma }(1+4iMm)}{\mathrm{\Gamma }(\frac{1}{2}+\mu \sigma )\mathrm{\Gamma }(K_{++})^2}}]\times (\mathrm{complex}\mathrm{conjugate})^1.`$
### B Intermediate tails

The effective contribution to the integral in Eq. (19) is claimed to be limited to the range $`|\omega m|=O(1/t)`$ or equivalently $`\varpi =O(\sqrt{m/t})`$, since the rapidly oscillating term $`e^{i\omega t}`$ which leads to a mutual cancellation between the positive and the negative parts of the integrand (see ). The intermediate tails becomes dominant at the late times in the range

$`mMmt{\displaystyle \frac{1}{(mM)^2}},`$ (36)
when the integral (19) should be estimated under the condition

$`\sigma {\displaystyle \frac{2m^2M^2}{ϵ}}=O(mM\sqrt{mt})1.`$ (37)
As was discussed in paper I, the small value of $`\sigma `$ represents that the backscattering due to the spacetime curvature is not effective at intermediate late times. Then, using the condition (37), we obtain

$`{\displaystyle \frac{\alpha }{\beta }}(|\omega |,ϵ){\displaystyle \frac{\alpha }{\beta }}(|\omega |,e^{i\pi }ϵ)|_{|\omega |m}i{\displaystyle \frac{2^{2l+1}l!^6}{(2l)!^2(2l+1)!^4}}mMϵ^{2l+1},`$ (38)
which is identical with Eq. (51) in paper I. Therefore, it is obvious that the same intermediate tails as Eq. (1) dominate at intermediate late times (36), which was numerically supported by .

### C Asymptotic tails

As discussed in paper I, the intermediate tails can not be an asymptotic behavior, and the long-term evolution from the intermediate behavior to the final one should occur. The asymptotic tail becomes dominant at very late times such that

$`mt{\displaystyle \frac{1}{m^2M^2}},`$ (39)
when the effective contribution to the integral (19) arises from the region

$`\sigma {\displaystyle \frac{2m^2M^2}{ϵ}}1,`$ (40)
which means backscattering effect due to curvature-induced potential dominates. Then, the ratio of $`\alpha `$ to $`\beta `$ is approximately given by

$`{\displaystyle \frac{\alpha }{\beta }}(|\omega |,ϵ)|_{|\omega |m}`$ $``$ $`{\displaystyle \frac{\eta ^{}e^{i\pi \sigma }+\gamma ^{}e^{i\pi \sigma }}{\eta e^{i\pi \sigma }+\gamma e^{i\pi \sigma }}},`$ (41)
where

$`\eta `$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }(2\mu )^2\mathrm{\Gamma }(14iMm)}{\mathrm{\Gamma }(K_{})^2}}(4m^2M^2)^\mu e^{i\pi \mu }+{\displaystyle \frac{\mathrm{\Gamma }(2\mu )^2\mathrm{\Gamma }(14iMm)}{\mathrm{\Gamma }(K_+)^2}}(4m^2M^2)^\mu e^{i\pi \mu },`$ (42)
and

$`\gamma `$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }(2\mu )^2\mathrm{\Gamma }(14iMm)}{\mathrm{\Gamma }(K_{})^2}}(4m^2M^2)^\mu e^{i\pi \mu }+{\displaystyle \frac{\mathrm{\Gamma }(2\mu )^2\mathrm{\Gamma }(14iMm)}{\mathrm{\Gamma }(K_+)^2}}(4m^2M^2)^\mu e^{i\pi \mu }.`$ (43)
As was shown in the paper I, the contribution from the Green’s function to the asymptotic tail part corresponds to the integral along the dashed line in Fig. 1 which is approximated by

$`G^C(r_{},r_{}^{};t)`$ $``$ $`{\displaystyle \frac{1}{4\pi mi}}\stackrel{~}{\psi }_1(r_{},m)\stackrel{~}{\psi }_1(r_{}^{},m){\displaystyle _{\mathrm{dashed}\mathrm{line}}}e^{i(2\pi \sigma \omega t)}e^{i\phi _s}𝑑\omega +(\mathrm{complex}\mathrm{conjugate}),`$ (44)
where the phase $`\phi _s`$ is defined by

$$e^{i\phi _s}=\frac{\eta ^{}+\gamma ^{}e^{2i\pi \sigma }}{\eta +\gamma e^{2i\pi \sigma }},$$
(45)
and it remains in the range $`0\phi _s2\pi `$, even if $`\sigma `$ becomes very large, since we have

$`|\eta |^2|\gamma |^2`$ $`=`$ $`{\displaystyle \frac{2\pi Mm}{\mu }}>\mathrm{\hspace{0.25em}0}.`$ (46)
Because the terms $`e^{i\omega t}`$ and $`e^{2i\pi \sigma }`$ in Eq.(44) are rapidly oscillating at very late times, the saddle-point integration allows us to evaluate accurately the asymptotic behaviors. Let us introduce the variable $`a`$ defined by

$`a\left(1{\displaystyle \frac{\omega }{m}}\right)^{1/2}(mt)^{1/3}.`$ (47)
Then, the oscillation terms $`e^{i(2\pi \sigma \omega t)}`$ in Eq.(44) can be approximately rewritten into the form

$`e^{i(mt)^{1/3}f_s(a)}e^{imt}`$ (48)
in the limits $`mt\mathrm{}`$ and $`\omega /m1`$, by keeping $`a`$ to be finite. Here we have

$`f_s(a)={\displaystyle \frac{\sqrt{2}\pi mM}{a}}+a^2.`$ (49)
The saddle point is found to exist at

$`a=a_0\left({\displaystyle \frac{\pi mM}{\sqrt{2}}}\right)^{1/3}.`$ (50)
At the saddle point the parameter $`ϵ`$ defined by Eq.(20) is given by

$`ϵ_02M\sqrt{m^2\omega _0^2}2\left({\displaystyle \frac{2\pi }{mt}}\right)^{1/3}(mM)^{4/3},`$ (51)
or equivalently

$`\sigma _0mM\left({\displaystyle \frac{2\pi mM}{mt}}\right)^{1/3}`$ (52)
which is equal to the value obtained in the nearly extreme case in paper I. Let us give the Taylor expansion of $`f`$ near the point $`a=a_0`$ as follows,

$`f_s(a)`$ $``$ $`f_s(a_0)+{\displaystyle \frac{1}{2}}(aa_0)^2f_s^{\prime \prime }(a_0)`$ (53)
$`=`$ $`3\left({\displaystyle \frac{\pi mM}{\sqrt{2}}}\right)^{2/3}+3(aa_0)^2.`$ (54)
Then, the integral in Eq.(44) can be approximately estimated to be

$`G^C(r_{},r_{}^{};t)`$ $``$ $`{\displaystyle \frac{1}{2\sqrt{3}}}(2\pi )^{\frac{5}{6}}(mM)^{\frac{1}{3}}(mt)^{\frac{5}{6}}\mathrm{sin}(mt3(2\pi mM)^{2/3}(mt)^{1/3}/2\phi _s(\sigma _0)+3\pi /4)`$ (56)
$`\times \stackrel{~}{\psi }_1(r_{},m)\stackrel{~}{\psi }_1(r_{}^{},m),`$
at very late times $`mt1/(mM)^2`$. The oscillation has the period $`2\pi /m`$, and is modulated by the two types of long-term phase shifts. First, the term $`(3/2)(2\pi mM)^{2/3}(mt)^{1/3}`$ represents a monotonously increasing phase shift. Second, the term $`\phi _s(\sigma _0)`$ represents a periodic phase shift with the amplitude $`|\phi _s|`$ smaller than $`2\pi `$, which will be studied in the following. Eq.(56) is similar to Eq.(69) in paper I, though the corresponding terms of the phase shift are omitted in Eq.(69) in paper I since we have focused our attention to the decay rate.

### D Phase shift

Now we study the periodic phase-shift effect caused by the term $`\phi _s(\omega _0)`$. From the equation

$`{\displaystyle \frac{d\phi _s}{d\sigma }}`$ $`=`$ $`0,`$ (57)
we find the maximum and minimum extremes, denoted by $`\phi _{s+}(\sigma _+)`$ and $`\phi _s(\sigma _{})`$ of $`\phi _s`$ at $`\sigma =\sigma _+`$ and $`\sigma =\sigma _{}`$, which given by

$`\phi _{s\pm }(\sigma _\pm )`$ $`=`$ $`\pm 2\theta +\phi _{s0},`$ (58)
and

$`2\pi \sigma _\pm `$ $`=`$ $`\pi {\displaystyle \frac{\pi }{2}}\pm \theta +const,`$ (59)
where

$`\mathrm{sin}\theta `$ $`=`$ $`\left|{\displaystyle \frac{\gamma }{\eta }}\right|`$ (60)
and

$`\mathrm{exp}(i\phi _{s0}){\displaystyle \frac{\eta ^{}}{\eta }}.`$ (61)
Expanding $`|\gamma /\eta |`$ with respect to $`mM1`$, we obtain

$`\left|{\displaystyle \frac{\gamma }{\eta }}\right|`$ $`=`$ $`1\delta ,`$ (62)
and

$`\delta `$ $``$ $`\pi \left({\displaystyle \frac{64}{15}}\right)^2(mM)^3(\mathrm{for}l=0),`$ (63)
$``$ $`{\displaystyle \frac{2^{4l+4}\pi (l!)^4}{(2l+1)(2l)!^4}}(mM)^{4l+3}(\mathrm{for}l0)`$ (64)
in Schwarzschild background. The phase-shift term $`\delta `$ is characteristic to the black hole. For example, for nearly extreme Reissner-Nordström background , the value of $`\delta `$ is given by

$`\delta _{\mathrm{nearly}\mathrm{extreme}}`$ $``$ $`\pi \left({\displaystyle \frac{50}{23}}\right)^2(mM)^3(\mathrm{for}l=0),`$ (65)
$``$ $`{\displaystyle \frac{2^{4l+3}\pi (l!)^2}{(2l+1)^2(2l)!^4}}(mM)^{6l+3}(\mathrm{for}l0).`$ (66)
From Eqs. (58) and (59) the oscillation amplitude of $`\phi _s`$ is given by

$`\mathrm{\Delta }\phi _s`$ $``$ $`\phi _{s+}\phi _s=4\theta `$ (67)
$``$ $`2\pi \sqrt{2\delta },`$ (68)
and the time interval between the two extremes is estimated to be

$`2\pi \mathrm{\Delta }\sigma `$ $``$ $`2\pi (\sigma _{}\sigma _+)=\pi 2\theta `$ (69)
$``$ $`\sqrt{8\delta }.`$ (70)
These equations mean that the smaller the field mass is, the more rapidly the change of the phase $`\phi _s`$ from the maxixum to the minimum occurs (see Fig. 2). The phase shift oscillates at the time scale of $`1/(m^3M^2)`$, which is much longer than the period of $`2\pi /m`$. So it is a long-term effect which modulates the basic oscillation.

## IV large mass case

### A Mode solutions

In this section we consider large-mass case such as $`mM1`$. It is enough to consider the narrow region (20) in the same manner with the small-mass case. Introducing the function $`\xi `$ defined as

$`\stackrel{~}{\psi }=\sqrt{{\displaystyle \frac{r}{r2M}}}\xi =\sqrt{{\displaystyle \frac{x}{x1}}}\xi ,`$ (71)
the mode equation (22) can be approximately given by

$`{\displaystyle \frac{d^2\xi ^2}{dx^2}}+\left[{\displaystyle \frac{4M^2\omega ^2x^2}{(x1)^2}}{\displaystyle \frac{4M^2m^2x}{x1}}\right]\xi `$ $`=`$ $`0,`$ (72)
since the other terms become negligibly small all over the region $`x1`$ in comparison with the two terms in the potential. Then we can write solutions for Eq.(72) using Whittaker’s functions. The solutions $`\psi _1`$ and $`\psi _2`$ satisfying the boundary condition (14) and (15), respectively, are

$`\stackrel{~}{\psi }_1`$ $`=`$ $`\sqrt{{\displaystyle \frac{x}{x1}}}M_{\sigma ,\rho }\left(2ϵ(x1)\right),`$ (73)
and

$`\stackrel{~}{\psi _2}`$ $`=`$ $`\sqrt{{\displaystyle \frac{x}{x1}}}W_{\sigma ,\rho }\left(2ϵ(x1)\right),`$ (74)
where

$`\rho `$ $`=`$ $`i\sqrt{4M^2\omega ^2{\displaystyle \frac{1}{4}}}.`$ (75)
From Eqs. (73) and (74), $`\psi _2`$ can be rewritten as

$$\stackrel{~}{\psi }_2=\frac{\mathrm{\Gamma }(2\rho )}{\mathrm{\Gamma }(\frac{1}{2}\rho \sigma )}\psi _1+\frac{\mathrm{\Gamma }(2\rho )}{\mathrm{\Gamma }(\frac{1}{2}+\rho \sigma )}\psi _1^{}$$
(76)
and the ratio of $`\alpha `$ to $`\beta `$ in Eq.(16) is

$`{\displaystyle \frac{\alpha }{\beta }}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }(2\rho )\mathrm{\Gamma }(\frac{1}{2}+\rho \sigma )}{\mathrm{\Gamma }(2\rho )\mathrm{\Gamma }(\frac{1}{2}\rho \sigma )}}.`$ (77)
### B Asymptotic tail

As was discussed in paper I, no intermediate tails appear for the scalar field with large mass, namely $`mM1`$. At late times given by Eq.(40), the ratio of $`\alpha `$ to $`\beta `$ is approximately written by

$`{\displaystyle \frac{\alpha }{\beta }}(|\omega |,ϵ)|_{|\omega |m}`$ $``$ $`{\displaystyle \frac{\zeta ^{}e^{i\pi \sigma }\chi ^{}e^{i\pi \sigma }}{\zeta e^{i\pi \sigma }\chi e^{i\pi \sigma }}}ϵ^{2\rho },`$ (78)
where

$`\zeta `$ $`=`$ $`\mathrm{\Gamma }(2\rho )(2m^2M^2)^\rho e^{i\pi \rho },`$ (79)
and

$`\chi `$ $`=`$ $`\mathrm{\Gamma }(2\rho )(2m^2M^2)^\rho e^{i\pi \rho }.`$ (80)
Substituting Eq.(78) into Eq.(19), we obtain the contribution from Green’s function to the asymptotic tail;

$`G^C(r_{},r_{}^{};t)`$ $``$ $`{\displaystyle \frac{1}{4\pi mi}}\stackrel{~}{\psi }_1(r_{},m)\stackrel{~}{\psi }_1(r_{}^{},m){\displaystyle _{\mathrm{dashed}\mathrm{line}}}e^{2i\pi \sigma 2\rho \mathrm{ln}ϵi\omega t}e^{i\phi _l}𝑑\omega +(\mathrm{complex}\mathrm{conjugate}),`$ (81)
where the phase $`\phi _l`$ is defined by

$$e^{i\phi _l}=\frac{\zeta ^{}+\chi ^{}e^{2i\pi \sigma }}{\zeta +\chi e^{2i\pi \sigma }},$$
(82)
and it remains in the range $`0\phi _l2\pi `$, even if $`\sigma `$ becomes very large, since we have

$`|\zeta |^2|\chi |^2`$ $`=`$ $`\mathrm{\Gamma }(2\rho )\mathrm{\Gamma }(2\rho )\left(e^{2i\pi \rho }e^{2i\pi \rho }\right)>\mathrm{\hspace{0.25em}0}.`$ (83)
Because the terms $`e^{i\omega t}`$, $`e^{2i\pi \sigma }`$ and $`e^{2\rho \mathrm{ln}ϵ}`$ in Eq.(81) are rapidly oscillating at very late times, the saddle-point integration allows us to evaluate accurately the asymptotic behaviors. Just like the small-mass case in the previous section, using the parameter $`a`$ defined as Eq. (47), the oscillation terms $`e^{i(2\pi \sigma \omega t2\rho \mathrm{ln}ϵ)}`$ in Eq.(44) can be rewritten into the form

$`e^{i(mt)^{1/3}f_l(a)}e^{imt}`$ (84)
in the limits $`mt\mathrm{}`$ and $`\omega /m1`$, by keeping $`a`$ to be finite. Here we have

$`f_l(a){\displaystyle \frac{\sqrt{2}\pi mM}{a}}+a^2+{\displaystyle \frac{4mM}{(mt)^{1/3}}}\mathrm{ln}\left[{\displaystyle \frac{2\sqrt{2}mMa}{(mt)^{1/3}}}\right].`$ (85)
The saddle point can be approximated by

$`a_0`$ $``$ $`\left({\displaystyle \frac{\pi mM}{\sqrt{2}}}\right)^{1/3}`$ (86)
for asymptotically late times

$`mtmM.`$ (87)
Note that Eq. (86) is identical with Eq. (50). Evaluating Eq. (81) using the saddle point integration, we obtain the asymptotic tail behavior

$`G^C(r_{},r_{}^{};t)`$ $``$ $`{\displaystyle \frac{1}{2\sqrt{3}}}(2\pi )^{\frac{5}{6}}(mM)^{\frac{1}{3}}(mt)^{\frac{5}{6}}`$ (90)
$`\times \mathrm{sin}(mt3/2(2\pi mM)^{2/3}(mt)^{1/3}4mM\mathrm{ln}ϵ_0\phi _l(\sigma _0)+3\pi /4)`$
$`\times \stackrel{~}{\psi }_1(r_{},m)\stackrel{~}{\psi }_1(r_{}^{},m).`$
The additional relation (87) is necessary for the asymptotic tail (90) to dominate. The oscillation has the period $`m/2\pi `$ and is modulated by the two types of phase shifts. First, the terms $`(3/2)(2\pi mM)^{2/3}(mt)^{1/3}+4mM\mathrm{ln}ϵ_0`$ represent a monotonous phase shift. Second, the term $`\phi _l(\sigma _0)`$ represents a periodic phase shift since it remains within $`0\phi _l2\pi `$. The decay rate of the asymptotic tail in Schwarzschild background is $`t^{5/6}`$, which remains identical with that in nearly extreme Reissner-Nordström background shown in the paper I.

### C Phase shift

In the same manner with the small mass case, we study the periodic phase-shift effect which is caused by the term $`e^{i\phi _l}`$ in Eq.(82). In the large mass limit $`mM1`$, the ratio of $`\chi `$ and $`\zeta `$ is approximated by

$`\left|{\displaystyle \frac{\chi }{\zeta }}\right|`$ $``$ $`e^{\pi 4mM}1`$ (91)
in the Schwarzschild case. On the other hand, in the nearly extreme Reissner-Nordström limit, we obtain the large mass limit from Eqs.(59) and (60) in paper I as follows,

$`\left|{\displaystyle \frac{\chi }{\zeta }}\right|_{\mathrm{nearly}\mathrm{extreme}}`$ $``$ $`e^{\pi \sqrt{5}mM}1.`$ (92)
The oscillation amplitude $`\mathrm{\Delta }\phi _l`$ of $`\phi _l`$ is

$`\mathrm{\Delta }\phi _l`$ $``$ $`4\left|{\displaystyle \frac{\chi }{\zeta }}\right|,`$ (93)
for which Eqs. (91) and (92) give different values, though it remains same that $`\mathrm{\Delta }\phi `$ vanishes in the large-mass limit.

## V Summary

In this paper we have investigated asymptotic tail behaviors of massive scalar fields in Schwarzschild background. We have shown that the asymptotic oscillatory tail with the decay rate of $`t^{5/6}`$ and with the period $`2\pi /m`$ dominates, which is identical with that in nearly extreme Reissner-Nordström background. In addition, we have noted that the relation between the value of $`mM`$ and the time-scale when the $`t^{5/6}`$ tail begins to dominate also holds; The smaller the value of $`mM`$ is, the later the $`t^{5/6}`$ tail begins to dominate ($`mt1/(mM)^2`$). The larger the value of $`mM`$ is, the later the $`t^{5/6}`$ tail begins to dominate ($`mtmM`$).

As far as the intermediate late-time behavior is concerned, our result agrees with . However they claimed ”SI perturbation fields decay at late times slower than any power law” in , which disagree with the expressions (56) and (81). We believe that their simulation was not carried out until enough late times for the asymptotic tail to dominate or there occured numerical errors in their simulation. Also, our result is consistent with , which was done numerically (both for a test field and for the fully nonlinear cases).

In terms of the phase shift we can understand more clearly the tail behavior to be a backscattering effect due to space-time curvature. In the small-mass case $`mM1`$, the intermediate tail with an oscillatory power-law dominates at the time scale $`1mt1/(mM)^2`$ before the asymptotic tail discussed here dominates  (see also ). At both the intermediate and asymptotic time scales the frequency of a wave which contributes to tails is very close to $`m`$. However, the oscillation of the asymptotic tail is modulated by phase shift effects which do not appear in the intermediate tail. There exist two types of the phase shift, which are monotonous and periodic with time. The periodic phase shift crucially depends on the field mass. The shift angle becomes nearly equal to $`2\pi `$ in small mass limit $`mM0`$, and the shift oscillates at the time scale of $`1/(m^3M^2)`$, which is a long-term effect if compared with the basic oscillation with the period of $`2\pi /m`$. In large mass limit $`mM\mathrm{}`$, the shift angle become very small. In scattering theory, it is well-known that the phase shift of wave is connected with the presence of the scatterer. In this case, of course, the scatterer is the effective potential of black-hole spacetime. Therefore it is clear that tails of massive scalar fields are generated from backscattering due to space-time curvature. The shift angle in Schwarzschild geometry is numerically different from that in nearly extreme Reissner-Nordström geometry, but the qualitative behavior depending on $`mM`$ is similar.

We can conclude that the asymptotic tail with the qualitatively same behavior dominates both in Schwarzschild and in nearly extreme Reissner-Nordström background. Therefore it is conjectured that the oscillatory $`t^{5/6}`$ tail caused by resonance backscattering at asymptotic late times is a general feature in arbitrary Reissner-Nordström background. It remains in a future work to investigate what kind of instability of Cauchy horizons is caused by the massive tails.

Finally we give a comment in comparison with the effects induced by the field mass $`m`$ on vacuum polarization of quantum massive scalar fields in the thermal state. It has been found that the amplitude of vacuum polarization is enhanced around $`mMO(1)`$ in the case of nearly extreme Reissner-Nordström geometry, while it is not seen in Schwarzschild case . The field-mass induced effect on vacuum polarization becomes clear only in the nearly extreme limit. The reason will be explained as follows; Vacuum polarization of massive scalar fields is due to both thermal excitation induced by black-hole temperature and mass induced excitation. In nearly extreme case, of which the black-hole temperature is nearly zero, the mass induced effect becomes significant because the thermal excitation is suppressed. On the other hand, in Schwarzschild case the mass induced effect is just hidden because the thermal effect by the black-hole temperature can dominate.

We expect that resonance behavior due to the mass of a field interacting with a black hole may exist in various processes as a generic feature of black-hole geometry. In order to confirm this conjecture, it is necessary to investigate various processes of massive fields in more general black hole models. This will be one of the interesting subjects in black-hole physics.

###### Acknowledgements.

The authors thank L. M. Burko and K. Nakamura for useful comments.


text:
# Low temperature kinetics of 2D exciton gas cooling in quantum well bilayer

We study the kinetics of 2D Bose gas cooling provided Bose particles interact with 3D phonons. At low temperatures phonon emission is prohibited by the energy and the momentum conservation. We show that both particle-particle scattering and impurity scattering assist the Bose gas cooling. The temporal relaxation of temperature follows the law $`T1/\sqrt{t}`$ above the Berezinski-Kosterlitz-Thouless phase transition point and $`T1/t`$ after a Bose-Einstein 2D quasi-condensate develops.

Exciton gas in GaAs bilayer quantum well represents a system where a 2D Bose-Einstein quasi-condensation (BEqC) is possible at low temperatures. Experimental efforts  have been directed to produce such an exciton gas and to cool it down to a BEqC temperature. A short laser photo-illumination pulse excites electron and hole pairs. The so-called indirect exciton technique is used where a perpendicular electric field drags electrons and holes apart into two spatially separated layers. Then an electron and a hole bind themselves into an indirect exciton particle. This experimental setup suppresses the electron-hole recombination giving rise to a relatively long exciton life time. During the initial photoexcitation pulse newly born excitons are hot and form a non-equilibrium state. In short time after the pulse ends the exciton gas reaches an equilibrium at some effective exciton temperature which is much higher than the lattice temperature of the cold GaAs crystal. Frequent exciton-exciton collisions ensures the exciton temperature to be uniform across the bilayer. Exciton gas then start to cool down slowly due to an emission of phonons into a crystal away from the bilayer. This is the longest phase of the experiment limited only by a decay time of excitons due to the electron-hole recombination. In order to reach a BEqC point one needs both low temperature and a high density of excitons. Hence a fast cooling is essential.

Recent calculation of energy losses in 2D ideal exciton system has predicted an extremely slow cooling at low temperatures with the temporal law: $`T(t)1/\mathrm{log}(t)`$, where $`t`$ is the time . This fact is intimately related to the energy and momentum conservation which prohibits an emission of phonons by an exciton moving slower than the velocity of sound in GaAs crystal: $`c`$. Thus, the exciton gas cooling in effect stops when the exciton temperature falls below a characteristic blocking temperature: $`T_b=mc^2/2`$, where $`m`$ is the mass of exciton, even if the crystal temperature is zero.

This kinetic bottleneck problem becomes especially acute when the exciton gas is subjected to a strong perpendicular magnetic field that quenches the motion of exciton to the lowest Landau level and, thus, is helping to bind electrons and holes into exciton pairs. In this case the effective mass of an exciton is determined by the Coulomb interaction and can be much larger than either the electron or hole mass . This results in a higher blocking temperature $`T_b`$ and makes it difficult to reach low temperatures in the end.

In this communication we supplement the analysis of exciton cooling of Ref. by an addition of exciton-exciton collisions and scattering on impurities. Both events substantially assist the phonon emission. We specialize to the case of exactly zero lattice temperature which allows us to neglect exciton-phonon scattering. Actually we are dealing with a general problem of 2D Bose gas cooling provided its particles interact with 3D phonons. The universal nature of 2D scattering at low energy of incoming particles makes these two assistance mechanisms to be robust to specific details of a particle-particle or impurity potential. The latter is only assumed to be short-ranged with the characteristic interaction radius $`r_0`$ being shorter than the DeBroigle wavelength. For example, an indirect exciton interacts with an impurity or another exciton via the electron-hole dipole moment: $`e\stackrel{}{d}`$, directed along the normal to the bilayer. Hence, the radius of such a dipole interaction equals to the spacing between the electron and the hole layers: $`r_0d`$.

The Hamiltonian of particle phonon interaction can be written generally as:

$$H_{xph}=\psi ^+(\stackrel{}{r})\psi (\stackrel{}{r})\delta (z)\mathrm{\Gamma }_i(\stackrel{}{\rho }\stackrel{}{\rho ^{}})u_i(\stackrel{}{\rho ^{}})d^3\stackrel{}{\rho ^{}}d^2\stackrel{}{r}𝑑z,$$
(1)
where $`\stackrel{}{\rho }=(\stackrel{}{r},z)`$, $`\psi ^+`$ and $`\psi `$ are the particle creation and annihilation operators and $`\stackrel{}{u}`$ is a crystal deformation induced by an acoustic phonon. Wavelengths of relevant phonons are much larger than the width of the bilayer: $`d`$. The lattice deformation $`\stackrel{}{u}`$ can be expanded into the normal phonon modes as:

$$u_i(\stackrel{}{\rho })=\underset{s,\stackrel{}{q}}{}\left(\frac{\mathrm{}}{2\rho \omega _s(\stackrel{}{q})}\right)^{1/2}(e_i^sb_s^+(\stackrel{}{q})+e_i^sb_s(\stackrel{}{q}))e^{i\stackrel{}{q}\stackrel{}{\rho }},$$
(2)
where $`b_s^+`$ and $`b_s`$ are the phonon creation and annihilation operators of polarization $`s`$, $`\rho `$ is the mass density of solid, $`\omega _s(\stackrel{}{q})=cq`$ is the phonon frequency dispersion, which we assume to be isotropic and independent of phonon polarization $`s`$.

The exciton phonon interaction in GaAs crystal can be separated into piezoelectric and deformation potential parts. Lattice deformation in a piezoelectric crystal induces a polarization density: $`P_i=\beta _{ijk}_ju_k`$ , where $`\beta _{ijk}`$ is the piezoelectric tensor. This polarization interacts with the exciton dipole moment. In the limit $`qd1`$, the deformation potential for an exciton: $`\mathrm{\Theta }`$, is a sum of the deformation potentials for an electron and a hole taken at the same spatial point. The latter represents a change of the semiconductor gap due to the local compression caused by a phonon deformation. Combining the piezoelectric and deformation parts and expanding the crystal lattice deformation in acoustic phonon modes we write the the exciton-phonon vertex in the Hamiltonian (1) as :

$$\mathrm{\Gamma }_i(\stackrel{}{q})=\left(\mathrm{\Theta }q_i+ed\beta _{ijk}\frac{4\pi q_zq_jq_k}{\stackrel{}{q}^2}\right).$$
(3)
For a cubic GaAs crystal without the inversion center $`\beta _{ijk}=\beta `$ if all $`i,j,k`$ are different and zero otherwise. In the limit of large $`d`$ the piezoelectric part dominates over the deformation potential but in the experiments  $`d50A`$ and both exciton-phonon interaction terms are of the same order of magnitude: $`\mathrm{\Theta }4\pi ed\beta 10eV`$.

Amplitude of phonon emission is given to the lowest order of the perturbation theory by the following matrix element:

$$M_{if}^s(\stackrel{}{q})=<f\stackrel{}{q}s|H_{ph}|i0>$$
(4)
between initial state of Bose gas $`|i0>`$ with no phonons and the final state of Bose gas $`|f\stackrel{}{q}s>`$ with just one phonon specified by the momentum: $`\stackrel{}{q}`$ and the polarization $`s`$. We assume that the thermalization of the Bose gas due to particle-particle scattering is much faster than the slow cooling due to phonon emission. Thus, at any given time $`t`$ the Bose gas is characterized by an effective temperature: $`T(t)`$ . This temperature defines the total gas energy: $`E=E(T)`$. The Fermi Golden Rule gives the probability of phonon emission per unit time and one needs to multiply it by the phonon energy: $`\omega _s(\stackrel{}{q})=cq`$, to find the total energy losses:

$$\frac{dE}{dt}=\frac{2\pi }{\mathrm{}}\underset{f\stackrel{}{q}s}{}cq|M_{if}^s(\stackrel{}{q})|^2\delta (E_iE_fcq).$$
(5)
Eq.(5) has to be averaged over Gibbs distribution of the initial states with the effective temperature $`T(t)`$. Both the initial and the final states of the Bose gas are calculated in the interaction representation (see e.g.). Particles are confined to the 2D layer and the energy losses are proportional to the area of this layer.

In the experiments  only a relatively low density of exciton gas has been achieved. Popov has shown  that for a 2D dilute Bose gas there is a 2D Berezinski-Kosterlitz-Thouless phase transition point:

$$T_c=\frac{2\pi n\mathrm{}^2}{gm\mathrm{log}L}$$
(6)
that separates high-T almost ideal Bose gas phase from the low-T superfluid phase. Actually, Popov theory is controlled by the large logarithm:

$$L\mathrm{log}(nr_0^2)\mathrm{log}\frac{E_0}{T_c}$$
(7)
where $`n`$ is the particle density, $`E_0=\mathrm{}^2/r_0^2m`$, $`g`$ is the particle internal degeneracy . For a Bose particle: $`g=2S+1`$, where $`S`$ is the spin of particle. It was shown in Ref. that electron and hole spins flip rapidly due to the spin-orbit interaction. Thus, $`g=4`$ for a GaAs exciton.

For 2D dilute non-ideal Bose gas one can distinguish three temperature regions. At high temperatures: $`TT_c\mathrm{log}L`$, the ideal Bose gas is a good approximation. At intermediate temperatures: $`gT_c/LTT_c/\mathrm{log}L`$, an overwhelming amount of particles constitute a 2D BEqC with the density:

$$n_s=n\left(1\frac{T}{T_c}\right),$$
(8)
whereas a small fraction of thermal particles have a bare dispersion and Bose distribution with the chemical potential: $`\mu gT_c/L`$ . At low temperatures: $`T\mu `$, a weak particle-particle interaction is crucial and the Bose gas becomes a Bose-liquid with quasiparticle excitations having Bogolubov sound like dispersion. The transfer of momentum to an impurity becomes here inefficient because sound-like quasiparticles ignore point-like impurities. Hot particles with a bare dispersion are rare due to an exponentially small Boltzmann factor. In the case of excitons in GaAs crystal the intermediate-T region hardly exists at all.

Fig.1 shows two scattering possibilities for a particle accompanied by an acoustic phonon emission. The left diagram shows a scattering on the second particle and the right diagram shows a scattering on impurity. First, we treat this problem in the high-T and intermediate-T regions where particles have the bare dispersion: $`ϵ=\stackrel{}{p}^2/2m`$. 2D scattering amplitudes in both cases are isotropic and depend only on the total kinetic energy in the center of mass frame: $``$ (in the impurity case $`=ϵ`$) in the limit $`0`$:

$$F()=\frac{2\pi \mathrm{}^2}{m^{}}\left(\mathrm{log}\frac{\mathrm{}^2}{r_0^2m}\right)^1,$$
(9)
where $`m^{}=m/2`$ for particle-particle scattering and $`m^{}=m`$ for the impurity scattering (see e.g. ). Both particle-particle interaction line and the impurity line on Fig.1 correspond to the scattering amplitude $`F()`$. The total amplitude of assisted phonon emission is universal in both cases and is given by the following matrix element:

$$M_{if}^s(\stackrel{}{q})=C\left(F()\frac{1}{cq}+F(cq)\frac{1}{cq}\right)e_i^s\mathrm{\Gamma }_i(\stackrel{}{q})\sqrt{\frac{q}{2\rho c}},$$
(10)
where $`C=2`$ for the particle-particle scattering and $`C=1`$ for the impurity scattering. We neglect the phonon momentum $`\stackrel{}{q}_{||}`$ transfer to the particle because in the low temperature limit: $`q_{||}p,p^{}cm`$.

Plugging Eqs.(10,3) into Eq.(5) and taking the integral over the final state of the Bose gas we obtain the total energy losses per unit time. In the high-T region we get:

$$\frac{dE}{dt}=\frac{2\pi c^3}{\mathrm{}T_{xph}^2}\left(\frac{1}{\mathrm{log}(E_0/ϵ)}\frac{1}{\mathrm{log}(E_0/ϵ^{})}\right)^2K(ϵ,T)\delta (ϵϵ^{}cq)𝑑ϵ𝑑ϵ^{}\frac{d^3\stackrel{}{q}}{(2\pi )^3},$$
(11)
where $`K(ϵ,T)=2AnN(ϵ)(1+1/g)`$ in the case of particle-particle scattering and $`K(ϵ,T)=An_{imp}N(ϵ)/2`$ in the case of impurity scattering. $`A`$ is the total area of bilayer, $`n_{imp}`$ is the areal density of impurities and

$$N(ϵ)=\frac{1}{\mathrm{exp}\left((ϵ\mu )/T\right)1}$$
(12)
is the Bose-Einstein occupation number. It is also convenient to define a characteristic exciton-phonon temperature in the case of GaAs:

$$T_{xph}=\sqrt{\frac{\rho c^5\mathrm{}^3}{\mathrm{\Theta }^2+(4\pi de\beta )^2/15}}5^{}K.$$
(13)
Combining the particle-particle and impurity contributions and using the ideal gas equation of state: $`E(T)=AnT`$, we find the total cooling rate:

$$\frac{dT}{dt}=(4(1+1/g)n+n_{imp})\frac{85}{9m}\frac{\mathrm{}}{L^4}\frac{T^3}{T_{xph}^2}.$$
(14)
Integrating the Eq.(13) we get the temperature relaxation law: $`T(t)1/\sqrt{t}`$.

In the intermediate-T region at $`g=1`$ the cooling rate is enhanced by the stimulated scattering into the BEqC final states:

$$\frac{dT}{dt}=n_s\left(64n_s(1\frac{\zeta (3)}{\zeta (2)})+n_{imp}\frac{\zeta (3)}{\zeta (2)}\right)\frac{2\pi \mathrm{}^3}{m^2}\frac{\mathrm{log}^2L}{L^4}\frac{T^2}{T_{xph}^2}.$$
(15)
In the low-T region the thermodynamic equation of state reads: $`E(T)=A\zeta (3)T^3/\pi s^2`$, where $`s`$ is the Bogolubov sound velocity: $`s=\sqrt{\mu /m}`$. In order to calculate the energy losses we apply a Bogolubov unitary transformation to the Hamiltonian (1):

$$H_{ph}=\underset{\stackrel{}{p}\stackrel{}{q}}{}\frac{\mu }{2ϵ(\stackrel{}{p})}\psi ^+(\stackrel{}{p})\psi ^+(\stackrel{}{p}+\stackrel{}{q})\mathrm{\Gamma }_i(\stackrel{}{q})u_i(\stackrel{}{q})+c.c.$$
(16)
This Hamiltonian allows an emission of a phonon. The cooling rate in this case is also enhanced by the condensate stimulation:

$$\frac{dT}{dt}=\frac{\mu ^2}{\mathrm{}\pi T_{xph}^2}\left(1\frac{\zeta (4)}{\zeta (3)}\right)T^2.$$
(17)
For experimental realization of an exciton cooling the most relevant is Eq.(14). Integrating it and inserting $`L6`$ and $`nn_{imp}`$ we find the overall cooling time $`t_c`$ required to reach the phase transition point. Here one can distinguish two cases: i) cooling of exciton gas with constant density e.g. sustained by photoexcitation and ii) cooling of decaying exciton gas with $`n(t)=n_0\mathrm{exp}(t/\tau _r)`$, where $`\tau _r`$ is the exciton recombination time. We find:

$$t_c=\frac{C\mathrm{}T_{xph}^2}{T_c(n)^3},$$
(18)
where $`T_c(n)`$ is the BKT temperature as a function of the exciton density $`n`$ (6), $`C10`$ in the case i) and $`C30`$ in the case ii). Note that $`C`$ does not depend on exciton mass $`m`$ and in the case ii) the best condition for reaching $`T_c`$ occurs at $`t=\tau _r/3`$. Eq.(18) defines a frontier in the bilayer parameter space: $`(n,\tau _r)`$ or equivalently $`(T_c(n),\tau _r)`$, separating the two kinetic phases - one that can condensate and the second that remains above $`T_c`$ during the exciton life time $`\tau _r`$. Fig.2 shows this borderline for the case of a thin GaAs bilayer.

This work was supported by RFFI and INTAS.