text:
# Formation of the First Supermassive Black Holes

## 1. INTRODUCTION

Supermassive black holes (SMBHs) are believed to provide the power source of quasars via the accretion of surrounding gas (Salpeter 1964; Zeldovich 1964; Lynden-Bell 1969; Rees 1984). They were able to form already at early cosmic times, as implied by the recent discovery of quasars at redshifs $`z6`$ (Becker et al. 2001; Djorgovski et al. 2001; Fan et al. 2002). The existence of SMBHs with inferred masses of $`10^9M_{}`$, less than a billion years after the big bang, provides important constraints on any formation scenario (Haiman & Loeb 2001).

The gas physics involved in the formation of SMBHs is still not well understood (Rees 1984, 2002; Loeb & Rasio 1994; Barkana & Loeb 2001). In this work, we attempt to simulate this process under the simplest and best prescribed set of initial conditions; those dictated by the early universe. The earliest galaxies are simple in that they are the first condensations of gas to grow out of the seed inhomogeneities in the early universe. The composition of the primordial gas is determined by big bang nucleosynthesis (Burles, Nollett, & Turner 2001, and references therein) and any primordial magnetic fields are not expected to be dynamically significant.

Previous numerical simulations of early luminous structures (Bromm, Coppi, & Larson 1999, 2002; Abel, Bryan, & Norman 2000, 2002; Nakamura & Umemura 2001) have focused on the formation of stars inside the very first gaseous objects with a mass, $`10^5M_{}`$, just above the cosmological Jeans mass. Due to their low virial temperature (hundreds of K), fragmentation into stars is possible inside these objects through the formation of molecular hydrogen, H2, which cools efficiently via rotational-vibrational transitions even at these low temperatures. However, a relatively modest UV flux is sufficient to photo-dissociate the fragile H2 molecules and suppress their role in cooling the gas (Haiman, Rees, & Loeb 1997; Ciardi et al. 2000; Haiman, Abel, & Rees 2000). A destructive flux of UV photons could be produced by a small, early population of stars.

Here we consider dwarf galaxies with virial temperatures $`10^4`$ K inside of which atomic cooling is effective (Oh & Haiman 2002). These galaxies comprise a substantial population at $`z10`$ (Barkana & Loeb 2001) that is expected to dominate the reionization of hydrogen at $`z6`$ (Wyithe & Loeb 2003a). We use smoothed particle hydrodynamics (SPH) simulations to describe the cooling and dynamics of the gas for different choices of its total angular momentum (or spin parameter), following the suggestion by Eisenstein & Loeb (1995) that low-spin systems would be more susceptible to the formation of a SMBH. If H2 cooling is suppressed inside these galaxies, then their gas will not cool below $`10^4`$K. When the temperature to which the gas cools is only somewhat lower than the virial temperature of the host galaxy (at which the pressure force balances gravity for the gas), we expect that fragmentation into small clumps will be avoided and the gas will tend to condense isothermally (at its temperature floor of $`10^4`$K) into large clumps. Such large clumps may then collapse to form a SMBH, possibly through an intermediate stage of a supermassive star. The viability of this scenario relies on the suppression of molecular H2 cooling, which when present is capable of cooling the gas to a temperature as low as 200 K. We therefore start our analysis in § 2 by considering the cosmic UV background expected from the first stars. Later on, in § 5 we carry out numerical simulations to determine the level of UV radiation necessary for the photodissociation of H2 molecules during the early collapse phase of the above dwarf galaxies.

The SPH code we use is the same as was used to simulate the formation of the first massive stars (Bromm et al. 1999, 2002), except that we apply it to more massive galaxies where atomic cooling is effective. Our goal is to find whether there is a direct route to the formation of SMBHs under these circumstances. Although our focus is on gas dynamics, we note that the first massive black holes (BHs) may have also formed out of the gravitational dynamics of clusters of stellar-mass BHs (Larson 2000, 2002; Madau & Rees 2001; Schneider et al. 2002; Islam, Taylor, & Silk 2003; Volonteri, Haardt, & Madau 2003).

Throughout this paper, we assume a standard $`\mathrm{\Lambda }`$CDM cosmology, with a total density parameter in matter of $`\mathrm{\Omega }_m=1\mathrm{\Omega }_\mathrm{\Lambda }=0.3`$, and in baryons of $`\mathrm{\Omega }_\mathrm{B}=0.045`$. The Hubble constant is $`h=H_0/(100`$ km s-1 Mpc-1)=0.7, and the present-day power-spectrum amplitude is $`\sigma _8=0.9`$ in spheres of radius $`8h^1`$Mpc.

## 2. THE COSMOLOGICAL CONTEXT

We first outline the basic physical reason for why small protogalactic systems of virial temperature $`10^4`$K and mass $`10^8M_{}`$, collapsing before the epoch of reionization at redshifts $`z10`$, provide intriguing sites for the formation of the first SMBHs.

The direct collapse of a primordial gas cloud into a central compact object is made difficult by fragmentation and consequent star formation (Loeb & Rasio 1994). The ability of the gas to fragment depends on the presence of an efficient cooling mechanism. By a redshift of $`z10`$, the intergalactic medium (IGM) is expected to be enriched with heavy elements from the first generation of stars to a level of $`Z10^{3.5}Z_{}`$ (e.g., Gnedin & Ostriker 1997; Mackey, Bromm, & Hernquist 2003). Recent numerical investigations have shown that such a metallicity is sufficient to enable the gas to efficiently cool, fragment, and subsequently form stars (Omukai 2000; Bromm et al. 2001a; Bromm & Clarke 2002). It is plausible that in some regions of the IGM, the gas that collapses into a dwarf galaxy has not been enriched with metals in excess of $`Z_{\mathrm{crit}}10^{3.5}Z_{}`$ (e.g., Thacker, Scannapieco, & Davis 2002; Furlanetto & Loeb 2003). In this paper, we investigate the collapse of a virtually metal-free system of this type.

It is well known that molecular hydrogen can effectively cool metal-free gas, and thus enable the formation of the first stars at $`z2030`$ (e.g., Haiman, Thoul, & Loeb 1996; Tegmark et al. 1997). Even in the absence of metals, star formation is therefore expected to occur, and the subsequent negative feedback due to supernova explosions would prevent the assembly of large quantities of gas in the center of the shallow dark matter (DM) potential well that characterizes the first dwarf galaxies (e.g., Dekel & Silk 1986). Molecular hydrogen, however, is fragile and can readily be destroyed by photons in the Lyman-Werner (LW) bands, within the energy range $`11.213.6`$ eV, via the two-step Solomon process (Stecher & Williams 1967)

$$\mathrm{H}_2+\gamma \mathrm{H}_2^{}2\mathrm{H}.$$

The intermediate stage involves an excited electronic state, $`\mathrm{H}_2^{}`$, from which a fraction of the subsequent decays end in the vibrational continuum of the ground state, resulting in the dissociation of the molecule.

The question then arises whether H2 cooling can indeed be suppressed in the pre-reionization dwarf galaxies considered in this paper. These systems are rather massive compared to the $`10^510^6M_{}`$ halos that host the formation of the first stars at $`z20`$. The gas might then be able to self-shield against the photo-dissociating LW background (e.g., Glover & Brand 2001; Machacek, Bryan, & Abel 2001). Close to the epoch of reionization, however, a significant flux in the LW bands ($`h\nu <13.6`$ eV) is expected. To estimate the LW flux, we first consider the flux of ionizing radiation just above the Lyman limit,

$$J_\nu ^+\frac{hc}{4\pi }\frac{N_\gamma X\rho _\mathrm{B}(z=10)}{m_\mathrm{H}}\text{ .}$$
(1)
Here, $`\rho _\mathrm{B}(z=10)`$ is the baryonic density at $`z10`$, $`X=0.76`$ the mass fraction in hydrogen, and $`m_\mathrm{H}`$ the mass of a hydrogen atom. We assume that $`N_\gamma 10`$ ionizing photons per baryon are required to reionize the universe (Wyithe & Loeb 2003a). Evaluating this expression for the cosmological parameters given above, we find: $`J_{21}^+40`$. Here, and in the remainder of the paper, UV fluxes are normalized as $`J_\nu =J_{21}\times 10^{21}`$ erg cm-2 s-1 Hz-1 sr-1.

The flux in the LW bands just below the Lyman limit, $`J_\nu ^{}`$ could, however, be much larger than $`J_\nu ^+`$. Assuming that only a fraction, $`f_{\mathrm{esc}}`$, of the ionizing photons can escape from the star forming halos, we have: $`J_{21}^{}J_{21}^+/f_{\mathrm{esc}}4\times 10^3(f_{\mathrm{esc}}/0.01)`$, where $`f_{\mathrm{esc}}`$ is expected to be low at high redshifts (Wood & Loeb 2000, and references therein)111A strong background of photo-dissociating photons could also result from star formation in the dwarf galaxy itself (e.g., Omukai & Nishi 1999; Oh & Haiman 2002). However, this would also result in the production of metals, probably leading to the rapid enrichment of the galaxy beyond $`Z_{\mathrm{crit}}`$..

The strong UV background flux estimated above makes it possible for H2 formation to be effectively suppressed close to the redshift of reionization. In the absence of molecular hydrogen, however, cooling can still proceed via atomic transitions in halos of mass (see Barkana & Loeb 2001)

$$M10^8M_{}\left(\frac{1+z}{10}\right)^{3/2}\text{ .}$$
(2)
The virial temperature, $`T_{\mathrm{vir}}10^4`$K, in these more massive halos allows for the very efficient cooling of the gas via lines of atomic hydrogen. Notice that in this case where the gas temperature is close to the virial temperature, the gas cloud as a whole can undergo collapse but it will not be able to fragment until high enough densities are reached so that the Jeans mass has declined sufficiently. A halo of total mass $`10^8M_{}`$ and a collapse redshift $`z10`$ corresponds to a 2$`\sigma `$ peak in the random field of primordial density fluctuations.

In § 4, we describe numerical simulations of such a 2 $`\sigma `$ peak. We consider the evolution with and without the presence of H2. To ascertain the importance of H2, we consider the limiting cases of a halo made of purely atomic hydrogen, and of a halo where H2 is allowed to form without any negative feedback. Subsequently, in § 5 we discuss simulations where an external LW background is included. We in particular assess the critical LW flux required to prevent molecules from forming. This required background flux may then be compared to the available flux at $`z10`$, as estimated in this section.

## 3. Numerical Methodology

### 3.1. Large-scale Simulations

The evolution of the dark matter and gas components is calculated with a version of TREESPH (Hernquist & Katz 1989), combining the SPH method with a hierarchical (tree) gravity solver (see Bromm et al. 2002 for further details). Here, we briefly describe the additions to the code which are necessary for the investigation of zero-metallicity gas. These include a method to treat the radiative cooling of the gas together with a multi-species chemical reaction network, and a technique to create sink particles. The thermal evolution of the gas is governed by the equation:

$$\frac{\mathrm{D}u}{\mathrm{D}t}=\frac{P}{\rho ^2}\frac{\mathrm{D}\rho }{\mathrm{Dt}}+\frac{\mathrm{\Gamma }\mathrm{\Lambda }}{\rho }$$
(3)
where $`\mathrm{D}/\mathrm{D}t`$ is the Lagrangian time derivative, $`P`$ and $`\rho `$ are the gas pressure and density, $`u`$ is the specific internal energy (in erg g-1), and $`\mathrm{\Gamma }`$ and $`\mathrm{\Lambda }`$ are the contributions from radiative heating and cooling, respectively (in erg cm-3 s-1). Since we are concerned with a system that formed prior to the epoch of reionization, we assume that there is no heating due to a photoionizing background ($`\mathrm{\Gamma }=0`$). Radiative cooling is due to lines of atomic hydrogen at gas temperatures $`T10^4`$ K, and lines of molecular hydrogen at lower temperatures. We have implemented the H2 cooling function given by Galli & Palla (1998). Finally, the first term on the right-hand side describes adiabatic cooling due to expansion or heating due to compression.

Since radiative cooling to temperatures below that of the cosmic microwave background (CMB), $`T_{\mathrm{CMB}}=2.7\text{ K}(1+z)`$, is thermodynamically not possible, we write the cooling term as

$$\mathrm{\Lambda }=\mathrm{\Lambda }(T)\mathrm{\Lambda }(T_{\mathrm{CMB}})\text{ .}$$
(4)
For $`T<T_{\mathrm{CMB}}`$, radiative cooling consequently turns into heating. This approximate treatment ensures that $`TT_{\mathrm{CMB}}`$, unless cooling proceeds via adiabatic expansion.

The chemical reaction network comprises the 9 species H, H+, H-, H2, H$`{}_{}{}^{+}{}_{2}{}^{}`$, $`e^{}`$, He, He+, and He++, including the reactions given in Haiman et al. (1996). Our implicit, backwards-differencing method to solve the coupled set of rate equations is fast and accurate (see Bromm et al. 2002 for test calculations).

We have devised an algorithm to merge SPH particles in high-density regions in order to overcome the otherwise prohibitive time-step limitation, as enforced by the Courant stability criterion. To follow the simulation for a few dynamical times, we allow SPH particles to merge into more massive ones, provided they exceed a predetermined density threshold. More details of the merging algorithm are given in Bromm et al. (2002).

### 3.2. Small-scale Simulations

The sink-particle technique allows one to study the collective dynamics of multiple centers of condensation, such as merging between clumps and ongoing accretion over many dynamical timescales (e.g., Bate, Bonnell, & Price 1995; Bate, Bonnell, & Bromm 2003). It is important, however, to also be able to follow the collapse to increasingly high density. To accomplish this, we have carried out a complementary simulation in which we do not create sink particles. Instead, we focus on the highest-density region in one of the large-scale simulations, resample the density field with an increased number of SPH particles, and adopt a rapidly decreasing timestep according to $`\mathrm{\Delta }t_{\mathrm{sys}}1/\sqrt{G\rho _{\mathrm{max}}}`$. Here, $`\rho _{\mathrm{max}}`$ is the maximum gas density at a given instant, and $`\mathrm{\Delta }t_{\mathrm{sys}}`$, the system timestep, is the maximum allowed time by which a particle is advanced within the multiple-timestep scheme employed in the simulations (see Hernquist & Katz 1989). In doing so, the overall system does hardly evolve at all, while the runaway collapse of one small region proceeds on a sufficiently short timescale. The mass resolution of a simulation is approximately

$$M_{\mathrm{res}}=\left(\frac{N_{\mathrm{neigh}}}{N_{\mathrm{SPH}}}\right)M_\mathrm{B}\text{ ,}$$
(5)
where $`N_{\mathrm{neigh}}32`$ is the number of particles within a given SPH smoothing kernel, $`N_{\mathrm{SPH}}`$ the total number of SPH particles, and $`M_\mathrm{B}`$ the total baryonic mass. To avoid numerical fragmentation, the Jeans scale has to be resolved (Bate & Burkert 1997): $`M_{\mathrm{res}}<M_J`$. After the onset of gravitational instability in the large-scale simulation, the rapid increase in density leads to a violation of this criterion. At this point, a sink particle is created in the simulation.

In our refined simulation, we follow a different computational strategy (see also Bromm 2000). Our starting configuration is the large-scale simulation, stopped at a moment briefly before the central region undergoes runaway collapse. We now apply the following resampling procedure to the fluid within this region. Every SPH particle in the original, unrefined simulation acts as a parent particle denoted by $`p`$, and spawns $`N_{\mathrm{ref}}`$ child particles denoted by $`k`$, where $`N_{\mathrm{ref}}=50`$ in this paper. The child particles are distributed according to the SPH smoothing kernel $`W(\stackrel{}{r}_k\stackrel{}{r}_p;h_p)`$, by employing a standard Monte-Carlo comparison-rejection method (see, e.g., Press et al. 1992). Here, $`h_p`$ is the smoothing length of the parent particle. The velocity, temperature, and fractional abundances of the 9 species included in the reaction network are directly inherited from the parent particle, $`\stackrel{}{v}_k=\stackrel{}{v}_p`$, $`T_k=T_p`$, and $`y_k=x_p`$ for all $`k`$, respectively. Finally, each child particle is assigned a mass of $`m_k=m_p/N_{\mathrm{ref}}`$. This procedure conserves linear and angular momentum well. Energy is also conserved well, although there arises a small artificial contribution to the gravitational potential energy due to the discreteness of the resampling. The resampling described in § 4.2 results in $`N_{\mathrm{SPH}}=`$ 65,500 within the central high-resolution region, and the mass resolution is now $`M_{\mathrm{res}}500M_{}`$, as compared to $`M_{\mathrm{res}}20,000M_{}`$ in the original simulation.

As we follow our simulation to densities in excess of $`n_\mathrm{H}10^7`$ cm-3, opacity effects begin to modify the thermal behavior of the gas. For simplicity, we continue to use the optically thin rate for cooling due to the collisional excitation of atomic hydrogen lines (Cen 1992). While this is clearly an oversimplification, we have verified that the resulting thermal evolution of the gas is very similar to what is found in a more sophisticated treatment of the cooling processes at high densities (Omukai 2001).

Our technique of refining a coarser, parent simulation, and following the further collapse with increased resolution, is conceptually similar to the Adaptive Mesh Refinement (AMR) method which was originally developed by Berger & Oliger (1984). The AMR method has already been successfully applied to astrophysical problems (e.g., Truelove et al. 1998; Norman & Bryan 1999), but our approach is one of the first attempts of implementing such a scheme within SPH (see also Kitsionas & Whitworth 2002).

### 3.3. External UV background

To investigate the effect of an external UV background on the evolution of the gas, we have included the relevant photo-reactions into our chemical network. The photo-rates (in units of s-1) are given by

$$k_{\mathrm{ph}}=4\pi _{\nu _{\mathrm{th}}}^{\mathrm{}}\frac{\sigma _\nu J_\nu \text{d}\nu }{h\nu }\text{ ,}$$
(6)
where $`h\nu _{\mathrm{th}}`$ is the threshold energy for a specific reaction. The cross-sections are given in Haiman et al. (1996). Following Omukai (2001), we consider both a power-law spectrum, $`J_\nu \nu ^1`$, appropriate for quasar sources, and a thermal Planck spectrum, $`J_\nu B_\nu (T_{})`$, appropriate for stellar sources. In the latter case, we evaluate the photo-rates for two different values of the radiation temperature: $`T_{}10^4`$ K and $`10^5`$ K, corresponding to ‘normal’ Population II stars and to very massive Population III stars, respectively (Bromm, Kudritzki, & Loeb 2001b). For each spectrum, we normalize the flux at the Lyman limit with a normalization constant $`J_{21}`$, and assume zero flux at $`h\nu >13.6`$ eV due to intergalactic HI absorption.

To take into account the effect of self-shielding, we write the H2 photodissociation rate as $`k_{\mathrm{diss}}J_{21}f_{\mathrm{shield}}`$, with a proportionality constant that depends on the chosen spectrum. For the shielding factor, we use the approximate expression $`f_{\mathrm{shield}}\mathrm{min}[1,(N_{\mathrm{H}_2}/10^{14}\mathrm{cm}^2)^{0.75}]`$ (Draine & Bertoldi 1996). This formula is accurate only for a static medium, and it will grossly overestimate the H2 line opacity in the presence of large-scale velocity gradients. In our simulations, we encounter moderately strong bulk flows, with velocities of order the sound speed. We therefore conservatively overestimate the effect of self-shielding. The H2 column density is estimated from local quantities only as follows: $`N_{\mathrm{H}_2}0.1f_{\mathrm{H}_2}n_\mathrm{H}L_{\mathrm{char}}`$. We define a local characteristic length, $`L_{\mathrm{char}}`$, such that the total baryonic mass would be contained in a sphere of uniform density $`n_\mathrm{H}`$, and radius $`L_{\mathrm{char}}`$. Thus, we have

$$L_{\mathrm{char}}=\left(\frac{3XM_\mathrm{B}}{4\pi m_\mathrm{H}n_\mathrm{H}}\right)^{1/3}\text{ .}$$
(7)
By comparing with the column densities that are actually realized in the simulations, we have verified the validity of this prescription.

We neglect any contributions to the heating of the gas due to photo-reactions, since these are always negligible in the system studied here. We also neglect, for simplicity, the sawtooth-modulation of the background spectrum due to the IGM line-opacity in the Lyman series (Haiman et al. 1997).

### 3.4. Initial Conditions

Within the hierarchical $`\mathrm{\Lambda }`$CDM model, the first luminous objects are expected to form out of high $`\sigma `$ peaks in the random field of primordial density fluctuations. The early (linear) evolution of such a peak, assumed to be an isolated and roughly spherical overdensity, can be described by the top-hat model (e.g., Padmanabhan 1993). We use the top-hat approximation to specify the initial DM configuration. In this paper, we investigate the fate of a 2 $`\sigma `$ peak of total mass $`10^8M_{}`$, corresponding to $`1.5\times 10^7M_{}`$ in baryons. On this mass scale, one can estimate the redshift of collapse (or virialization) to be $`z_{\mathrm{vir}}10`$.

Our simulation is initialized at $`z_i=100`$, by performing the following steps. The collisionless DM particles are placed on a cubical Cartesian grid, and are then perturbed according to a given power spectrum $`P(k)=Ak^n`$, by applying the Zeldovich (1970) approximation which also allows to self-consistently assign initial velocities. The power-law index is set to $`n=2.5`$ which approximately describes the spectral behavior on the scale $`10^8M_{}`$. To fix the amplitude $`A`$, we specify the initial variance of the fluctuation amplitude

$$\sigma _i^2=Ak^n\text{ .}$$
(8)
The summation is over all contributing modes, where the minimum wavenumber is given by the overall size of the Cartesian box, and the maximum wavenumber by the Nyquist frequency. Choosing $`\sigma _i^20.01`$, the rms amplitude of fluctuations at the moment of collapse is approximately

$$\sigma (z=10)\left(\frac{1+z_i}{1+z}\right)\sigma _i1\text{ .}$$
(9)
This choice ensures that the substructure develops on a similar timescale as the overall collapse of the background medium. Next, particles within a (proper) radius of $`R_i=`$ 776 pc are selected for the simulation. The resulting number of DM particles is $`N_{\mathrm{DM}}=17074`$. Finally, the particles are set into rigid rotation and are endowed with a uniform Hubble expansion (see also Katz 1991). Angular momentum is added by assuming a spin parameter $`\lambda =L|E|^{1/2}/(GM^{5/2})=0`$ and 0.05, where $`L`$, $`E`$, and $`M`$ are the total angular momentum, energy, and mass, respectively. The spin parameter is a measure of the degree of rotational support, such that the ratio of centrifugal to gravitational acceleration is given by $`\lambda ^2`$ at virialization. The second value corresponds to the average spin parameter found in cosmological simulations (e.g., Barnes & Efstathiou 1987; Jang-Condell & Hernquist 2001).

The collisional SPH particles ($`N_{\mathrm{SPH}}=32768`$) are randomly placed to approximate a uniform initial density. The SPH particles follow the same Hubble expansion and rigid rotation as the DM particles. For the initial gas temperature at $`z_i=100`$ we adopt the value of $`200`$ K (Tegmark et al. 1997). The fractional free-electron abundance is initialized as $`x_e=4.6\times 10^4`$, and the hydrogen molecule abundance, in the runs where H2 formation is allowed, as $`f_{\mathrm{H}_2}=2\times 10^6`$ (Anninos & Norman 1996).

In Table 1, we summarize the parameters of the large-scale simulations. We will discuss the initial state of the refined, small-scale simulation in § 4.2. We empasize that the initial particle setup, and in particular the realization of the DM fluctuations, is identical for all runs which only differ in the choice of spin, in the presence or absence of a soft UV background, and in whether H2 cooling is allowed or not. The runs with no H2 cooling correspond to the limiting case in which all the H2 has been radiatively destroyed by a soft UV background. In § 2, we have discussed the possible emergence of such an efficient feedback following the build-up of the cosmic UV background close to the epoch of reionization.

## 4. The Simulations

In this section, we only consider simulations where no external UV background is included, and defer the discussion of its effect to § 5. We carry out these simulations in two separate steps. First, we simulate the large-scale evolution on the scale of the host halo, $`1`$ kpc, ending with the formation of high-density clumps on a scale of $`1`$ pc. Once the mass resolution of the simulation becomes larger than the local Jeans mass, $`M_JM_{\mathrm{res}}`$, a sink particle is created. Consequently, the internal dynamics of a clump cannot be studied any further in these coarse-grain simulations. Second, we follow the evolution of a clump to higher densities and smaller spatial scales, by refining the resolution in the vicinity of the clump to again resolve the Jeans mass. We begin with a discussion of the large-scale evolution.

### 4.1. Large-scale evolution

We first describe the run with no spin and with suppressed H2 formation (Run A) in more detail, and subsequently explore what happens under variations of these assumptions. In Figure 1, we show the initial configuration at $`z_i=100`$, which is identical for all runs. The overdense region initially expands with the Hubble flow at a reduced rate, subsequently turns around at $`z_{\mathrm{ta}}16`$, and eventually collapses.

Figure 2 shows the system at $`z13`$, briefly after turnaround. The DM component has developed a marked substructure in response to the imprinted perturbations, and the baryons have begun to fall into the deepest DM potential wells. Eventually, close to the redshift of virialization, $`z_{\mathrm{vir}}10`$, the DM is undergoing violent relaxation (Lynden-Bell 1967), resulting in an approximate balance between kinetic and gravitational potential energy.

The baryons, on the other hand, dissipatively settle into the center of the DM halo (see Figure 3). Focusing on the innermost $`200`$ pc of the simulation box, Figure 4 displays the central distribution of the gas density. At this stage in the simulation, a high-density sink particle of mass $`3\times 10^6M_{}`$ and radius $`1`$ pc has been created. Subsequently, the sink particle grows in mass by continuous accretion from the surrounding diffuse gas.

To gain further insight into the physics of the simulation, we consider the thermal and chemical properties of the gas, as illustrated in Figure 5. We show the free-electron abundance, temperature, and Jeans mass as a function of gas density for every SPH particle. This manner of presentation contains an additional dimension of information, concerning the overall timescale of evolution: when the evolution proceeds slowly, particles tend to pile-up in the respective plots, whereas only a few particles populate regions in parameter space where the evolution is fast. As can be seen in panel (d), a gas density of $`n10^2\mathrm{cm}^3`$ and a Jeans mass of $`M_J10^6M_{}`$ approximately mark the transition between a phase of slow and fast evolution. Physically, this transition corresponds to the onset of runaway collapse when the clump mass is $`MM_J`$.

In Figure 6, we show the central gas configuration for the simulation where the initial spin of the system is nonzero (Run B). In this case, a binary system of clumps has formed with a separation of $`10`$ pc. These two clumps individually have properties similar to the one clump formed in Run A. The binary is rapidly drawn into the center of the halo by dynamical friction which operates on a timescale $`t_{\mathrm{df}}0.1r^2v_c/(GM_{\mathrm{Cl}})10^7`$ yr (Binney & Tremaine 1987). Here, $`v_c`$ is the circular velocity at radius $`r`$, and $`M_{\mathrm{Cl}}`$ is the clump mass. This timescale is short compared to the Hubble time at $`z10`$, $`t_H10^8`$ yr. Such a system of two compact objects is expected to efficiently radiate gravitational waves that could be detected with the planned Laser Interferometer Space Antenna222See http://lisa.jpl.nasa.gov/ (LISA) (Wyithe & Loeb 2003b).

Massive clumps, either in isolation or as a binary, are able to form in both runs A and B, which correspond to different values of initial spin. The result that such a compact object can form even in a halo with average initial spin (run B), is seemingly at odds with the earlier investigation by Eisenstein & Loeb (1995), who have argued that only the lowest-spin cosmological perturbations can harbor seed BHs. The gas that is eventually incorporated into the two massive clumps in run B must have been able to efficiently lose angular momentum. It is crucial that the evolution in the $`\lambda =0.05`$ case leads to the formation of a binary system. Tidal torques can then transfer much of the angular momentum from the gas around each forming clump to the orbital motion of the system, similar to the case of present-day star formation in a clustered environment (e.g., Larson 2002; Bate et al. 2003). We plan to more fully address the complex physics of angular momentum transport in virializing DM halos in future work.

This efficient mechanism of assembling a large amount of gas into a compact ($`1`$ pc) configuration crucially depends on the suppresion of star formation in the collapsing halo. In Run C, which otherwise has the same initial conditions as Run A, we now allow for the formation of hydrogen molecules. Figure 7 (which should be compared to Fig. 5 of Run A), shows that the thermal evolution of the gas in this case proceeds very differently from Runs A and B. Here, the gas is able to cool to temperatures $`T200`$ K due to the presence of H2. Provided that H2 formation is not suppressed, the properties of primordial, metal-free gas are rather similar in systems of mass $`10^8M_{}`$, and in the smaller halos of mass $`10^6M_{}`$ that are predicted to host the very first stars at $`z2030`$ (see Fig. 10 in Bromm et al. 2002).

The thermal evolution is clearly reflected in the morphology. Figure 8 shows the projected gas density. Due to efficient cooling, the gas can readily fragment and subsequently undergo runaway collapse to form stars already shortly after turnaround. This star formation activity will tend to prevent the assembly of gas in compact clumps in two ways. First, the Population III stars are expected to be short-lived, with lifetimes of $`3\times 10^6`$ yr, and they will exert a strong negative feedback on their surroundings upon their death as supernovae (SNe). Second, the SNe will disperse the nucleosynthetic products from the first generation of stars into the remaining gas of the halo. Henceforth, the enriched gas will be able to cool and form stars regardless of whether H2 is present or not. Since the photodissociation of H2 is now no longer able to limit the gas condensation (see § 2), the star formation rate is expected to increase substantially (e.g., Nishi & Tashiro 2000; Bromm & Clarke 2002).

The massive clumps formed in runs A and B are evidently not yet BHs. To examine the subsequent dynamics of the clump, we have resimulated the evolution of a typical clump with higher resolution in order to reach smaller spatial scales. We report on this simulation next.

### 4.2. Small-scale evolution

As a representative example, we have selected Run A for initializing the refined simulation. We focus on the central region, as shown in Figure 4, and resample the gas within a radius of $`5`$ pc from the density maximum. The initial moment of the resampled run is chosen to be somewhat earlier than in Figure 4, corresponding to $`z10.6`$, slightly before a sink particle was created in the coarse-grain simulation. In the resampled region, the baryons comprise a mass of $`6\times 10^5M_{}`$, and are now the dominant component, compared to only $`10^5M_{}`$ in dark matter.

After $`10^4`$ yr, we reach the end-state of our refined simulation, at which point the Jeans mass again becomes comparable to the resolution limit. The central gas cloud is in a state of free-fall, and we do not see any signs of further sub-fragmentation. We have verified that fragmentation is not artificially suppressed in the resampling process. To this extent, we have carried out a fiducial simulation that does lead to fragmentation at some point in the evolution, and a comparison calculation where we resample the fluid before fragmentation occurs. We find that the resulting fragmentation pattern is very similar in the two cases, as is physically expected.

The innermost region, of size $`0.1`$ pc, comprises $`10^4M_{}`$ in gas with densities in excess of $`10^9`$ cm-3. If we were to continue the simulation further in time, the amount of gas residing inside the central $`0.1`$ pc would rapidly increase. Indeed, the mass of the sink particle formed in the large-scale simulation of Run A, of order a few $`10^6M_{}`$, is indicative of the total baryonic mass that will end up in the central compact object. At this stage, this object is characterized by a ratio $`E_{\mathrm{rot}}/|E_{\mathrm{grav}}|0.5`$, where $`E_{\mathrm{rot}}`$ and $`E_{\mathrm{grav}}`$ are the rotational and gravitational energies of the central compact object, respectively. The angular momentum present in the central clump must have arisen through torques from the clumpy DM distribution during the relaxation process, as run A has zero initial spin. Again, we plan to revisit this issue in our planned work on angular momentum transport in collapsing DM halos.

The crucial question now is: What is the further fate of the central object? In particular, do we expect the cloud to fragment, and eventually form a central stellar cluster? Or else, will fragmentation be inhibited, as in the earlier, large-scale evolution described in the previous section? In the latter case, the additional question arises (Loeb & Rasio 1994) whether the gas cloud will continue to collapse directly to a black hole or else settle temporarily into a pressure-supported configuration of a supermassive star (see §17 in Shapiro & Teukolsky 1983)? We can gain important physical insight into these questions by comparing three important timescales in the problem, namely: the free-fall time, $`t_{\mathrm{ff}}`$, during which dynamical equilibrium is established; the cooling time, $`t_{\mathrm{cool}}`$, during which the gas may radiate its thermal energy as it contracts and heats to higher temperatures; and the viscous time, $`t_{\mathrm{vis}}`$, during which angular momentum is transported. The cooling time remains relatively short and allows continued collapse until the cloud becomes optically thick to the radiation it produces. Spherically-symmetric calculations indicate that when the particle density in the contracting gas clump rises above $`10^{17}\mathrm{cm}^3`$, the ionization fraction increases sharply by 7 orders of magnitude to a value $`x_e0.1`$ (see Fig. 4 in Omukai 2001). A clump of mass $`2\times 10^6M_{}`$ achieves this dense phase when its radius $`R`$ shrinks to $`(R/10^{16}\mathrm{cm})R_{16}0.5`$, about two orders of magnitude smaller than the value at the end of our simulations. When $`R_{16}1`$, the optical depth of the cloud to Thomson scattering by free electrons is $`\tau 10^7R_{16}^2`$ and so the cooling time is given by the photon diffusion time $`t_{\mathrm{cool}}\tau R/c10^5\mathrm{yr}R_{16}^1`$. The free fall time, $`t_{\mathrm{ff}}(G\rho )^{0.5}1\mathrm{yr}R_{16}^{1.5}`$, is much shorter than this cooling time and so the gas may settle into a radiation-pressure supported configuration resembling a rotating supermassive star (SMS). The viscous time is shorter than the cooling time, $`t_{\mathrm{vis}}4\times 10^4\mathrm{yr}R_{16}^{0.5}(P_{\mathrm{gas}}/P_{\mathrm{tot}})\alpha ^1`$, where $`\alpha 1`$ is a dimensionless measure of the viscosity coefficient (Shakura & Sunyaev 1973; §14.5 in Shapiro & Teukolsky 1983) in terms of the total (gas$`+`$radiation) pressure, $`P_{\mathrm{tot}}=(P_{\mathrm{gas}}+P_{\mathrm{rad}})P_{\mathrm{gas}}`$. The resulting viscous dissipation will likely heat the optically-thick gas close to its virial temperature since the rotational energy is substantial at the end of our simulation. Combined with additional heating from adiabatic compression, the collapse may be halted as soon as the gas becomes highly optically-thick and is unable to cool efficiently. We therefore conclude that if the simulated gas clouds continue to shrink during the early optically-thin (efficient cooling) regime by about two orders of magnitude in radius (as they do earlier in the simulations), then they will likely form a rotating SMS.

Recent fully-relativistic calculations of the evolution of rotating SMSs have predicted the collapse into a massive BH (Baumgarte & Shapiro 1999; Shibata & Shapiro 2002). Under a wide range of initial conditions, a substantial fraction of the mass of the SMS ($`90\%`$), is expected to end up in the BH.

## 5. EFFECT OF A PHOTODISSOCIATING BACKGROUND

We now return to the question: What is the required level of the LW background to prevent H2 from forming when self-shielding is taken into account? To answer this question, we have carried out a series of simulations where we allow H2 to form self-consistently in the presence of a soft UV background. We have determined $`J_{\mathrm{crit}}`$, the critical flux above which H2 is not able to form throughout the simulation volume, for the three spectra described in § 3.3. For both the power-law spectrum, as well as the thermal spectrum with $`T_{}=10^5`$ K, we find: $`J_{\mathrm{crit}}10^5`$. This value is larger than the expected background close to reionization (see § 2), and it would be difficult to prevent H2 cooling for these two spectra.

A much lower background level, however, suffices in the case of a thermal spectrum with $`T_{}=10^4`$ K, for which we find: $`J_{\mathrm{crit}}10^3`$, close to the predicted UV background at $`z10`$. In Figure 9 we show the timescales that are relevant to understand the abundance of H2 in the simulation with $`J_{21}=10^3`$. It is evident that even when self-shielding is included, the H2 formation time, $`t_{\mathrm{form}}`$, is longer than the H2 destruction time, $`t_{\mathrm{dest}}`$ throughout the simulation. Once the gas has reached a density of $`n_\mathrm{H}10^3\mathrm{cm}^3`$, collisional dissociation becomes very effective in destroying H2 at a gas temperature of $`10^4`$K. The ever-increasing column densities, and the resulting level of self-shielding, have therefore no effect on the molecule abundance at densities high enough for collisions to become important.

Using an idealized one-zone model to describe the dynamics of a collapsing primordial cloud, but otherwise implementing the relevant radiative processes in a sophisticated way, Omukai (2001) has determined the critical UV flux in the case of a Planck spectrum with $`T_{}=10^4`$ K to be: $`J_{21}10^3`$, in good agreement with our numerical simulations.

Both these estimates crucially depend on the presence of a stellar-like radiation background. The formation of the first quasars, therefore, might well require an earlier epoch of star formation. The old puzzle of whether quasars or stars were the first luminous objects in the universe (e.g., Silk & Rees 1998) would then be answered in favor of the latter.

## 6. SUMMARY AND CONCLUSIONS

Our numerical simulations show that metal-free dwarf galaxies at redshifts $`z10`$ whose cooling is dominated by atomic transitions (with a virial temperature just above $`10^4`$K) tend to form a central condensation consisting of one or two clumps which contain $`10\%`$ of the total baryonic mass of the galaxy. As long as H2 formation is suppressed, these massive clumps do not fragment but rather cool and continue to collapse isothermally at a temperature of $`10^4`$K. We have simulated the collapse to the point where $`10^6M_{}`$ of gas condense to a scale $`1`$pc. At the end of our simulation the clump maintains a nearly free-fall collapse and shows no signs of fragmentation. We expect the collapse to continue until the gas cloud becomes optically thick to Thomson scattering and its cooling time is much longer than the infall time. At this stage the collapse may be halted due to viscous dissipation of rotational energy, and a rotating SMS is likely to form. After shedding a small fraction of its mass, the SMS will collapse to a SMBH (Baumgarte & Shapiro 1999; Shibata & Shapiro 2002).

The possible existence of a supermassive star as an intermediate stage in the formation of the first quasars can be tested by direct observations. The SMS spectrum is expected to be close to thermal at an effective temperature of $`10^5`$ K (Bromm et al. 2001b). The SMS mass can be inferred from its bolometric luminosity, given that the latter equals the Eddington value, $`L_{\mathrm{Edd}}=1.4\times 10^{44}\mathrm{erg}\mathrm{s}^1(M/10^6M_{})`$. The emission redshift can be easily derived from the Gunn-Peterson trough produced by Ly$`\alpha `$ absorption of the neutral IGM (Gunn & Peterson 1965). Our model predicts that the emission spectrum of the SMS would show only hydrogen and helium lines but no metal lines (such as the broad emission lines observed for quasars at lower redshifts), since metal enrichment would signal small scale fragmentation into stars and would in turn enable the production of molecules on dust grains. However, all the above characteristics are common to clusters of Population III stars (Bromm et al. 2001b). To distinguish between a SMS and the extended image of a star cluster requires high angular resolution, possibly achievable with the James Webb Space Telescope333See http://ngst.gsfc.nasa.gov/. Alternatively, it is possible to make use of the high probability of gravitational lensing by intervening galaxies out to $`z10`$ (Barkana & Loeb 2000; Wyithe & Loeb 2002a,b; Comerford, Haiman, & Schaye 2002). In particular, gravitational microlensing by stars within the lens galaxies (Wyithe & Loeb 2002b) is only possible for a SMS and not for a star cluster whose extent is much larger than the Einstein radius of a solar-mass star at a cosmological distance, $`10^2\mathrm{pc}`$.

Our calculations show that fragmentation into stars is suppressed inside the above dwarf galaxies only if H2 cooling is negligible. The photodissociation of H2 molecules requires a background flux of UV photons below 13.6eV (to which the neutral IGM is nearly transparent) that could possibly be produced by stars at $`z10`$, prior to reionization (see §2).

What is the comoving density of the first generation of massive BHs? At $`z10`$, the comoving density of halos with a virial temperature $`10^4`$K is $`5\mathrm{Mpc}^3`$ (Sheth & Tormen 1999). Requiring that the total mass density of massive BHs would not exceed the value observed in the local universe (Yu & Tremaine 2002), $`2.5\times 10^5M_{}\mathrm{Mpc}^3`$, we infer that the average BH mass per halo must be $`5\times 10^4M_{}`$, much smaller than the characteristic mass of the gas clumps in our simulations. This follows naturally from our expectation that only a small fraction of all halos are both metal poor and exposed to a high level of UV flux during their entire history (see §2 and §5). Even though a substantial fraction of the IGM volume is not enriched with metals at $`z10`$ (Thacker et al. 2002; Furlanetto & Loeb 2003), most $`2\sigma `$ halos may have been enriched by stars that were formed inside their hierarchical building blocks at an earlier cosmic time, when the UV background was much lower. The metals may then facilitate H2 formation on dust particles as well as efficient cooling through metal emission lines (both atomic and molecular), resulting in further fragmentation. In such a case, only a small minority of all halos forms by direct accretion of pristine gas from the IGM, as we have assumed in our simulations.

Alternatively, the final BH may contain a small fraction of the initial mass of its parent clump of gas. Since we did not follow the collapse of each clump down to its Schwarzschild radius, we cannot be certain of the final BH formation efficiency. It is possible that hydrodynamic or radiative feedback from a growing BH seed at the center limits the final mass that it obtains, so that eventually most of the mass of its parent clump is expelled in a wind due to energy release near the center (e.g., Haehnelt, Natarajan, & Rees 1998; Silk & Rees 1998). The correlation between BH mass and the velocity dispersion of its host galaxy as observed in the local universe (Merritt & Ferrarese 2001; Tremaine et al. 2002) would imply an exceedingly low efficiency ($`10^4`$) for the BH formation process in the dwarf galaxies under consideration (Wyithe & Loeb 2002c). However, this phenomenological correlation is measured in metal-rich galaxies where molecular cooling and star formation are abundant, and it may not hold under the unusual physical conditions inside the first generation of metal-free galaxies where the seeds for the present-day population of BHs were produced.

An early formation of SMBHs would lead to a considerable rate of detectable bursts of gravitational radiation from the coalescence of BH binaries at high redshifts (Wyithe & Loeb 2003b). In addition to the commonly discussed origin of binaries in galaxy mergers (a possible example of which has recently been observed in the nearby galaxy NGC 6240 by Komossa et al. 2003), we have demonstrated that BHs may form in binaries to start with, just as stars do.

We are indebted to Lars Hernquist for making available to us a version of TREESPH. We also thank Rennan Barkana and the anonymous referee for valuable comments on the manuscript. This work has been supported in part by NSF grants AST-0071019, AST-0204514. VB thanks the Institute for Advanced Study for its hospitality during the completion of this work. AL acknowledges support from the Institute for Advanced Study at Princeton and the John Simon Guggenheim Memorial Fellowship.


format: text

text:
# Semiclassical quantization of maps with a variable time scale

## I Introduction

It is known for many applications  that classical dynamics near hyperbolic points can be described by a map, $`\widehat{𝒯}`$ in the energy–time $`(E,t)`$ canonical variables $`\widehat{𝒯}(E,t)(E,t)`$. This map is known as the separatrix map, and defines a motion in the vicinity of a separatrix, where a period of the unperturbed (twist) map is an arbitrary function of energy $`T(E)`$. Its explicit form is defined by $`n+2`$ turning points as follows

$$T(E)\{\begin{array}{cc}\mathrm{log}\frac{1}{|E|}\hfill & \text{if }n=1\hfill \\ E^{(n1)/2}\hfill & \text{if }n>1\text{,}\hfill \end{array}$$
(1.1)
where the energy of the separatrix without loss of generality is taken to be zero and it is assumed that $`|E|1`$. A perturbation that is a periodic function of time, for instance, $`ϵ\mathrm{sin}\nu t`$, that is relevant for a variety of applications, may be considered. In this case, $`\widehat{𝒯}`$ is an energy balance equation that describes the energy change over the period $`T(E)`$:

$`E_{n+1}=E_n+ϵ\mathrm{sin}\nu t_n`$ (1.2)
$`t_{n+1}=t_n+T(E_{n+1}),`$ (1.3)
where $`\nu `$ is the frequency of the perturbation of the strength $`ϵ`$, and it is assumed that $`E_n=E(t_n0)`$. The period of the unperturbed nonlinear motion $`T(E_n)`$ describes a wide class of nonlinear systems with variety of applications, including: the celestial mechanics of the perturbed Kepler system , charge particles in a field of a wave packet  leading to a separatrix mesh phenomenon in non-KAM (Kolmogorov-Arnold-Moser) systems  ( with possible applications for atom cooling traps ), and of electron dynamics in superconducting Josephson junctions . Maps like (1.2) are also related to a description of Rydberg atoms in a microwave field by the Kepler map . Also Bloch electrons in external fields are described by a separatrix map  and similar systems.

The map (1.2) can be derived from the Hamiltonian,

$$H=H_0(E)+\stackrel{~}{V}(t,\tau )=H_0(E)+(ϵ/\nu )\mathrm{cos}(\nu t)\delta _{2\pi }(\tau ),$$
(1.4)
where the unperturbed Hamiltonian

$$H_0(E)=\frac{1}{2\pi }^ET(E^{})𝑑E^{}$$
(1.5)
depends only on the energy $`E`$ that is conjugate to the time $`t`$, while $`\tau `$ is the formal time parameter, and $`\delta _{2\pi }(\tau )_{n=\mathrm{}}^{\mathrm{}}\delta (\tau 2\pi n)`$ is the periodic $`\delta `$-function with the period $`2\pi `$. Quantization of the Hamiltonian (1.4) in the framework of the energy-time $`(E,t)`$ canonical pair has been presented in previous studies (see for instance ). The main deficiency of quantization of the Hamiltonian (1.4) is an appearance of the unphysical time for a wave function and this fact has been pointed out earlier, in publications . Classically it is always possible to establish a link between the formal time parameter $`\tau `$ and the real time $`t`$ for any individual trajectory by the solution $`t=t(\tau )`$ (or by canonical transformation). In quantum mechanics (for a wave function), it is no longer possible. Therefore, this approach of quantization in the framework of $`(E,t)`$ variables was used to study time independent characteristics, such as localization . On the other hand, properties which are explicitly time dependent, such as dynamical correlation functions or dynamical fluctuations, cannot be studied by this approach. Moreover, there is no rigorous proof of the validity of these asymptotic solutions. It is also believed that the main deficiency is related to the semiclassically approximate definition of the energy operator as $`\widehat{E}=i\stackrel{~}{h}/t`$, where $`\stackrel{~}{h}`$ is the effective Planck’s constant. There is no physical justification of this semiclassical definition. Therefore, to overcome these shortcomings of the quantum description, it is reasonable to rewrite the system in such a form, where the time parameter appears as the physical time. It is convenient to rewrite the system (1.2) and (1.4) in terms of the action-angle variables $`(I,\theta )`$, related by a canonical transformation $`H(E,t,\tau )(I,\theta ,t)`$, where $`t`$ is the real time. This canonical transformation must be such that under this variable change one transforms $`H_0(E)_0(I)`$ and $`\stackrel{~}{V}(t,\tau )V(\theta ,t)`$. It should be stressed that for kicked systems this change of variables has some specific properties, namely, the new potential $`V(\theta ,t)`$ must be independent of $`I`$ and it must be a function only of the phase $`\theta `$ and the time $`t`$. Otherwise a shift in action at a kick is an ambiguous procedure due to the discontinuity of the action $`I(t=t_n)`$ at the $`n`$-th kick.

In what follows, for the semiclassical quantization, the map (1.2) will be rewritten in the action–angle variables. In this case the semiclassical quantization procedure is standard . Semiclassical quantization of area–preserving maps was subject of earlier studies. In particular quantization of monotonic twist maps , with interpolating flows  was considered. We show first, in Sec. 2, that the map (1.2) can be rewritten in terms of action–angle variables. Semiclassical quantization and an exactly solvable problem are presented in Sec. 3. The results are summarized in Sec. 4.

## II Exact Canonical Transformation $`(E,t)(I,\theta )`$

Hamiltonian equations of motion that produce the map (1.2) have the following formal form

$`(a):`$ $`dE/d\theta =H/t=ϵ\mathrm{sin}\nu t{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}\delta (\theta 2\pi n)`$ (2.1)
$`(b):`$ $`dt/d\theta =H/E=T(E)/2\pi ,`$ (2.2)
where $`H`$ is given by (1.4) and $`\frac{d}{d\theta }`$ is the formal time-derivative ($`\tau \theta `$). Integration over $`\theta `$ in the limits $`(2\pi n0,2\pi n+2\pi 0)`$ gives (1.2). We obtain from the equations (2.1a,b) the Hamiltonian equations for the action–angle variables $`(I,\theta )`$ and the real time parameter $`t`$. Let us start from (2.1b). Inverting this equation, one obtains

$$\frac{d\theta }{dt}=\left[T(E)/2\pi \right]^1=\mathrm{\Omega }(I),$$
(2.3)
where it is supposed that $`T(E)0`$, that means that $`t`$ is a single valued function of $`\theta `$. The frequency $`\mathrm{\Omega }(I)`$ corresponds to some parametric variable change $`E=E(I)`$ along a trajectory, such that

$$dE/dI=\mathrm{\Omega }(I).$$
(2.4)
In this case the energy $`E=_0(I)`$ is an unperturbed Hamiltonian with the new action variable $`I`$, that is defined by (2.4). From (2.1)–(2.4) we perform the following chain of parametric changes

$$\frac{dE}{d\theta }=\frac{dt}{d\theta }\frac{dE}{dt}=\frac{dt}{d\theta }\frac{dE}{dI}\frac{dI}{dt}=\dot{I}.$$
(2.5)
It follows from (2.1a) and (2.5) that

$$\dot{I}=ϵ\mathrm{sin}\nu t\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\delta (\theta 2\pi n)=ϵ\mathrm{sin}\nu t\frac{}{\theta }\mathrm{\Theta }(\theta ),$$
(2.6)
where $`\mathrm{\Theta }(\theta )`$ is a step function of the following form

$$\mathrm{\Theta }(\theta )=\frac{\theta }{2\pi }+\frac{1}{\pi }\underset{l=1}{\overset{\mathrm{}}{}}\frac{1}{l}\mathrm{sin}l\theta ,$$
(2.7)
and it is a result of integration of the periodic $`\delta `$-function written in the form

$$\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\delta (\theta 2\pi n)=\frac{1}{2\pi }+\frac{1}{\pi }\underset{l=1}{\overset{\mathrm{}}{}}\mathrm{cos}l\theta .$$
(2.8)
Equations (2.3) and (2.6) are the Hamilton equations generated by the Hamiltonian

$$=_0(I)2ϵ^{}\mathrm{sin}\nu t\left[\theta +\underset{l=1}{\overset{\mathrm{}}{}}\frac{1}{l}\mathrm{sin}l\theta \right].$$
(2.9)
where $`ϵ^{}=\frac{ϵ}{2\pi }`$.

To restore periodicity of the Hamiltonian in $`\theta `$, the following gauge–like transformation is carried out

$$J=I+\frac{ϵ^{}}{\nu }\mathrm{cos}\nu t.$$
(2.10)
The corresponding Hamiltonian is

$$=_0(J\frac{ϵ^{}}{\nu }\mathrm{cos}\nu t)2ϵ^{}\mathrm{sin}\nu t\underset{l=1}{\overset{\mathrm{}}{}}\frac{1}{l}\mathrm{sin}l\theta _0+V.$$
(2.11)
The equations of motion (2.3) and (2.6) take the form

$`\dot{J}=2ϵ^{}\mathrm{sin}\nu t{\displaystyle \underset{l=1}{\overset{\mathrm{}}{}}}\mathrm{cos}l\theta `$ (2.12)
$`\dot{\theta }=\mathrm{\Omega }(J{\displaystyle \frac{ϵ^{}}{\nu }}\mathrm{cos}\nu t).`$ (2.13)
In this form, the equations were derived from the Hamiltonian $``$, therefore the transformation from $`H`$ to $``$, that was presented here, is a canonical transformation.

## III Semiclassical quantization

### A Floquet operator

The system (2.12),(2.11) can be simply quantized semiclassically with

$`J\widehat{J}=i\stackrel{~}{h}/\theta =\stackrel{~}{h}\widehat{n},`$

where $`\stackrel{~}{h}`$, in our dimensionless units, is a dimensionless parameter, that plays the role of Planck’s constant, while $`\widehat{n}`$ is a quantum numbers operator. In what follows $`\stackrel{~}{h}`$ will be used for the semiclassical quantization. Since the Hamiltonian (2.11) is periodic both in time and in the angle $`\theta `$, the Floquet theory can be used for the analysis. Therefore, a solution of the Schrödinger equation

$$i\stackrel{~}{h}\frac{}{t}\psi (\theta ,t)=\widehat{}(\widehat{n},\theta ,t)\psi (\theta ,t)$$
(3.1)
can be considered due to the Floquet theorem in the following form

$$\psi (\theta ,t)=e^{i\lambda t}e^{i\kappa \theta }\psi _{\lambda ,\kappa }(\theta ,t),$$
(3.2)
where periodicity of $`\widehat{}`$ in both $`\theta `$ and $`t`$ is taken into account, and the functions

$$\psi _{\lambda ,\kappa }\psi _{\lambda ,\kappa }(\theta ,t)=\psi _{\lambda ,\kappa }(\theta +2\pi ,t+2\pi /\nu )$$
(3.3)
are periodic in time with the periods of the perturbation $`2\pi /\nu `$ and periodic in $`\theta `$ with the period $`2\pi `$, while $`0\kappa <1`$ and $`\lambda `$ are a “quasimomentum” and a quasienergy correspondingly. Taking into account that commutation of the Hamiltonian with the exponential $`e^{i\kappa \theta }`$ leads to the shift on $`\kappa `$ for $`\widehat{n}`$,

$$\widehat{}(\widehat{n},\theta ,t)e^{i\kappa \theta }\psi _{\lambda ,\kappa }(\theta ,t)=e^{i\kappa \theta }\widehat{}(\widehat{n}+\kappa ,\theta ,t)\psi _{\lambda ,\kappa }(\theta ,t),$$
(3.4)
we can rewrite the Schrödinger equation (3.1) in the following form

$$\widehat{F}(\kappa )\psi _{\lambda ,\kappa }=\stackrel{~}{h}\lambda \psi _{\lambda ,\kappa },$$
(3.5)
where $`FF(\kappa )`$ is the so–called Floquet operator

$$\widehat{F}(\kappa )=i\stackrel{~}{h}\frac{}{t}+\widehat{}(\widehat{n}+\kappa ,\theta ,t),$$
(3.6)
with $`\lambda `$ and $`\psi _{\lambda ,\kappa }`$ as the eigenvalues and eigenfunctions correspondingly. The solution is cast in the form of the Fourier expansion

$$|\psi _{\lambda ,\kappa }=\underset{n,j}{}\varphi _{n,j}|n|j,$$
(3.7)
moreover the Fourier harmonics

$$|n,j|n|j=\frac{\sqrt{\nu }}{2\pi }e^{in\theta }e^{ij\nu t}$$
(3.8)
are the eigenfunctions of the unperturbed system with $`ϵ=0`$, and $`\varphi _{n,j}=j|n|\psi _{\lambda ,\kappa }`$ are the coefficients of the expansion. Matrix elements of the Floquet operator

$$F_{n,n^{}}^{j,j^{}}=j|n|\widehat{F}|n^{}|j^{}$$
(3.9)
specify the equations for the coefficients $`\varphi _{n,j}`$. First we calculate the matrix elements for $`\widehat{}_0`$. For this purpose we rewrite it formally as

$$\widehat{}_0(\widehat{J}+\kappa \frac{ϵ^{}}{\nu }\mathrm{cos}\nu t)=e^{i\frac{ϵ^{}}{\nu }\mathrm{cos}\nu t(i_\kappa )}\widehat{}_0(\widehat{J}+\kappa ),$$
(3.10)
where $`_\kappa /\kappa `$ and it is assumed that $`_0`$ does not operate on functions of $`\kappa `$. Calculation of the matrix elements (3.9) yields for (3.10)

$$_{0}^{}{}_{n,n^{}}{}^{j,j^{}}=j|e^{i\frac{ϵ^{}}{\nu }\mathrm{cos}\nu t(i_\kappa )}|j^{}n|\widehat{}_0(\widehat{J}+\kappa )|n^{}=\underset{m}{}i^mJ_m(i\frac{ϵ^{}}{\nu }_\kappa )\delta _{j^{},j+m}_0(n^{}+\kappa )\delta _{n^{},n},$$
(3.11)
where $`J_m(x)`$ is the Bessel function. Correspondingly, the matrix elements for the additive part of the perturbation in (2.11) are

$$V_{n,n^{}}^{j,j^{}}=\{\begin{array}{cc}\hfill \frac{1}{2}\frac{ϵ^{}}{nn^{}}\left[\delta _{j^{},j+1}\delta _{j^{},j1}\right]:& nn^{}\hfill \\ \hfill 0:& n=n^{}\hfill \end{array}.$$
(3.12)
These matrix elements of (3.11) and (3.12) lead to the following equation for the expansion coefficients

$$\underset{m}{}i^mJ_m(i\frac{ϵ^{}}{\nu }_\kappa )_0(n+\kappa )\varphi _{n,j+m}+\frac{ϵ^{}}{2}\underset{n^{}n}{}\frac{1}{nn^{}}[\varphi _{n^{},j+1}\varphi _{n^{},j1}]=\stackrel{~}{h}(\lambda +\nu j)\varphi _{n,j}.$$
(3.13)
It is similar to a quasi-1D Anderson-like chain for the cases when $`_0(n+\kappa )`$ corresponds to a random potential. It is also the generalization of a specific application in . The classical motion is chaotic in many types of $`T(E)`$ in (1.2). Therefore, the obtained equation establishes a relation between a wide class of energy balance equations of the form (1.2), (1.4) and dynamical localization of the classical chaotic diffusion . In the general case this equation is useful for numerical studies. It is analytically tractable in perturbation theory for the parameter $`\frac{ϵ}{\nu }1`$. Therefore, it is instructive to present an exactly solvable example where (1.5) is linear, with a constant period $`T(E)=2\pi /\omega `$. In this case a solution can be obtained in an explicit form analytically.

### B An exact model

An exact solution can be obtained for the harmonic oscillator $`_0(I)=\omega I`$ with constant frequency $`\mathrm{\Omega }(I)=\omega `$. The linear driven system described by the Hamiltonian $`(I,\theta ,t)`$ is integrable. The equation (3.13) takes the following simple form

$$\stackrel{~}{h}[\omega n\nu j\lambda +\omega \kappa ]\varphi _{n,j}\frac{ϵ^{}\omega }{2\nu }[\varphi _{n,j+1}+\varphi _{n,j1}]+\frac{ϵ^{}}{2}\underset{n^{}n}{}\frac{1}{nn^{}}[\varphi _{n^{},j+1}\varphi _{n^{},j1}]=0.$$
(3.14)
For the specific case of the harmonic oscillator $`\kappa =0`$. For other linear models, for example models of an appropriate band structure in solids, $`\kappa `$ does not have to vanish. A solution of equation (3.14) is cast in the form of a sinc function by the following substitution

$$\varphi _{n,j}=Z_j\text{sinc}[\pi (n\lambda /\omega +\kappa )],$$
(3.15)
where the sinc function is $`\text{sinc}x=\frac{\mathrm{sin}x}{x}`$ and $`Z_j`$ has to be determined. In this case the quasienergy spectrum turns out to be

$$\lambda =\omega \kappa +\omega n.$$
(3.16)
Inserting (3.6) in (3.14) and taking into account that

$$\underset{n^{}n}{}\frac{1}{nn^{}}\text{sinc}\pi (n^{}n)=0,$$
(3.17)
since $`\mathrm{sin}\pi (n^{}n)=0`$, one obtains from (3.14) the following relation for the $`Z_j`$

$$2jZ_j+z[Z_{j+1}+Z_{j1}]=0,$$
(3.18)
where $`z=ϵ^{}\omega /\stackrel{~}{h}\nu ^2`$. The solution is

$$Z_jZ_j(z)=(1)^jJ_j(z).$$
(3.19)
The result (3.16) and (3.19) is only a specific solution. The general solution is found by replacing in (3.14) $`j`$ by $`jj_0`$ and $`\lambda `$ by $`\lambda \nu j_0`$. The resulting solution of the eigenvalue problem is

$$\lambda _{j_0}=\nu j_0+\omega \kappa +\omega n$$
(3.20)
and the corresponding eigenfunction is

$$Z_j=(1)^{jj_0}J_{jj_0}(z).$$
(3.21)
Floquet theory implies that only $`\lambda _{j_0}`$ in an interval of size $`\nu `$, say $`[0,\nu )`$ should be used (see (3.2)). In particular for the harmonic oscillator

$$\lambda _{j_0}=\nu j_0+\omega n,$$
(3.22)
if it is in the interval $`[0,\nu )`$, as expected.

## IV Summary

Quantization of energy balance equations, in the form of the map (1.2), is presented. This procedure of quantization consists of the two steps. The first one is a transformation of the Hamiltonian (1.4) to the form (2.11) by the canonical transformation $`(E,t)(I,\theta )`$. The second stage is semiclassical quantization. The Hamiltonian formulation of the problem allows semiclassical quantization of action–angle variables in the framework of the standard procedure: $`I\widehat{I}=i\stackrel{~}{h}/\theta \stackrel{~}{h}\widehat{n}`$, where $`\stackrel{~}{h}`$ is a dimensionless semiclassical parameter such that the classical limit is $`\stackrel{~}{h}0`$. It should be noted that the semiclassical approximation requires that the width of the potential $`V`$ is larger than the de-Broglie wavelength. Therefore, it is necessary to truncate the Fourier expansion in (2.7), corresponding to replacement of the $`\delta `$-function by a function of finite width. The results obtained in Sec. 3.2 in the semiclassical approximation are correct only up to the leading correction in $`\stackrel{~}{h}`$ to the classical ones . The advantage of the present analysis is that the transition from (1.2) to (2.11) is exact, and wave functions depend on physical time, while the main deficiency of straightforward quantization of (1.2) is an appearance of an unphysical time parameter .

## Acknowledgments

We thank V. Afraimovich, I. Guarneri, A. Neishtadt and V. Rom-Kedar for fruitful discussions. We also thank U. Smilansky for drawing our attention to the results presented in . This research is supported by U.S. Department of Navy Grant No N00014-02-1-0056, the US–Israel Binational Science Foundation (BSF) and by the Minerva Center of Nonlinear Physics of Complex Systems.