text:
# Two-color stabilization of atomic hydrogen in circularly polarized laser fields

## I Introduction

For sufficiently high laser frequency and intensity, the ionization probability of an atom decreases by further increasing the laser intensity. This, at first sight, counterintuitive effect is called dynamic or adiabatic stabilization, depending on the pulse length pont89 ; pontgav90 ; pontwalet ; gavrevI ; eberlykul ; mittlemanbook ; faisal ; gavrilarevII . So far, it has been observed in numerical simulations only (see kulander  and, for very recent results, dondera ; choi ). Its theoretical explanation is most elucidating in the framework of high-frequency Floquet theory gavrevI : the time-averaged Kramers-Henneberger potential, as it is seen by a free electron moving in the laser field, provides bound states. The electronic probability density of these states is spread over spatial regions with dimensions of the order of the classical excursion amplitude which may be tens of atomic units or more. The laser frequency $`\omega `$, required for dynamic stabilization where the electron starts from the field-free ground state with binding energy $`_0`$, has to exceed the field-dressed binding energy, $`\mathrm{}\omega >|_0^{}|`$. Since such high frequencies are not yet accessible to nowadays intense laser systems the experimental verification of dynamic stabilization is still to come. However, with the intense, short-wavelength, and coherent light sources under construction worldwide (e.g., at DESY in Hamburg) ground state stabilization should become experimentally accessible soon. Signatures of dynamic stabilization of Rydberg atoms, where $`|_0^{}|`$ is smaller, have been measured druten .

Dynamic stabilization has been studied extensively during the last decade. Ionization rates and energy shifts were obtained numerically by Floquet theory doerr  or approximate, analytical high-frequency Floquet theory pontgav90 ; gavrevI . Full numerical ab initio solutions of the time-dependent, single-electron Schrödinger equation were obtained for low-dimensional model systems (see, e.g., su ; volkova ; patelI ; patelII ; kwon ; ryab ) and, after the early work in kulander  very recently, in three spatial dimensions for linear polarization dondera  and circularly polarized laser pulses choi . Although the Floquet ionization rates decrease monotonically for laser intensity $`I\mathrm{}`$ the ionization probability in full ab initio solutions of the TDSE for finite laser pulses assumes a minimum and increases again thereafter. In the limit of laser intensity $`I\mathrm{}`$ but otherwise fixed pulse parameters the survival probability of an atom will thus approach zero. Moreover, it has been found that the magnetic field of the laser—usually neglected because of the dipole approximation—leads to increased ionization kylstra . This is due to the ponderomotive force which pushes the electrons in propagation direction. In multi-electron atoms electron correlation also hinders stabilization (see bauer , and references therein). Two-color interference effects on stabilization in linear polarization was studied in cheng ; potvl .

In this paper the dynamic stabilization of H in circularly polarized single- and two-color fields is investigated. The article is organized as follows: in Sec. II the three-dimensional, time-dependent Schrödinger equation to be solved numerically is introduced, and the form of the applied two-color fields is explained. In Sec. III the results concerning the single-color, circular polarization case are presented. Besides snapshots of the probability density and photoelectron spectra, particular emphasis is put on the comparison of ionization rates and energy shifts with results from Floquet theory. In Sec. IV the main two-color results are presented. Probability density structures are observed which display well separated maxima, the number of which is related to the frequency ratio of the two fields. The structure moves as a whole in the laser field and can be explained with the help of two-color, time-averaged Kramers-Henneberger potentials. Survival probabilities after the two-color pulses are presented as a function of the two laser’s peak field strengths. The two-color survival probabilities and above-threshold ionization spectra are compared with the corresponding single-color results where the second field was absent. It is good news that adding a second color is possible without a drastic inhibition of stabilization even if the second laser’s intensity lies in the “death valley” of high ionization. Finally, we conclude in Sec. V.

## II Theory

The time-dependent Schrödinger equation (TDSE) for the electron of atomic hydrogen in the electric field of a laser

$$\text{i}\frac{}{t}\mathrm{\Psi }(𝒓,t)=\left(\frac{1}{2}[\text{i}\mathbf{}+𝑨(t)]^2\frac{1}{r}\right)\mathrm{\Psi }(𝒓,t)$$
(1)
is solved numerically in three spatial dimensions (atomic units are used unless noted otherwise). The intense laser field can be treated classically. Since the laser wavelength is large compared to all other relevant length scales and the electron dynamics is non-relativistic for our laser parameters, the dipole approximation can be applied so that the vector potential $`𝑨(t)`$ does not depend on spatial coordinates. Hence, the magnetic component of the laser field is neglected. If we restrict ourselves to study coplanar two-color fields (of, however, arbitrary ellipticity) we may choose the coordinate system in such a way that the $`z`$-component of the electric field $`𝑬(t)=_t𝑨(t)`$ vanishes: $`𝑬(t)=𝑬^{(1)}(t)+𝑬^{(2)}(t)=[E_x^{(1)}(t)+E_x^{(2)}(t)]𝒆_x+[E_y^{(1)}(t)+E_y^{(2)}(t)]𝒆_y`$.

The wave function $`\mathrm{\Psi }(𝒓,t)`$ is expanded in spherical harmonics,

$$\mathrm{\Psi }(r,\vartheta ,\phi ,t)=\frac{1}{r}\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{\Phi }_\mathrm{}m(r,t)Y_\mathrm{}m(\vartheta ,\phi ),$$

which leads to a set of coupled TDSEs for the radial wave functions $`\mathrm{\Phi }_\mathrm{}m(r,t)`$,

$`\text{i}{\displaystyle \frac{}{t}}\mathrm{\Phi }_\mathrm{}m(r,t)`$ $`=`$ $`\left({\displaystyle \frac{1}{2}}{\displaystyle \frac{^2}{r^2}}+V_{\mathrm{}}^{\text{eff}}(r)\right)\mathrm{\Phi }_\mathrm{}m(r,t)`$
$`+\text{i}\sqrt{{\displaystyle \frac{2\pi }{3}}}{\displaystyle \underset{\mathrm{}^{}m^{}}{}}\{\stackrel{~}{A}^{}(t)\mathrm{}m|11|\mathrm{}^{}m^{}{\displaystyle \frac{}{r}}\stackrel{~}{A}(t)\mathrm{}m|11|\mathrm{}^{}m^{}{\displaystyle \frac{}{r}}`$
$`{\displaystyle \frac{\stackrel{~}{A}^{}(t)}{r}}\mathrm{}m|11|\mathrm{}^{}m^{}(1+m^{})+{\displaystyle \frac{\stackrel{~}{A}(t)}{r}}\mathrm{}m|11|\mathrm{}^{}m^{}(1m^{})`$
$`{\displaystyle \frac{\stackrel{~}{A}^{}(t)}{r}}\mathrm{}m|10|\mathrm{}^{}m^{}+1N_\mathrm{}^{}m^{}^++{\displaystyle \frac{\stackrel{~}{A}(t)}{r}}\mathrm{}m|10|\mathrm{}^{}m^{}1N_\mathrm{}^{}m^{}^{}\}\mathrm{\Phi }_\mathrm{}^{}m^{}(r,t)`$
where $`V_{\mathrm{}}^{\text{eff}}(r)=\frac{1}{r}+\frac{\mathrm{}(\mathrm{}+1)}{2r^2}`$ is the effective atomic potential including the centrifugal barrier, $`\stackrel{~}{A}(t)=A_x(t)+\text{i}A_y(t)`$, $`N_\mathrm{}m^\pm =\sqrt{(\mathrm{}m)(\mathrm{}\pm m+1)/2}`$, and the purely time-dependent $`𝑨^2`$-term has been transformed away. The matrix elements $`\mathrm{}m|LM|\mathrm{}^{}m^{}:=\text{d}(\mathrm{cos}\vartheta )\text{d}\phi Y_\mathrm{}m^{}(\vartheta ,\phi )Y_{LM}(\vartheta ,\phi )Y_\mathrm{}^{}m^{}(\vartheta ,\phi )`$ may be expressed in terms of Clebsch-Gordan coefficients.

It is worth noticing the coupling of different $`\mathrm{}`$ and $`m`$ quantum numbers in the TDSE (II) due to the laser field: there is a $`\mathrm{}m\mathrm{}\pm 1,m\pm 1`$ coupling as well as a $`\mathrm{}m\mathrm{}1,m\pm 1`$ coupling (where compatible with the condition $`\mathrm{}m\mathrm{}`$). From the knowledge of these couplings one can immediately deduce non-perturbative selection rules, e.g., if the initial wave function is the 1s groundstate, states with ($`\mathrm{}=1`$, $`m=0`$), with ($`\mathrm{}=2`$, $`m=\pm 1`$), with ($`\mathrm{}=3`$, $`m=\pm 2`$) etc. are ‘unreachable’ and thus remain unoccupied.

In the actual implementation, the radial coordinate $`r`$ was discretized directly in position space. The details of the implementation may be found in bauerhabil  and is similar to what has been described in mullermethod  for the simpler case of linear polarization.

The radius of our numerical grid, including the region where the non-vanishing imaginary potential damps away probability density approaching the grid boundary, was between 240 a.u. and 450 a.u. The maximum $`\mathrm{}`$ quantum number was typically $`\mathrm{}_{\mathrm{max}}=30`$. For testing the convergence of our results runs up to $`\mathrm{}_{\mathrm{max}}=50`$ were also performed—with no significant changes in the observables of interest. For obtaining the results presented in this paper, the run times of our Schrödinger-solver were always less than 18 hours on 1 GHz PCs.

We focus on two-color fields where both fields are circularly polarized. We assume a $`\mathrm{sin}^2`$-shaped up- and down-ramp of the $`i`$th field, $`i=1,2`$ over $`N_1^{(i)}=N_{23}^{(i)}`$ cycles and a constant part of $`N_{12}^{(i)}`$ cycles duration. In the following we will refer to such a pulse as a $`(N_1^{(i)},N_{12}^{(i)},N_{23}^{(i)})`$-pulse. For the $`i`$th electric field in, say, $`x`$-direction

$$E_x^{(i)}(t)=\{\begin{array}{ccc}E_{\text{ud}}^{(i)}(t)\hfill & \text{for}& 0t<T_1,\hfill \\ \widehat{E}^{(i)}\mathrm{cos}\omega _it\hfill & \text{for}& T_1t<T_2,\hfill \\ E_{\text{ud}}^{(i)}(tT_2+T_1)\hfill & \text{for}& T_2t<T_3,\hfill \\ 0\hfill & \text{otherwise}& \end{array}$$
(3)
is set where $`E_{\text{ud}}^{(i)}(t)=\widehat{E}^{(i)}\mathrm{sin}^2(\pi t/2T_1)\mathrm{cos}\omega _it`$. The absolute ramping time as well as the duration of the constant part of the laser pulse are chosen the same for both colors so that the times $`T_1=2\pi N_1^{(i)}/\omega _i`$, $`T_2=T_1+2\pi N_{12}^{(i)}/\omega _i`$, etc. are independent of $`i`$. Since we are interested in opposite circular polarizations

$$E_y^{(1)}(t)=E_x^{(1)}(t3\pi /2\omega _1),E_y^{(2)}(t)=E_x^{(2)}(t\pi /2\omega _2)$$
(4)
are taken as the $`y`$-components of the two colors. The corresponding vector potential can be calculated from $`𝑨(t)=_0^t\text{d}t^{}𝑬(t^{})`$. In such two-color pulses, as long as $`N_1^{(i)}=N_{23}^{(i)}`$, $`N_{12}^{(i)}`$ are integer numbers, a free electron (initially at rest) returns to its starting point at the end of the laser pulse—without residual drift momentum.

## III The circular single-color case revisited

Before turning to the two-color study let us analyze the simpler single-color case. For the latter, results from ab initio solutions of the three-dimensional TDSE were obtained only recently dondera  whereas two-dimensional model solutions were presented earlier patelI ; patelII ; kwon ; ryab . Ionization rates as well as energy shifts for hydrogen in the stabilization regime (with dipole approximation) are available for several years doerr  so that a comparison with our finite pulse results is possible and will be presented in this paper. Since all vector potential amplitudes $`\widehat{A}`$ considered in this work are $`10`$ a.u. the dipole approximation is applicable, as has been shown in kylstra .

### III.1 Survival probability vs. time

In order to calculate the probability for the hydrogen atom to “survive” non-ionized till the end of the pulse the probability density is integrated over a spherical volume of radius $`R_b=60`$ atomic units. The radius $`R_b`$ was chosen several times greater than the biggest excursion amplitude $`\widehat{\alpha }=10`$ a.u. considered in this paper. Since the free part of the wavefunction vanishes inside that sphere at sufficiently large times after the laser pulse has passed by, the amount of probability density inside the sphere leads then with high accuracy to the same probability one would obtain by projecting on all bound states.

The results are presented in Fig. 1. In the upper panel (a) the results for $`\omega =1`$ are presented, in the lower panel (b) those for $`\omega =2`$. All pulses were up- and down-ramped over 4 cycles. The constant part of the pulse lasted 16 cycles for $`\omega =1`$ and 32 cycles for $`\omega =2`$. By Fig. 1 the very existence of dynamic stabilization is proven: the survival probability at the end of the pulse increases with increasing peak field strength (attached to the curves) up to a maximum value. Increasing the peak field strength further yields a decreasing survival probability (not shown in Fig. 1). The highest survival probability in plot (b) amounts to $`47\%`$. It should be noted, however, that the survival probability is rather sensitive to the laser pulse shape and duration. This was also found in the case of linear polarization and is discussed at length in Ref. dondera .

Whenever the slope of the probability curves vs. time is linear on the logarithmic scale an ionization rate can be determined. While this is straightforward for the curves in the $`\omega =2`$ panel it turns out to be more difficult for the $`\omega =1`$-results in (a). Longer pulses of constant intensity are then needed to obtain unambiguous rates.

### III.2 Snapshots of the probability density

For sufficiently large excursions $`\widehat{\alpha }`$ and adiabatic turn-ons one expects the system to occupy essentially the ground state in the time-averaged Kramers-Henneberger (KH) potential. This corresponds to lowest-order high-frequency Floquet theory gavrevI . The KH potential is the ionic potential a free electron, just oscillating in the laser field, would “see”, i.e., $`V_{\text{KH}}(𝒓)=V(𝒓+𝜶(t))`$ where $`𝜶(t)`$ is the trajectory of a free electron in the field

$$𝜶(t)=_0^t\text{d}t^{}𝑨(t^{}).$$
(5)
In the case of a single-color, circularly polarized laser the KH potential, averaged over one laser cycle, has azimuthal symmetry with respect to the laser propagation direction (the $`z`$-axis in our case) and, for sufficiently big $`\widehat{\alpha }`$, its minimum along a circle of approximate radius $`\widehat{\alpha }`$ in a plane perpendicular to it (the $`xy`$-plane in our case). Consequently, the KH ground state probability density has its maximum along that circle. By transforming back to the lab frame, where the ionic potential is assumed centered at the origin, this ring-shaped, stabilized probability density orbits with its center along the free electron trajectory which itself is a circle of radius $`\widehat{\alpha }`$. What results is a “hula-hoop”-like dynamics choi  of a probability density ring about the ion which assumes the role of the “hula-hoop”-dancer’s hip.

In Fig. 2 we present snaphots of the probability density in the $`xy`$-plane at $`z=0`$ during the 9th laser cycle of the (2,8,2)-pulse with $`\omega =1`$. The peak field strength was $`\widehat{E}=10`$, yielding an excursion radius $`\widehat{\alpha }=10`$. Although a probability density ring, describing the previously explained “hula-hoop”-dynamics, is clearly visible, the probability density is not uniformly distributed along that ring but has rather a sickle-like shape. Hence, the system is obviously not in the KH ground state but in a superposition of states with different azimuthal quantum numbers. This leads to Rabi floppings on time scales that are large compared to the laser period. Note that there is always a probability density maximum at the position of the ion. For lower values of $`\widehat{\alpha }`$ we observed more complex situations where no clear ring structure but two or more wave packets swirl about the ion.

### III.3 Ionization rates

For the single-color case, accurate ionization rates were numerically obtained many years ago by means of non-perturbative Floquet theory doerr . Those rates are accurate in the sense that if the system somehow manages to occupy a single Floquet state it will decay precisely with the corresponding Floquet rate. Usually an adiabatic turn-on is needed to transfer the system from the field-free groundstate to a single Floquet state. The problem is that an adiabatic ramping of the laser pulse counteracts stabilization because the system has to spend too much time in the so-called “death valley” of medium laser intensities where ionization is strong, and thus nothing would survive to be stabilized. Hence, short turn-on times are desirable. However, after such rapid turn-ons the system might be “shaken-up” into a superposition of many Floquet states. Under such conditions only the full solution of the TDSE yields reliable ionization probabilities and, if existing, ionization rates.

In Fig. 3 our ionization rates for $`\omega =1`$ and $`\omega =2`$ (determined from plots like those in Fig. 1) are compared with Floquet results doerr  and the analytical results derived by Pont and Gavrila pontgav90 ; gavrevI . The Pont-Gavrila rate for circular polarization

$$\mathrm{\Gamma }_{\text{PG}}\frac{0.223}{\widehat{\alpha }^4\omega ^2}$$
(6)
was derived from high-frequency Floquet theory combined with the first-order Born approximation. In order to meet the plane wave criterion of the Born approximation the excursion amplitude $`\widehat{\alpha }`$ must not be too small. For high frequencies and $`\widehat{\alpha }1`$ the Pont-Gavrila rate should merge with the full Floquet result. However, in doerr  Floquet rates were calculated only for moderate values of $`\widehat{\alpha }`$ \[note, that, by definition, the excursion amplitude $`\widehat{\alpha }`$ for circular polarization in doerr  differs by a factor $`\sqrt{2}`$ with Pont’s and Gavrila’s, and ours, i.e., $`\mathrm{\Gamma }_{\text{PG}}=\mathrm{\Gamma }_{\text{Fl}}/\sqrt{2}`$\]. For $`\omega =2`$ the Pont-Gavrila rate comes close to the full Floquet result already for $`\widehat{\alpha }=1`$. For increasing $`\widehat{\alpha }`$ the full Floquet result oscillates around the Gavrila rate. In the lower frequency case $`\omega =1`$ the Gavrila rate overestimates ionization. Our numerical rates were obtained from simulations of (4,16,4)-pulses for $`\omega =1`$ and (4,32,4)-pulses for $`\omega =2`$. The results for the highest $`\widehat{\alpha }`$ values shown were difficult to determine because the ionization probability curves vs. time displayed no clear exponential decrease. In order to determine unambiguous ionization rates also for higher $`\widehat{\alpha }`$ a more adiabatic up-ramp of the laser field and a longer pulse duration would be required. However, since it is not the goal of this paper to determine rates which are probably easier to obtain via Floquet theory, we do not proceed further in this direction. Our numerical results agree well with the full Floquet rates, where available. For $`\omega =1`$ and $`\widehat{\alpha }3.5`$ the Gavrila rate agrees quite well with our numerical data whereas for $`\omega =2`$ it underestimates ionization.

### III.4 Above-threshold ionization (ATI)

In Fig. 4 photoelectron spectra after (4,8,4)-pulses (for $`\omega =1`$) and (8,16,8)-pulses (for $`\omega =2`$) of different peak field strength are presented. They were calculated from the wave function immediately after the pulse with the window-operator method proposed in windowop . The three lowest-order ATI peaks $`n=1,2,3`$, separated by $`\mathrm{}\omega `$, are clearly visible although they broaden with increasing laser intensity. From lowest order perturbation theory (LOPT) an exponential decrease of the peak heights $`\widehat{E}^{2n}`$ with increasing $`n`$ is expected. Within high-frequency Floquet theory, for $`\widehat{\alpha }1`$, $`(_0^{}+n\omega )\widehat{\alpha }1`$ the contributions of higher orders $`n`$ do not vanish exponentially but only $`n^2`$ pontgav90 . By virtue of the first two peaks in the highest-$`\widehat{E}`$ spectra of Fig. 4 this trend of increasing importance of higher order peaks is confirmed \[note that in the numerically obtained wavefunction after the pulse the ATI peak $`n=3`$ might be slightly suppressed because of the absorbing boundaries\].

With increasing population of excited states the left shoulder of the ATI peaks becomes more pronounced. For certain laser intensities a substructure of smaller peaks appears. Since these peaks are separated by the energy difference of $`n=2`$ and $`n=3`$-states, that is $`\mathrm{\Delta }=(4^19^1)/20.07`$, this substructure is likely connected to the population of these states. However, we did not perform yet a detailed study of this substructure. Here we want to focus on the ATI peak positions. As can be easily inferred from Fig. 4 the ATI peaks move with increasing laser intensity to higher photoelectron energies. This is contrary to what happens for low frequencies and high intensities. In general, the positions of the ATI peaks follow

$$_{\text{low}}^{(n)}=_0^{}\mathrm{\Delta }_{\text{cont}}+n\omega .$$
(7)
Here, $`_0^{}=_0+\mathrm{\Delta }`$ is the field-dressed ground state energy which is ac Stark-shifted by $`\mathrm{\Delta }`$ compared to the unperturbed ground state energy $`_0`$ ($`=0.5`$ for H), and $`\mathrm{\Delta }_{\text{cont}}`$ is the energy up-shift of the continuum threshold. $`n`$ is the number of absorbed photons $`n_{\mathrm{min}}`$ where $`n_{\mathrm{min}}`$ is the smallest integer that gives $`_{\text{low}}^{(n_{\mathrm{min}})}>0`$. For low frequencies one has $`\mathrm{\Delta }_{\text{cont}}=U_p=f\frac{\widehat{E}^2}{2\omega ^2}`$ where $`U_p`$ is the ponderomotive potential, i.e., the time-averaged quiver energy of a free electron in the laser pulse of peak field strength $`\widehat{E}`$. The pre-factor $`f`$ is (for our definition of $`\widehat{E}`$) $`1`$ for circular and $`1/2`$ for linear polarization. Since in the case of long wavelengths the up-shift of the continuum threshold by $`U_p`$ is bigger than the shift $`\mathrm{\Delta }`$ of the ground state, the ATI peaks move towards lower energies with increasing laser intensities. The “drowning” of a certain peak no. $`n`$ below $`=0`$ is the celebrated “channel closing” in ATI physics.

We shall see that in the stabilization regime, where $`\omega >|_0^{}|`$ must hold, the shift of ATI peaks is dominated by the ground state up-shift $`\mathrm{\Delta }`$ so that

$$_{\text{high}}^{(n)}_0^{}+n\omega .$$
(8)
Accurate values of the field dressed energies $`_0^{}`$ in the stabilization regime can be obtained from Floquet theory doerr . Solving the set of coupled Floquet equations yields complex energies the real part of which is the field-dressed energy. The complex part is half the ionization rate, $`_{0,\text{Fl}}=_0^{}\text{i}\mathrm{\Gamma }_{\text{Fl}}/2`$. From the high-frequency Floquet theory point-of-view $`_0^{}`$ is the groundstate energy of the time-averaged KH potential. In Ref. pontgav90  Pont and Gavrila gave KH ground state energies for H in a circularly polarized laser field. Because $`|_0^{}|<|_0|`$ the ATI peaks move to higher energies with increasing laser intensity—contrary to the high-intensity, low-frequency-case (7).

In Fig. 5 our numerically determined positions for the $`n=1`$ ATI peak are compared with Eq. 8 where the value for $`_0^{}`$ has been determined (i) by the full Floquet solution doerr , (ii) as the time-averaged KH potential ground state energy pontgav90 , and (iii) through Pont’s approximate formula $`_0^{}=(\mathrm{ln}\widehat{\alpha }+2.654284)/(2\pi \widehat{\alpha })`$ pont89 . The agreement may be considered acceptable although not excellent. The differences are likely due to the population of more than a single Floquet state. Our numerical result for $`\omega =1`$ and $`\widehat{\alpha }=5`$ lies above the analytical curves because the substructure in the ATI spectrum (cf. Fig. 4) also alters the absolute maximum of an ATI peak so that its determination becomes ambiguous.

## IV Two-color case with opposite circular polarizations

### IV.1 Time-averaged KH potentials

It is useful for the analysis of our two-color results to get familiar with the free electron trajectories in such fields and the corresponding time-averaged KH potentials. In the two-color case the KH potentials were averaged over the low-frequency laser period $`2\pi /\omega _1`$. Let us consider the two laser frequencies $`\omega _1=1`$ and $`\omega _2=2`$. If the second field is absent the free electron trajectory $`𝜶(t)`$ describes a circle in the $`xy`$-plane of radius $`\widehat{\alpha }_0=\widehat{E}^{(1)}/\omega _1^2`$. If the second field amplitude $`\widehat{E}^{(2)}`$ equals $`\widehat{E}^{(1)}`$, i.e., both laser fields have the same intensity, the trajectory resembles an equilateral triangle with rounded edges, and $`\mathrm{max}|\widehat{\alpha }|=5\widehat{\alpha }_0/4`$. If, instead, the vector potential amplitudes of both fields are equal, $`\widehat{E}^{(1)}/\omega _1=\widehat{E}^{(2)}/\omega _2=`$, the electron follows a convex triangle with sharp edges, and $`\mathrm{max}|\widehat{\alpha }|=3\widehat{\alpha }_0/2`$. Finally, equal excursion amplitudes, i.e., $`\widehat{E}^{(1)}/\omega _1^2=\widehat{E}^{(2)}/\omega _2^2`$, yields a rosette-like trajectory, and $`\mathrm{max}|\widehat{\alpha }|=2\widehat{\alpha }_0`$. The time-averaged KH potentials look accordingly. They are presented in Fig. 6. In the two-color cases (b)–(d) the three-fold symmetry leads to three potential minima near the turning points. If the stabilized electron occupies the KH ground state the probability density will show “trichotomy” instead of the well-known “dichotomy” in the case of a single-color linearly polarized laser where two probability peaks, separated by $`2\widehat{\alpha }`$ were observed kulander . In general, when $`\omega _2=n\omega _1`$ the trajectories display $`(n+1)`$-fold symmetry.

Provided the peak field strengths are properly chosen, for $`n=3,4,\mathrm{}`$ the $`n+1`$ minima of the time-averaged KH potentials will lead to “quatrochotomy,” “pentachotomy,” etc. The electron dynamics in the lab frame should then resemble “hula-hoop” with a triangle, a square, a pentagon, and so forth. We would like to remark that the discrete symmetry of the trajectory and the time-averaged KH potential is closely related to selection rules for harmonic generation in two-color fields (see cecch , and references therein). If the two fields with frequencies $`\omega _1`$ and $`\omega _2=n\omega _1`$ are equally polarized instead of oppositely, the KH potentials have ($`n1`$)-fold symmetry. Hence, to obtain “trichotomy” $`\omega _2=4\omega _1`$ has to be chosen in this case. If the two frequencies are not commensurable the free electron trajectories are not closed, and the two-color, time-average KH potential is ill-defined. However, we observed stabilization also in this case although no distinct probability density peaks were obtained.

### IV.2 Snapshots of the probability density

In Fig. 7 snapshots of the trichotomous “hula-hoop”-dynamics for $`\omega _1=1`$, $`\omega _2=2`$ during one laser cycle (with respect to the $`\omega _1=1`$-field) is presented. The two laser intensities were chosen equal so that the free-electron trajectory is a triangle with rounded edges. This triangle is drawn in black in each of the plots of Fig. 7. As expected, the center of the trichotomous probability density (aligned along the time-averaged KH potential) moves clockwise along the black triangle. Note that the stabilized structure moves as a whole and does not rotate about its center. In the comoving frame, i.e., the KH frame, the ion swirls along the triangle, forming another probability density maximum.

In Fig. 8 three other examples for stabilized probability density are presented. In panel (a), where the excursion amplitudes for both fields were chosen equal, no three distinct peaks are visible because $`\widehat{\alpha }=5`$ is too small. In (b) a clear example of equal intensity quatrochotomy (i.e., $`\omega _2=3`$) is shown and in (c) a case of pentachotomy ($`\omega _2=4`$) where both fields had equal vector potential amplitudes. The density of the $`n+1`$ peaks varies during the course of the laser pulse. This implies that several Floquet states are involved.

### IV.3 Survival probability vs. time

So far the two-color survival probability was not addressed at all. If, for instance, $`\widehat{E}^{(2)}`$ is chosen so that the $`\omega _2=2`$-field alone yields optimal stabilization, and a $`\omega _1=1`$-field with $`\widehat{E}^{(1)}`$ in the “death valley” is added—what is the survival probability? In the worst case ionization would be dominated by the “death valley”-field, making the two-color stabilization not very attractive. Fortunately, it is not the field which would lead to highest ionization, if applied alone, that dominates two-color dynamic stabilization. In Fig. 9 the temporal evolution of the survival probability in the two-color case is compared with the two corresponding single-color results. From panel (a) it is clearly seen that if the two laser intensities are tuned to yield equal excursion amplitudes $`\widehat{\alpha }`$ the two-color survival probability is closer to the single-color high-frequency result with relatively low ionization probability. The low-frequency field alone has a higher ionization probability. In the equal intensity plot (b) it is the other way round: the two-color result follows closely the low-frequency curve whereas the high-frequency field alone yields higher ionization because $`\widehat{E}^{(2)}`$ is still close to the “death valley.” Since in plot (a) all curves cross each other it is apparent that the final survival probability is sensitive to the pulse length. However, it should be stressed that, fortunately, the two-color survival probability is not determined by the frequency which yields highest ionization. This means that the two-color “probability density shaping,” as discussed in the previous subsection, can be performed without drastic reduction of the survival probability.

In Fig. 10 the survival probabilities at the end of (2,8,2)- and (4,16,4)-pulses for $`\omega _1=1`$ and $`\omega _2=2`$, respectively, are plotted vs. the field amplitudes $`\widehat{E}^{(1)}`$ and $`\widehat{E}^{(2)}`$. Because of the probability density representing free but slow electrons it takes a long time until the probability density integrated over the sphere of radius $`R_b=60`$ assumes a constant value. Hence, to obtain Fig. 10 we have chosen to measure the survival probability as the maximum probability to find the electron in a smaller sphere of radius $`R_s=7.5`$ within a time period of 85 a.u. after the pulse. The maximum probability has to be taken because of oscillations due to excitation of the 2s and the 3s state. Note that this method for measuring the survival probability slightly overestimates ionization because probability density which represents slow recombining electrons, slowly spirals back to the ion and does not enter the small sphere within the observation time. However, here we are not interested in absolute numbers (which depend sensitively on the pulse duration and shape anyway). From Fig. 10 one infers that the highest survival probability is obtained, as expected, in the single-color high-frequency case where $`\widehat{E}^{(1)}=0`$. In general, the overall behavior is an increasing stabilization probability up to a maximum value, followed by increasing ionization if the laser intensities are increased further. This is different from high-frequency Floquet theory where the ionization rate monotonically decreases with increasing intensity because the “shake-off” during the pulse turn-ons and -offs is not taken into account.

In our results, the field strength regime where $`\widehat{E}^{(1)}/\omega _1,\widehat{E}^{(2)}/\omega _2<1`$ in Fig. 10 may be identified as the “death valley” of high ionization probability. Of course, the survival probability increases towards unity as $`\widehat{E}^{(1)},\widehat{E}^{(2)}0`$. However, this has been suppressed in the plot because it blocked the visibility of the results in the stabilization regime in which we are interested. Given a certain $`\widehat{E}^{(1)}`$, adding $`\widehat{E}^{(2)}`$ is, from a stabilization point of view, almost always beneficial, provided that $`\widehat{E}^{(2)}`$ is not too strong so that shake-off dominates. If, on the other hand, we choose a certain $`\widehat{E}^{(2)}`$ and add a $`\omega _1`$-field the survival probability may be enhanced or diminished, depending on whether $`\widehat{E}^{(2)}`$ lies in its single-color stabilization regime or near its single-color “death valley.” The case of a rosette-like stabilized probability density where the excursion amplitudes $`\widehat{\alpha }`$ are equal leads only to a relatively small reduction in the survival probability with respect to the optimally stabilized $`\omega _2`$-field alone (around $`\widehat{E}^{(2)}=12`$ in Fig. 10). In contrast, the equal intensity result remains close to the low-frequency single-color result where $`\widehat{E}^{(2)}=0`$. Note that maximum stabilization is achieved for modest values of $`\widehat{\alpha }=\widehat{E}/\omega ^2`$ where the probability density splitting, i.e., trichotomy, quatrochotomy, etc., is not yet developed.

### IV.4 ATI in circular two-color stabilization

Finally, let us briefly discuss the two-color photoelectron spectra. The single color case has been analyzed in subsection III.4.

In Fig. 11 the two-color spectra are compared with the corresponding single-color results where either the $`\omega _1=1`$-field or the $`\omega _2=2`$-field was present. The laser intensities were chosen not deep in the stabilization regime because at such high laser intensities the clear ATI-structure gets corrupted by peak broadening and modulations (cf. Fig. 4).

The ATI peaks in the two-color spectra are located at higher energies than those of the corresponding single-color results. This is due to the fact that the ground state in the time-averaged KH potential experiences a stronger up-shift as compared to the single-color KH potentials. It would be interesting to analyze the ATI peak positions in the two-color case in the same way as it was done in Sec. III.4 for single-color fields. However, to our knowledge two-color Floquet energies for high-frequency, circular polarization are not yet published.

As far as the peak strengths are concerned the two-color spectrum is, of course, not simply the sum of the two single-color spectra: the first ATI peak in the single-color $`\omega _2`$-spectrum is stronger than both, the second peak of the two-color and the single-color $`\omega _1`$-field. This is compatible with Fig. 10 from which one infers that adding to the $`\omega _2`$-field (whose intensity lies close to the “death valley”) the lower frequency $`\omega _1`$-field suppresses ionization. In the equal-$`\widehat{A}`$-result the two-color spectrum is close to the single-color, low-frequency spectrum. Instead, the two-color spectrum for equal $`\widehat{\alpha }`$ is qualitatively different from both two single-color spectra. Not only modulations appear in the two-color result which are absent in the $`\omega _1`$-spectrum but also the height of the peaks differ significantly from each other.

## V Conclusion

We have studied the dynamic stabilization of atomic hydrogen against ionization in single- and two-color, circularly polarized laser pulses by solving numerically the three-dimensional, time-dependent Schrödinger equation. In the single-color case we confirmed the “hula-hoop”-like dynamics for sufficiently high excursions $`\widehat{\alpha }`$. Ionization rates were compared with those from (high-frequency) Floquet theory, and good agreement was found. The positions of above-threshold ionization (ATI) peaks in our numerical photoelectron spectra were analyzed. It was found that the ATI peaks move with increasing laser intensity towards higher energies, in contrast to what happens in intense, low-frequency fields. This behavior can be explained by the laser field-induced up-shift of the initially populated, field-free ground state.

For two-color laser fields of opposite circular polarization we presented stabilized probability density structures. For laser frequencies $`\omega _1`$ and $`\omega _2=n\omega _1`$, $`n=2,3,\mathrm{}`$ and sufficiently large excursion amplitudes $`(n+1)`$ distinct probability density peaks were observed, in accordance to what one expects by virtue of the two-color, time-averaged Kramers-Henneberger potentials. As the generalization of the well-known “dichotomy” in linearly polarized laser fields, we called this multiple probability density peak-splitting “trichotomy,” “quatrochotomy,” “pentachotomy” etc.

The survival probability after two-color pulses of frequencies $`\omega _1=1`$ and $`\omega _2=2`$ was calculated as a function of the two peak field strengths. As expected, highest stabilization is achieved with the high-frequency field alone for a certain, optimal laser intensity. However, adding the second, low-frequency field, even if its intensity falls into the “death valley,” fortunately does not drastically reduce stabilization, as one might expect from a naive superposition point-of-view.

###### Acknowledgements.

The authors thank Dr. R. Potvliege for providing electronically the Floquet rates and energy shifts to which comparison was made in Sec. III. The permission to run our codes on the Linux cluster at PC<sup>2</sup> in Paderborn, Germany, is gratefully acknowledged. This work was supported in part by the INFM through the Advanced Research Project CLUSTERS.


format: text

text:
# Molecular gas in spiral galaxies: a new warm phase at large galactocentric distances?

## 1 Introduction

The use of 12CO, the most abundant molecule after H2 itself, to trace H2 gas is a widely used practice since the discovery of bright 12CO J=1–0 emission in Orion (Wilson, Jefferts, & Penzias 1970). The most easily excited H2 transition, the quadrupole S(0): $`\mathrm{J}_\mathrm{u}\mathrm{J}_\mathrm{l}=20`$, corresponds to $`\mathrm{\Delta }\mathrm{E}_{20}/\mathrm{k}510`$ K, much too high to be significantly populated in the cool ISM ($`\mathrm{T}_{\mathrm{kin}}1050`$ K) revealed by the 12CO J=1–0 transition ($`\mathrm{\Delta }\mathrm{E}_{10}/\mathrm{k}5.5`$ K). Furthermore the large optical depths of the latter ($`\tau _{10}510`$) reduce its critical density of $`\mathrm{n}_{\mathrm{cr}}^{(10)}380`$ cm-3 (for $`\mathrm{T}_{\mathrm{kin}}=50`$ K; Kamp & Zadelhoff 2001) to $`\mathrm{n}_{\mathrm{cr}}^{(10)}(1\mathrm{e}^{\tau _{10}})\tau _{10}^14080\mathrm{cm}^3`$, similar to the lowest average densities observed in Giant Molecular Clouds (GMCs). Hence the 12CO J=1–0 transition is expected to be excited even in the coldest and most diffuse regions of the molecular ISM.

However its large optical depths (e.g. Martin, Sanders, & Hills 1984; Sage & Isbell 1991; Falgarone 1998), do not allow a straightforward use of its luminosity together with a $`[^{12}\mathrm{CO}/\mathrm{H}_2]`$ abundance ratio for H2 mass estimates. The widely used $`\mathrm{X}_{\mathrm{CO}}=\mathrm{M}(\mathrm{H}_2)/\mathrm{L}_{\mathrm{CO}}`$ factor (where $`\mathrm{L}_{\mathrm{CO}}`$ is the velocity/area-integrated brightness temperature of 12CO J=1–0) is based on the notion of an effectively optically thin medium, where the 12CO line emission arises from an ensemble of small, radiatively de-coupled cells not overlapping in space or velocity. Then the high optical depths arise locally within the individual cells and the observed wide line profiles result from their macroscopic motions rather than their intrinsic, much narrower linewidths (e.g. Martin et al. 1984; Young & Scoville 1991 and references therein). This picture, along with the observed virialization of molecular clouds (see e.g. Larson 1981) over a wide range of scales ($`0.1100`$ pc) allows a statistical derivation of $`\mathrm{X}_{\mathrm{CO}}`$ as an ensemble-average that is relatively insensitive to local molecular cloud conditions (e.g. Dickman, Snell, & Schloerb 1986; Young & Scoville 1991; Bryant & Scoville 1996) and, for scales $`10`$ pc, is expected to yield the correct H2 mass to within a factor of two.

The existence of such cells is supported by high angular and velocity resolution studies of galactic clouds (e.g. Falgarone et al. 1998; Tauber 1996) while the statistical robustness of $`\mathrm{X}_{\mathrm{CO}}`$ for a variety of conditions has been verified using both observational (e.g. Young & Scoville 1991 and references therein), and theoretical methods (e.g. Kutner & Leung 1985; Dickman, Snell & Schloerb 1986; Maloney & Black 1988; Wolfire, Hollenbach, & Tielens 1993; Sakamoto 1996; Bryant & Scoville 1996). However these studies also suggest a likely failure of this method in two particular environments, namely a) galactic and/or starburst nuclei where non-virial motions cause the luminous 12CO emission to overestimate the H2 mass (e.g. Dahmen et al. 1998; Downes & Solomon 1998), b) metal-poor regions where it seriously underestimates the H2 mass (e.g. Maloney & Black 1988; Arimoto, Sofue, & Tsujimoto 1996; Israel 1988, 1997, 1999). The last case is the most serious one since then 12CO emission is usually faint and thus precludes the observations of more optically thin isotopes like 13CO, C18O that could yield the H2 mass without using the $`\mathrm{X}_{\mathrm{CO}}`$ factor.

In this paper we show that the decreasing metallicity with galactocentric radius in disks and a larger H2 formation rate open up the possibility of large amounts of H2 gas being there hitherto undetected by the conventional method. This phase will be warm, diffuse and devoid of CO, and may contribute significantly to the pressure and total mass of the disk.

## 2 Spiral disks: H2 formation in the metal-poor parts

Molecular hydrogen in interstellar conditions forms from HI association onto dust grains (e.g. Gould & Salpeter 1963) and deep sub-mm imaging shows the grain distribution extending beyond the optical and well into the HI disk (Nelson, Zaritsky, & Cutri 1998; Xylouris et al. 1999), and possibly out to $`2\mathrm{R}_{25}`$ (Thomas et al. 2001, 2002; Guillandre et al. 2001). Thus the basic “ingredients” for H2 formation can be found to the outermost parts of the disk as defined by the HI distribution itself. Using the results of Federman, Glassgold, & Kwan (1979), with an exponent of m=1.5 for the H2 self-shielding function, we find that an HI cloud of radius $`\mathrm{R}_{\mathrm{cl}}`$ (pc), average volume density n and temperature T, illuminated by a FUV field strength $`\mathrm{G}_{}`$ ($`\mathrm{G}_{}`$(solar)=1) starts turning molecular when $`s>1`$, where

$$s=\frac{\mathrm{n}}{\mathrm{n}_{\mathrm{min}}}=\mathrm{R}_{\mathrm{cl}}^{2/5}\left(\frac{\mathrm{n}}{65\mathrm{c}\mathrm{m}^3}\right)\left(\frac{\mathrm{T}}{100\mathrm{K}}\right)^{3/10}\left(\frac{\mu \mathrm{z}}{\mathrm{G}_{}}\right)^{3/5}.$$
(1)
The density $`\mathrm{n}_{\mathrm{min}}`$ is the minimum required for the onset of an HI-to-H2 transition and $`\mathrm{z}=\mathrm{Z}/\mathrm{Z}_{}`$ is the normalized metallicity. We kept the explicit temperature dependence of the H2 formation rate $`\mathrm{R}_\mathrm{f}=\mathrm{S}_\mathrm{f}\mathrm{T}^{1/2}`$, assumed that it scales linearly with total grain surface and thus metallicity, and included the factor $`\mu =\mathrm{S}_\mathrm{f}/\mathrm{S}_\mathrm{f}^{(\mathrm{o})}`$, to parametrize for a value different than the currently adopted $`\mathrm{S}_\mathrm{f}^{(\mathrm{o})}=3\times 10^{18}`$ cm3 s-1 K-1/2 (Jura 1975a).

Over an HI disk it is $`\mathrm{T}(\mathrm{HI})80150`$ K, with a slight increase with galactocentric distance $`\mathrm{R}_{\mathrm{gal}}`$ (e.g. Braun 1997), a temperature range typical of the warm neutral halos around molecular clouds (e.g. Andersson & Wannier 1993) where the HI-H2 transition zone lies. Thus no significant variation of $`s`$ is expected because of a strongly varying $`\mathrm{T}(\mathrm{R}_{\mathrm{gal}})`$. The influence of the $`\mathrm{z}/\mathrm{G}_{}`$ ratio is less straightforward since significant metallicity gradients exist in spiral disks (e.g. Henry 1998; Garnett 1998 and references therein) along with the well known exponential optical (and hence FUV) light profiles. Their effects have been thoroughly modeled by Elmegreen (1989, 1993) and applied to spiral galaxies by Honma, Sofue, & Arimoto (1995). The latter find the metallicity gradient to be the main cause of the so called molecular front, the steep H2-to-HI transition at a particular $`\mathrm{R}_{\mathrm{tr}}`$ deduced for most spirals from CO imaging. Evidently a larger H2 formation rate will alter this picture by “pushing” this front to encompass a larger portion of a typical disk. Rates that are $`48`$ times higher than the standard value are compatible with observationally deduced upper limits (Jura 1974), since they always contain an $`\mathrm{n}\mathrm{S}_\mathrm{f}`$ “degeneracy” that can bias $`\mathrm{S}_\mathrm{f}`$ estimates towards lower values (e.g. Jura 1975a,b; van Dishoeck & Black 1986). Recent results from ISO observations of H2 rotational lines from photodissociation regions (Habart et al. 2000; Li et al. 2002) strongly favor $`\mu =\mathrm{S}_\mathrm{f}/\mathrm{S}_\mathrm{f}^{(\mathrm{o})}5`$ (see also Sternberg 1988 for early indications of such high values), which will raise the value of $`s`$ at all $`\mathrm{R}_{\mathrm{gal}}`$.

This rescaling has no effect for those regions of the disk where the old values are already $`s>1`$ and which CO imaging shows them to be indeed H2-rich, but now a larger portion of the disk is expected to have $`s>1`$ and thus conditions favorable for an HI-to-H2 transition. A rough estimate of this new radius can be obtained from the results of Honma et al. (1995) that show this transition to occur very sharply at a particular $`\mathrm{R}_{\mathrm{tr}}`$. The expected FUV volume emissivity $`\mathrm{j}_{}`$, its dust absorption coefficient $`\mathrm{k}_\mathrm{d}`$, and the metallicity follow a similar spatial distribution (Honma et al. 1995) hence, neglecting the H2 self-shielding contributions to the total absorption coefficient corresponds to $`\mathrm{z}/\mathrm{G}_{}=\mathrm{z}/(\mathrm{j}_{}/\mathrm{k}_\mathrm{d})\mathrm{z}`$. For an exponential disk profile and $`\mathrm{s}(\mathrm{R}_{\mathrm{tr}})1`$ it is

$$\mathrm{R}_{\mathrm{tr}}(\mu )=\mathrm{R}_{\mathrm{tr}}(1)+(\mathrm{ln}\mu )\mathrm{R}_\mathrm{e},$$
(2)
where $`\mathrm{R}_\mathrm{e}`$ is the scale-length of the exponential distribution of optical light and metallicity. Thus the molecular front will now reside at least $`(\mathrm{ln5})\mathrm{R}_\mathrm{e}1.6\mathrm{R}_\mathrm{e}`$ further, past that inferred by CO imaging. In typical disks the later is at $`\mathrm{R}_{\mathrm{tr}}(1)=(0.51)\mathrm{R}_\mathrm{e}`$ (e.g. Young & Scoville 1991; Regan et al. 2001), thus $`\mathrm{R}_{\mathrm{tr}}(\mu )(22.5)\mathrm{R}_\mathrm{e}`$ places the molecular front well inside a typical HI disk. A flattening of the abundance gradient, and thus of dust grain surface per H nuclei, at large $`\mathrm{R}_{\mathrm{gal}}`$ (Henry 1998) and H2 self-shielding will only strengthen the above arguments by yielding still larger values of $`s`$ in the outer parts of the disk.

Therefore it now seems that the molecular front in spirals, deduced from CO imaging, is in reality a (low) metallicity-induced 12CO detection threshold, and the necessary conditions for H2 formation exist over a much larger part of a typical spiral disk. This is not surprising since it is known that $`\mathrm{X}_{\mathrm{CO}}\mathrm{z}^1\mathrm{z}^2`$ (e.g. Arnault et al. 1988; Israel 1997), which alone implies a $`\mathrm{R}_{\mathrm{tr}}`$ larger than that found by using 12CO imaging and a constant value of $`\mathrm{X}_{\mathrm{CO}}`$.

### 2.1 The role of pressure in the HI/H2 transition

In obtaining the previous result we assumed no significant changes of $`(\mathrm{n},\mathrm{R}_{\mathrm{cl}})`$ of the various clouds with $`\mathrm{R}_{\mathrm{gal}}`$, which leaves only the $`\mathrm{z}/\mathrm{G}_{}`$ variable to define $`\mathrm{R}_{\mathrm{tr}}`$. In a forthcoming paper we will present a more detailed analysis of the HI-H2 transition without such simplifications, yet a re-formulation of $`s`$ in terms of cloud mass M and boundary pressure $`\mathrm{P}_\mathrm{e}`$ suggests that the aforementioned assumptions do not subtract from the generality of our result. Using the results of Elmegreen (1989) (his S function and $`s`$ are related simply as $`s=\mathrm{S}^{3/5}`$), we obtain

$$s1.1\left(\frac{\mathrm{P}_\mathrm{e}}{10^4\mathrm{cm}^3\mathrm{K}}\right)^{13/20}\left(\frac{\mathrm{M}}{10^5\mathrm{M}_{}}\right)^{3/10}\left(\frac{\mathrm{T}}{100\mathrm{K}}\right)^{3/10}\left(\frac{\mu \mathrm{z}}{\mathrm{G}_{}}\right)^{3/5}.$$
(3)
It is clear that unless M is a strong function of $`\mathrm{R}_{\mathrm{gal}}`$, a radially declining $`\mathrm{P}_\mathrm{e}`$ is the only other important factor besides $`\mu \mathrm{z}/\mathrm{G}_{}`$ in determining $`\mathrm{R}_{\mathrm{tr}}`$. The increase of s with decreasing M may seem counterintuitive but it is due to the fact that small clouds have higher average densities than big ones with the same boundary pressure. This stems from the widespread virialization of most of the clouds, which in combination with the observed velocity dispersion vs cloud size relation (e.g. Larson 1981), yields $`\mathrm{n}\mathrm{R}_{\mathrm{cl}}^\alpha `$ ($`\alpha 1`$). In the metal-rich part of a typical spiral, where CO can be relied upon as an H2 gas mass tracer, the bulk of the H2 mass is distributed in a cloud hierarchy with $`\mathrm{M}10^5\mathrm{M}_{}`$ being typically the upper limit. The value of the latter is not observed to change with galactocentric distance (e.g. Elmegreen & Elmegreen 1987) except in galactic centers where tidal effects may truncate such a cloud hierarchy at higher minimum densities and cloud masses (e.g. Maloney & Black 1988).

The measurement of $`\mathrm{P}_\mathrm{e}`$, especially as a function of $`\mathrm{R}_{\mathrm{gal}}`$, is not easy and deducing it by assuming equipartition with the random cloud motions (e.g. Honma et al. 1995) relies on a known total gas surface density, which is what needs to be determined. This task is easier for the Galaxy whose disk we consider typical in that respect. In this case a pressure probe sampling its outer parts, namely diffuse CO-luminous molecular clouds, suggests pressures $`\mathrm{P}_\mathrm{e}(12\mathrm{k}\mathrm{p}\mathrm{c})(12)\times 10^4\mathrm{cm}^3\mathrm{K}`$ (Heyer, Carpenter, & Snell 2001), which are higher than expected and similar to the Solar neighborhood value. Further out pressures of $`\mathrm{P}_\mathrm{e}(17\mathrm{k}\mathrm{p}\mathrm{c})5\times 10^3\mathrm{cm}^3\mathrm{K}`$ are typical, with few clouds reaching up to $`10^5\mathrm{cm}^3\mathrm{K}`$ (Brand, & Wouterloot 1995). It is straightforward to show that since $`0.25\mathrm{z}/\mathrm{G}_{}1`$ (at $`\mathrm{R}_{\mathrm{gal}}17`$ kpc), for $`\mu =5`$, such $`\mathrm{P}_\mathrm{e}`$ values yield $`s1`$. Thus the H2 formation conditions are fulfilled out to $`2.5\mathrm{R}_\mathrm{e}(17\mathrm{kpc})`$. This radius is well past the location of the 12CO J=1–0 brightness peak at $`5`$ kpc or any significant 12CO emission (e.g. Blitz 1997) and well into the HI disk. Interestingly HII regions, another potential pressure probe, may also hint at larger than expected pressures in the outer Galaxy (Rudolph et al. 1996) and dwarf irregulars (Elmegreen, & Hunter 2000). In both of these cases low metallicities make 12CO a poor indicator of H2 gas which may thus be responsible for the high pressures.

## 3 Molecular gas at large galactocentric radii: CO-deficient and warm

The arguments for favorable conditions for an HI-to-H2 phase transition presented here do not and could not yield a mass estimate for the associated H2 gas. Nevertheless the simple fact that this gas resides untraceable by CO in the metal-poor parts of the disk can offer significant insights regarding its state. An obvious consequence of the low metallicities is reduced dust shielding which allows even the weaker FUV radiation field at large $`\mathrm{R}_{\mathrm{gal}}`$’s to dissociate the mainly dust-shielded CO and confine it to much smaller volumes while leaving the largely self-shielding H2 intact. A typical range of $`\mathrm{G}_{}0.30.6`$ and an $`\mu 5`$ yields a minimum H2 formation density of $`\mathrm{n}_\mathrm{f}57\mathrm{cm}^3`$ for cloud sizes of $`1050`$ pc. Observations of CO-luminous molecular clouds yield average densities of $`50500\mathrm{cm}^3`$ irrespective of galactocentric distance, hence $`\mathrm{n}5500\mathrm{cm}^3`$ and $`\mathrm{z}0.15`$ are appropriate inputs for Photodissociation Region models (PDR) of the molecular gas at large $`\mathrm{R}_{\mathrm{gal}}`$. Similar modeling of metal-poor gas has been performed before but at much higher UV fields and gas densities than here (e.g. Pak et al. 1998).

We investigate the gas abundance of warm H2 using a PDR code that models the chemistry of a diffuse cloud represented as a semi-infinite slab, with a pseudo-time-dependent multi-point chemical simulation (Viti et al. 2001). The UMIST chemical database (Millar et al. 1997) was used with some changes (see Viti et al. 2001). The heating, the cooling and CO emission is as described by Taylor et al (1993); we have however made the following changes: i) the thermal balance is achieved using Ridder’s method, a fast non-linear convergence algorithm, ii) the heating mechanisms have been revised and now include also large grains and PAHs (Weingartner & Draine 2001), as well as UV pumping (not important for low G’s and densities) and chemical heating, iii) the collisional rates have been updated to include all the dominant collisional partners for such an ISM phase (H, H2, He). All our runs were performed assuming a microturbulent line formation with a typical width of $``$3 km s-1, and following Tielens & Hollenbach 1985 the level populations are updated from the outside in. The heating caused by the decay of the turbulent gas is negligible.

The results are summarized in Figure 1 where it becomes apparent that in diffuse gas of sub-solar metallicity the CO-deficient H2 gas dominates the cloud. This radiation-driven spatial segregation also maintains a temperature difference in which the pure H2 phase stays always warmer than the CO-luminous one. This simply reflects the different dominant coolants present, namely 12CO and its isotopes for the CO-luminous phase, and CII, H2 for the CO-deficient one. Average temperatures of $`\mathrm{T}_\mathrm{k}1520`$ K are observed for most quiescent CO-luminous molecular clouds in the disk (with $`\mathrm{T}_\mathrm{k}3060`$ K found only in star-forming warm “spots”) and remain almost constant with $`\mathrm{R}_{\mathrm{gal}}`$. In the Galaxy this is beautifully demonstrated by COBE observations of CO J+1$``$J (J=0–7) lines that show invariant relative strengths suggesting $`\mathrm{T}_\mathrm{k}20`$ K over the entire disk except the Galactic center (Fixsen, Bennett, & Mather 1999). This temperature is similar to those of several individual molecular clouds, some as far as $`\mathrm{R}_{\mathrm{gal}}30`$ kpc (Digel, de Geus, & Thaddeus 1994), and also similar to that measured by COBE for the bulk of the dust in the Galaxy (Sodroski et al. 1994; Lagache et al. 1998) and by sub-mm imaging in other spirals (e.g. Alton et al. 1998). Such gas temperatures are indeed expected for the dust-shielded ($`\mathrm{A}_\mathrm{v}4`$) CO-rich H2 clouds, where CO cooling (Goldsmith, & Langer 1978; Neufeld, Lepp, & Melnick 1995) lowers $`\mathrm{T}_\mathrm{k}`$ to that of the concomitant dust whose temperature is weakly dependent to the ambient ISRF and thus to $`\mathrm{R}_{\mathrm{gal}}`$.

For the unshielded regions of the molecular ISM, now expected to contain the bulk of the H2 gas mass in the metal-poor parts of spiral disks, lack of the CO coolant will keep $`\mathrm{T}_\mathrm{k}>\mathrm{T}_{\mathrm{dust}}`$ (e.g. Tielens, & Hollenbach 1985). Its temperature range can be found from

$$\mathrm{\Lambda }_{\mathrm{CII}}+\mathrm{\Lambda }_{\mathrm{H}_2}=\mathrm{\Gamma }_{\mathrm{ph}}+\mathrm{\Gamma }_{\mathrm{CR}},$$
(4)
where the cooling $`(\mathrm{\Lambda })`$ due to the CII and H2 lines is balanced by the photoelectric heating from dust grains $`\mathrm{\Gamma }_{\mathrm{ph}}`$ (e.g. Tielens & Hollenbach 1985) and cosmic rays. Unlike solar metallicity environments, the cooling due to H2 lines may now be important and is included.

The CII emission emanates from a 2-level system, hence following Hollenbach, Takahashi, & Tielens (1991) and for optically thin emission it is

$$\mathrm{\Lambda }_{\mathrm{CII}}=1.82\times 10^{23}\mathrm{zn}\left[1+\frac{1}{2}\mathrm{e}^{92/\mathrm{T}}\left(1+\frac{\mathrm{n}_{\mathrm{cr}}}{\mathrm{n}}\right)\right]^1\mathrm{ergs}\mathrm{s}^1,$$
(5)
where $`\mathrm{n}_{\mathrm{cr}}3\times 10^3\mathrm{cm}^3`$ is the transition’s critical density, and all carbon is assumed in the form of CII, mixed with a fully molecular gas phase at an abundance of $`[\mathrm{C}/\mathrm{H}]=3\times 10^4\mathrm{z}`$.

For the expected range of densities and $`\mathrm{T}_\mathrm{k}200`$ K, (a typical upper limit for the Cold Neutral Medium (CNM) in spiral disks), only the S(0) and S(1) lines rotational lines of H2 can be appreciably excited and thus contribute to the cooling. Moreover, since no radiative or collisional processes “mix” ortho-H2 (J=odd) and para-H2 (J=even) states we can consider S(0) and S(1) practically arising from two different molecules, each approximated by a 2-level system, therefore

$$\mathrm{\Lambda }_{\mathrm{H}_2}=2.06\times 10^{24}\frac{\mathrm{n}}{1+\mathrm{r}_{\mathrm{op}}}\left[1+\frac{1}{5}\mathrm{e}^{510/\mathrm{T}}\left(1+\frac{\mathrm{n}_{\mathrm{cr}}(0)}{\mathrm{n}}\right)\right]^1\left(1+\mathrm{R}_{10}\right)\mathrm{ergs}\mathrm{s}^1,$$
(6)
where $`\mathrm{r}_{\mathrm{op}}=(\mathrm{n}_3+\mathrm{n}_1)/(\mathrm{n}_2+\mathrm{n}_0)`$ is the ortho-to-para H2 ratio, $`\mathrm{n}_{\mathrm{cr}}(0)54\mathrm{cm}^3`$ is the critical density of the S(0) transition and $`\mathrm{R}_{10}=\mathrm{\Lambda }_{\mathrm{S}(1)}/\mathrm{\Lambda }_{\mathrm{S}(0)}`$, namely

$$\mathrm{R}_{10}=26.8\mathrm{r}_{\mathrm{op}}\left[\frac{1+\frac{1}{5}\mathrm{e}^{510/\mathrm{T}}\left(1+\frac{\mathrm{n}_{\mathrm{cr}}(0)}{\mathrm{n}}\right)}{1+\frac{3}{7}\mathrm{e}^{845/\mathrm{T}}\left(1+\frac{\mathrm{n}_{\mathrm{cr}}(1)}{\mathrm{n}}\right)}\right].$$
(7)
The critical density of the S(1) transition is $`\mathrm{n}_{\mathrm{cr}}(1)10^3\mathrm{cm}^3`$, and this term dominates for the warmest and densest gas in the chosen parameter space.

We adopt the latest photoelectric heating model of Weingartner & Draine (2001) with $`\mathrm{\Gamma }_{\mathrm{ph}}`$ proportionally scaled by the normalized metallicity factor z and parameters corresponding to carbonaceous grains (the last two models in their Table 2). The assumed electron abundance is the average for the CNM, under similar set of $`(\mathrm{z},\mathrm{G}_{})`$, i.e. $`\mathrm{n}_\mathrm{e}/2\mathrm{n}5\times 10^4`$ (Wolfire et al. 1995). Finally we include a heating factor from cosmic rays $`\mathrm{\Gamma }_{\mathrm{CR}}=6.4\times 10^{28}\mathrm{n}`$ ergs s-1 (Goldsmith & Langer 1978), which is $`10\%`$ of the photelectric term.

For $`\mathrm{G}_{}=0.30.6`$, $`\mathrm{z}=0.15`$ and $`\mathrm{r}_{\mathrm{op}}=13`$, both a diffuse ($`520`$ cm-3) and a dense ($`100500`$ cm-3) H2 gas phase satisfy (4) with the former being warmer ($`50100`$ K) than the latter ($`2030`$ K). At least for the Galaxy, these solutions can be further constrained by the values of the total pressure $`\mathrm{P}_\mathrm{e}(520)\times 10^3`$ cm-3 K, at $`\mathrm{R}_{\mathrm{gal}}8`$ kpc, where

$$\mathrm{P}_\mathrm{e}=\mathrm{nk}_\mathrm{B}\mathrm{T}+56.8\mathrm{n}\sigma _\mathrm{z}^2,$$
(8)
(e.g. van Dishoeck & Black 1986). The second term expresses the non-thermal part of the pressure due to the macroscopic motions of the gas, which are assumed of turbulent origin with $`\sigma _z`$ being the one-dimensional FWHM of a gaussian distribution. The velocity dispersions of CO-bright molecular clouds or HI in face-on spirals yield $`\sigma _z5`$ km s-1 (Combes 1999 and references therein) hence the aforementioned constraints on the total disk pressure correspond to $`\mathrm{n}515`$ cm-3 (for $`\mathrm{T}100`$ K). Thus a metal-poor, warm and diffuse H2 gas in the outer parts of spiral disks is likely to exist, and may indeed be the type of gas detected through its S(0) line with ISO in the edge-on spiral NGC 891 (Valentijn & van der Werf 1999), with a reported temperature of $`\mathrm{T}90`$ K and a mass that outweighs that of HI by factors of $`515`$.

Several issues regarding such a gas phase remain open, namely its thermal stability, low star formation rates, its exact relation to the HI gas distribution and the dynamical evolution of such potentially massive gaseous disks. Moreover an enhanced H2 formation rate, if indeed at work, may alter relevant aspects of the “standard” ISM picture (e.g. gas heating) that must be taken into account in a detailed analysis. These and other issues will be the subject of a future paper.

## 4 Conclusions

In the light of mounting evidence that the standard method of using rotational line emission of the 12CO molecule to trace H2 underestimates its mass at low metallicities we examined the state of such gas in the disks of spiral galaxies at large galactocentric distances. Our results can be summarized as follows,

A high H2 formation rate along with the presence of dust throughout a typical HI disk, both suggested by recent observations, raise the possibility of an extended H2 gas phase, well past the one inferred by the 12CO brightness distribution.

This phase is likely to be warm and diffuse and it may be responsible for the observed high total pressures at large galactocentric distances in the Galaxy. In the case of NGC 891 such a massive gaseous reservoir may have already been detected through its S(0) line emission.

We thank the anonymous referee for useful comments that improved the present work. P. P. P. would also like to thank Paula Gutierrez Ruiz for many inspiring conversations that helped this work through.


text:
# Shear localization in a model glass

## Abstract

Using molecular dynamics simulations, we show that a simple model of a glassy material exhibits the shear localization phenomenon observed in many complex fluids. At low shear rates, the system separates into a fluidized shear-band and an unsheared part. The two bands are characterized by a very different dynamics probed by a local intermediate scattering function. Furthermore, a stick-slip motion is observed at very small shear rates. Our results, which open the possibility of exploring complex rheological behavior using simulations, are compared to recent experiments on various soft glasses.

Shear localization is a commonly observed phenomenon in the rheology of complex fluids. Over some range of shear rates, a fluid undergoing simple shear flow in, say, a Couette cell tends to separate into bands parallel to the flow direction, with high shear rate regions coexisting with smaller shear rate regions. In some cases olmsted , this shear-banding phenomenon can be understood in terms of underlying structural changes in the fluid, analogous to a first order phase transition. In other systems, however, no such changes are evident, and coexistence appears between a completely steady region (zero shear rate) and a sheared, fluid region. This second type of behavior has been observed Chen ; Pignon::JRheo40::1996 ; Losert::PRL85::2000 ; Debregeas::PRL87::2001 ; Coussot-Raynaud-et-al::PRL88::2002 , in particular in systems of the so-called ‘soft glass’ type Sollich . Such systems include dense colloidal pastes, granular materials, emulsions, etc. Their rheological behavior is essentially determined by the competition between an intrinsic slow dynamics and the acceleration caused by the external flow Sollich ; BB::PRE61::2000 ; BB2 . The large time scales inherent to the glassy state manifest themselves in the non-linear character of the rheological properties as a function of the shear rate $`\gamma `$: existence of a yield stress (the system does no flow until the stress $`\sigma `$ exceeds a threshold value $`\sigma _\text{y}`$) and non-linear flow curves $`\sigma =\sigma (\gamma )`$, leading to a shear rate dependent viscosity.

The observation of strong heterogeneities in the flow of such systems suggests that a global flow curve is not sufficient to fully characterize the flow behavior. Indeed, a simple shear-thinning behavior would in general imply homogeneous shear flow in a planar Couette cell, since the shear stress is constant across the cell. More generally, it remains to clarify if these observations in systems with very different microscopic interactions are intrinsic to glassy dynamics, that is, if a generic scenario for inhomogeneous shear flow can be proposed, as attempted in several recent studies Dhont ; Picard .

In order to investigate this issue, we performed molecular dynamics simulations of a simple and generic glass forming system, consisting of a binary Lennard-Jones mixture. Since our work is done on a simple liquid, it can be seen as a link between experiments on complex materials and pure phenomenology.

The present model system has been extensively studied BB::PRE61::2000 ; BB2 ; Kob-Andersen . It is known to exhibit, in the bulk state, a computer glass transition (in the sense that the relaxation time becomes larger than typical simulation times) at a temperature $`T_\text{c}0.435`$ for a 80:20 mixture at a density $`\rho =1.2`$ (in Lennard-Jones units Kob-Andersen ). All our simulation parameters are chosen equal to those in Refs. BB2 ; Kob-Andersen . Since our aim is to analyse possible flow heterogeneities, we do not impose a constant velocity gradient over the system as was done in Ref. BB2 , where a homogeneous shear flow was imposed through the use of Lee-Edwards boundary conditions. Rather we confine the system between two solid walls, which will be driven at constant velocity. By doing so, we mimic an experimental shear cell, without imposing an homogeneous flow.

We first equilibrate a large simulation box with periodic boundary conditions in all directions, at $`T=0.5`$. The system is then quenched to a temperature below $`T_\text{c}`$. Most of the results discussed in this work correspond to $`T=0.2`$. At this temperature, the structural relaxation times are orders of magnitude larger than at $`T_\text{c}`$. On the time scale of computer simulations, the system is in a glassy state, in which its properties slowly evolve with time (aging). After a time of $`t=4.10^4`$ \[$`2.10^6`$ MD steps\], we create 2 parallel solid boundaries by freezing all the particles outside two parallel $`xy`$-planes at positions $`z_{\text{wall}}=\pm L_z/2`$. For each computer experiment, 10 independent samples (each containing 4800 particles) are prepared using this procedure. The cell containing the fluid particles has dimensions $`L_x=L_y=10`$ and $`L_z=40`$.

An overall shear rate is imposed by moving in the $`x`$-direction all the atoms of, say, the left wall ($`z_{\text{wall}}=20`$) with a strictly constant velocity of $`U_{\text{wall}}`$. This defines the total shear rate $`\gamma _{\text{tot}}=U_{\text{wall}}/L_z`$. Velocity profiles are recorded discarding transients (of typical duration $`1/\gamma `$) associated with the start of the shear motion. In order to remove the viscous heat created by the shear, the system is coupled to a heat bath. For this purpose, we divide the system into parallel layers of thickness $`dz=0.25`$ and rescale (once every 10 integration steps) the $`y`$-component of the particle velocities within the layer, so as to impose the desired temperature $`T`$. Such a local treatment is necessary to keep a homogeneous temperature profile when flow profiles are heterogeneous. To check for a possible influence of the thermostat, we compared, for an extremely low shear rate ($`\gamma _{\text{tot}}=U_{\text{wall}}/L_z=10^5`$), these results with the output of a simulation where the inner part of the system was unperturbed and the walls were thermostatted instead. Both methods give identical results, indicating that the results are not affected by the thermostat.

Our main observation is reported in Fig. 1. We find that, at low overall shear rates $`\gamma _{\text{tot}}`$, the system separates into a homogeneously sheared band and a part which is essentially unsheared. Note that, the velocity profiles are not symmetric with respect to the midplane. Rather, the shear band is localized close to one of the walls, so that a single interface between sheared and steady region is formed. The symmetry between the two walls, which is a result of Galilean invariance, is restored only on average, the shear band occurring equally likely on both sides of the simulation cell (see Fig. 1). A similar behavior has been observed in various experiments Chen ; Pignon::JRheo40::1996 ; Losert::PRL85::2000 ; Debregeas::PRL87::2001 ; Coussot-Raynaud-et-al::PRL88::2002 .

In some cases, we observed oscillations of the velocity profile between the two mentioned solutions. This effect is possibly related to the finite size of the simulation box leading to a finite probability for the band to oscillate. In contrast, this probability would be zero for (very large) experimental systems, thus stabilizing one of the solutions. A discussion of this interesting aspect is beyond the scope of this paper, and we postpone it to future work.

As shown in Fig. 1, the thickness $`h`$ of the sheared region depends on the wall velocity $`U_{\text{wall}}`$. For very small $`U_{\text{wall}}`$, $`h`$ is of the order of a few atomic diameters and varies only slightly with the wall velocity. As a consequence, the shear rate inside the sheared region, of order $`U_{\text{wall}}/h`$, decreases as $`U_{\text{wall}}`$ is reduced. As $`U_{\text{wall}}`$ is further increased, $`h`$ increases to reach the full slab size, $`h=L_z`$, at a given wall velocity, $`U_c`$. For $`U_{\text{wall}}>U_c`$ the velocity profile is therefore linear, as is oberved in simple fluids. These findings are in qualitative agreement with experimental results on clay suspensions Pignon::JRheo40::1996 . The rather slow variation of $`h`$ with respect to a change of $`\gamma _{\text{t}ot}`$ at small wall velocities should, however, be contrasted to reports in Coussot-Raynaud-et-al::PRL88::2002 . However, the small size of our system compared to the width of interfacial regions makes a more quantitative analysis rather difficult for the moment.

At a given global shear rate $`\gamma _{\text{tot}}=10^4(U_{\text{wall}}=0.004)`$, we also studied the influence of temperature on shear localization. We find that increasing the temperature is qualitatively similar to increasing the global shear rate: the thickness of the sheared region is an increasing function of the temperature. For example, $`h15`$ \[see Fig. 1\] at $`T=0.2`$ whereas $`h25`$ at $`T=0.4`$ and $`hL_z(=40)`$ at $`T=0.5`$.

We find no obvious structural differences between the two regions of the sample. As shown in Fig. 2 (lower panel), static properties like the density profile, shear stress, and normal pressure are found to be constant across the film shearstress .

The distinction between the two bands is therefore purely dynamical. This can be seen for instance on the layer-resolved intermediate scattering function, $`\varphi _q(t;z)=_i\mathrm{exp}[\text{i}q_y(y_i(t)y_i(0))]\delta \left(z_iz\right)`$. We choose $`q_y=7.1`$, which corresponds to the first maximum of the static structure factor. The upper panel of Fig. 2 depicts $`\varphi _q(t;z)`$ obtained from a simulation with $`\gamma _{\text{tot}}=0.83\times 10^4`$. One may first note the existence of two limiting behaviors of the correlation function, corresponding to the sheared and jammed regions, with a rather rapid change from one behavior to the other within a few layers. While in the jammed region (close to the right wall in Fig. 2) $`\varphi _q(t;z)`$ barely relaxes, the correlation function in the sheared region (close to the left wall) exhibits a two-step relaxation as observed in homogeneously sheared systems. This observation reflects the acceleration of the structural relaxation due to the flow BB2 . In the jammed region, the system behaves as a glassy solid, and $`\varphi _q(t;z)`$ does not relax to zero on the simulation time scale.

A way to quantify the relation between the structural relaxation and the variation of the shear rate across the system is to look at the quantity $`\varphi _q(\tau _0;z)`$, where $`\tau _0`$ is of order of $`1/\gamma _{\text{tot}}`$. This quantity reflects the way the system has relaxed on the time scale imposed by the global shear rate. Results for $`\varphi _q(\tau _0;z)`$ are shown in the lower panel of Fig. 2. The change of the velocity profile, when going from the sheared towards the unsheared region, is accompanied by a sharp jump in $`\varphi _q(\tau _0;z)`$ at the interface between these two regions. The profile of $`\varphi _q(\tau _0;z)`$ is very similar to that of an order parameter across an interface harrowell . Other order parameters could be proposed to characterize the local dynamics: for example, the local relaxation time of $`\varphi _q(t;z)`$, defined as the time after which $`\varphi _q(t;z)`$ goes below a fraction of unity or the local diffusion coefficient parallel to the walls. All these definitions lead to a two-phase picture, with well defined and spatially constant local characteristics in both phases.

In Fig. 3, we summarize the results by plotting the $`\sigma (\gamma )`$ flow curve of the model at $`T=0.2`$. In this figure, we indicate the flow curve obtained for a system undergoing homogeneous shear flow, taken from Ref. BB2 . We also show the points obtained for the same system driven by the boundaries, as described in the present paper. Finally, we indicate the static yield stress of the system, $`\sigma _\text{y}`$.

To obtain $`\sigma _\text{y}`$, we apply a small tangential force, $`F_\text{T}`$, acting on the left wall (treated as a rigid object with overdamped dynamics) for a certain amount of time, during which the velocity profile, $`u(z)`$, is sampled. The force on the wall is then slightly increased for a new measurement before going over to the next higher value. The static yield stress is then defined as the smallest force (per unit area) for which the average streaming velocity in the left half of the system, $`u_{\text{av, left}}=_{20}^0u(z)𝑑z/20`$, exeeds $`u_{\text{m}in}`$, a minimum value, a few times larger than the typical statistical error of $`u_{\text{av, left}}`$. We empirically find that, for a measurement time of $`4.10^3`$ \[$`5.10^4`$MD steps\], $`u_{\text{m}in}=4.10^4`$ is a reasonable choice. The result on $`\sigma _\text{y}`$ is then averaged over $`10`$ independent runs. Again, for each initial configuration, the system is first equilibrated at $`T=0.5`$ before being quenched to a temperature below $`T_\text{c}=0.435`$. This ensures the equivalence of all individual runs (no history dependence). The system is then propagated with $`F_\text{T}=0`$ for a time $`t_\text{w}`$ before the very first increment of the tangential force. At $`T=0.2`$ and using an increment of $`dF_\text{T}=0.02`$ (once in $`5.10^4`$MD steps) we have obtained $`\sigma _\text{y}=0.596\pm 0.022,\mathrm{\hspace{0.33em}0.658}\pm 0.009`$ and $`0.652\pm 0.015`$ corresponding to waiting times of $`t_\text{w}=10^3,4.10^3`$ and $`4.10^4`$ respectively. Thus, already at a waiting time of $`t_\text{w}=4.10^3`$, aging effects on $`\sigma _\text{y}`$ are negligible at the temperature studied. Note that, due to aging effects, $`\sigma _\text{y}`$ might also depend on the speed at which the threshold value of the tangential force is reached. Therefore, we have determined $`\sigma _\text{y}`$ for a higher value of $`dF_\text{T}=0.05`$ in the case of $`T=0.2`$ and $`t_\text{w}=4.10^4`$. This gives $`\sigma _\text{y}=0.66\pm 0.015`$ which agrees well with $`0.652\pm 0.015`$ within the error bars.

It can be seen in Fig. 3 that $`\sigma _\text{y}>\sigma (\gamma _{\text{tot}}0)`$. Therefore, shear-banding can be expected in the region limited by the vertical dotted line, which corresponds to $`\sigma (\gamma _{\text{tot}})<\sigma _\text{y}`$. Once the yield stress $`\sigma _\text{y}`$ is added to the flow curve, the shear rate becomes multivalued in a range of shear stress, a situation encountered in several complex fluids olmsted . This is the very origin of the shear-banding we observe. As a consequence, this phenomenon should be generic for many soft glassy materials.

Finally, in the very low shear rate region, we observe a time dependence of the shear stress characteristic of stick-slip behavior, as demonstrated in Fig. 4. This is reminiscent to what is observed in friction simulations robbins . Due to numerical limitations, the limit between continuum sliding and stick-slip behavior could not be precisely located. Qualitatively, this behavior is obtained when the thickness of the sheared layer $`h`$ becomes of the order of a few particle diameters, which also corresponds to the width of the interface separating the sheared and jammed regions.

Our results lead to the following picture. (i) For global shear rates $`\gamma _{\text{tot}}`$ smaller than a critical value, $`\gamma _c=U_c/L_z`$, the system separates spatially into two regions, one with a finite, approximatively uniform, shear rate, the other being jammed; (ii) when the total shear rate is increased, the thickness of the sheared region increases; (iii) for $`\gamma _{\text{tot}}>\gamma _c`$, $`h`$ reaches the thickness of the slab, so that the flow is homogeneous with a linear velocity profile; (iv) at very small total shear rates, a stick-slip phenomenon is obtained; (v) the flow profile closely follows the local dynamics of the material, the strongly sheared region corresponding to a rapid structural relaxation, the jammed region to a glassy solid.

The qualitative agreement of these phenomena with experimental observations in very different systems is remarkable: Ref. Pignon::JRheo40::1996  describes explicitly the same scenario for the flow behavior of a Laponite clay suspension, including the stick-slip at very small shear rates. A more quantitative study of the flow profiles has been performed in Ref. Coussot-Raynaud-et-al::PRL88::2002  using a Magnetic Resonance Imaging technique, confirming the existence of well defined jammed and flow regions, whose relative width depends on the global shear rate. In contrast to molecular simulations, an experimental determination of velocity profiles together with a local probe of the dynamics is difficult, which makes the observation of such effects in simulations particularly promising.

Our results do also constrain possible phenomenological descriptions of the flow scenario. It is indeed important to remark that models relating the local shear stress to the local shear rate, as $`\sigma =\sigma _c+\alpha \gamma ^n`$ ($`\alpha `$ being a constant and $`n=1`$ corresponding to Bingham fluids, and $`n>1`$ to Herschel-Buckley fluids), are unable to account for these results. As emphasized above, the shear stress is constant throughout the Couette cell, so that any of the afore-mentioned models would predict a constant shear rate in the cell, in contrast to the present results. It appears therefore necessary to include explicitly non-local terms in phenomenological descriptions olmsted ; Dhont . The existence of a close connection between velocity profile and local dynamics supports some of the assumptions of recent approaches treating the local ‘fluidity’ of the system as an ‘order parameter’ Picard . More detailed simulations, with systematic investigation of size effects, should allow a direct comparison with the predictions of such models.

F.V. thanks the Deutsche Forschungsgemeinschaft (DFG) and the CECAM for financial support. L.B. is supported by the EPSRC. Generous grants of simulation time by the computer centers at the University of Mainz (ZDV) and at “l’Ecole Normale Supérieure de Lyon” (PSMN) are also acknowledged.


text:
# Surface properties at the Kosterlitz-Thouless transition

## 1 Introduction

The two-dimensional classical $`XY`$ model exhibits a non standard behaviour. Long-range order of spins with continuous symmetries is indeed forbidden according to the Mermin-Wagner theorem , and thus there is no spontaneous magnetisation at finite temperature. The transition is governed by unbinding of topological defects. The low-temperature regime was investigated by Berezinskii  and the mechanism of unbinding of vortices was studied by Kosterlitz and Thouless  using renormalization group methods. For reviews, see e.g. Refs. .

The notations are specified below: we consider a $`2d`$ $`XY`$-model on a square lattice with two-components spin variables $`\stackrel{}{\varphi }_w=(\mathrm{cos}\theta _w,\mathrm{sin}\theta _w)`$ located at the sites $`w`$ of a lattice $`\mathrm{\Lambda }`$ of linear extent $`L`$. The spins interact through the usual nearest-neighbour ferromagnetic interaction

$$\frac{H}{k_BT}=K\underset{w}{}\underset{\mu }{}\stackrel{}{\varphi }_w\stackrel{}{\varphi }_{w+\widehat{\mu }},$$
(1)
where $`K=J/k_BT`$ and $`\widehat{\mu }`$ is a unit vector in the $`\mu `$direction.

At low temperature, the system is partially ordered, apart from the existence of vortices, which appear in pairs in increasing number with increasing temperature. In the low-temperature limit, the effect of vortices can be neglected and the behaviour is governed by spin-wave excitations, obtained after expanding the cosine. This harmonic approximation is justified, provided that the spin disorientation remains small, i.e. at sufficiently low temperature. Within this approximation, the quadratic energy in the Boltzmann factor leads to a Gaussian equilibrium distribution and the two-point correlation function becomes

$$\stackrel{}{\varphi }_{w_1}\stackrel{}{\varphi }_{w_2}|w_1w_2|^{1/2\pi K},$$
(2)
hence

$$\eta (T)=k_BT/2\pi J.$$
(3)
When the temperature increases, the influence of vortices becomes more prominent, producing a deviation from the linear spin-wave contribution in equation (3), but the order parameter correlation function still decays algebraically with an exponent $`\eta (T)`$ which depends on the temperature. The existence of such a scale-invariant power-law decay of the correlation function is the signature of a continuous line of fixed points at low temperatures. In this paper we are interested in the end point of this critical line where the topological transition takes place: at the transition temperature $`T_{\mathrm{KT}}`$ (usually called Berezinskii-Kosterlitz-Thouless critical temperature), the pairs are broken and the system becomes completely disordered. This very peculiar topological transition is characterised by essential singularities when approaching the critical point from the high temperature phase, and at $`T_{\mathrm{KT}}`$ the correlation function exponent takes the value $`\eta (T_{\mathrm{KT}})=1/4`$ . In order to recover this result, many numerical studies were performed in the last decade using Monte Carlo simulations (see e.g. Refs. ), but the analysis was made difficult by the existence of logarithmic corrections, e.g.

$$\stackrel{}{\varphi }_{w_1}\stackrel{}{\varphi }_{w_2}\underset{L\mathrm{}}{}|w_{12}|^\eta (\mathrm{ln}|w_{12}|)^{\eta /2},$$
(4)
at $`T_{\mathrm{KT}}`$ . Due to these logarithmic corrections which make the fits quite difficult, the value of $`\eta `$ at $`T_{\mathrm{KT}}`$ was a bit controversial as shown in table 1 of reference . The resort to large-scale simulations, up to systems as large as $`L=1200`$, was then needed in order to confirm this picture .

Recently, we proposed a rather different approach  which is easily implemented. We use the covariance law of correlation functions under the mapping of a two-dimensional system confined inside a square onto the half-infinite plane. The scaling dimensions are then obtained through a simple power-law fit when the correlators are expressed as functions of conveniently rescaled variables. The shape effects are explicitly encoded in the conformal mapping. This is the crucial point, since then even very small systems lead to promising results. We may then determine accurately the bulk correlation function exponent $`\eta `$ at $`T_{\mathrm{KT}}`$ but also, by choosing convenient boundary conditions, the surface critical exponent which describes the decay of the correlation function parallel to a free surface, $`\eta _{}`$.

## 2 Critical system confined inside a square with open boundaries

For a scale-invariant system at its critical point (which also exhibits the properties of isotropic scaling and short-range interactions), conformal invariance enable to include explicitly the shape dependence in the functional expression of the correlators through the conformal covariance transformation under a mapping $`w(z)`$:

$$\stackrel{}{\varphi }_{w_1}\stackrel{}{\varphi }_{w_2}=|w^{}(z_1)|^{x_\sigma }|w^{}(z_2)|^{x_\sigma }\stackrel{}{\varphi }_{z_1}\stackrel{}{\varphi }_{z_2}.$$
(5)
The functional expression of the correlation function inside the square geometry $`w`$ simply follows the mapping $`w(z)`$ which realizes the conformal transformation of the half-plane $`z=x+iy`$ ($`0y<\mathrm{}`$) inside a square $`w=u+iv`$ of size $`L\times L`$ ($`L/2uL/2`$, $`0vL`$) with free boundary conditions along the four edges. This is realized by a Schwarz-Christoffel transformation

$$w(z)=\frac{L}{2\mathrm{K}}\mathrm{F}(z,k),z=\mathrm{sn}\left(\frac{2\mathrm{K}w}{L}\right).$$
(6)
Here, $`F(z,k)`$ is the elliptic integral of the first kind, $`\mathrm{sn}(2\mathrm{K}w/L)`$ the Jacobian elliptic sine, $`\mathrm{K}=K(k)`$ the complete elliptic integral of the first kind, and the modulus $`k`$ depends on the aspect ratio of $`\mathrm{\Lambda }`$ and is here solution of $`K(k)/K(\sqrt{1k^2})=\frac{1}{2}`$. Using the mapping (6), one obtains the local rescaling factor in equation (5), $`w^{}(z)=\frac{L}{2\mathrm{K}}[(1z^2)(1k^2z^2)]^{1/2}`$. One simply has now to include in equation (5) the expected behaviour of $`\stackrel{}{\varphi }_{z_1}\stackrel{}{\varphi }_{z_2}`$ in the half-infinite geometry in order to get the functional expression of the corresponding quantity in the confined system. For example in the semi-infinite geometry $`z=x+iy`$ (the free surface being defined by the $`x`$ axis), the two-point correlation function is fixed up to an unknown scaling function: fixing one point $`z_1`$ close to the free surface ($`z_1=i`$) of the half-infinite plane, and leaving the second point $`z_2`$ explore the rest of the lattice, the following behaviour is expected:

$$\stackrel{}{\varphi }_{z_1}\stackrel{}{\varphi }_{z_2}_{\frac{1}{2}\mathrm{}}(y_1y_2)^{x_\sigma }\psi (\omega ),$$
(7)
with $`\eta =2x_\sigma `$. The dependence on $`\omega =\frac{y_1y_2}{z_1z_2^2}`$ of the universal scaling function $`\psi `$ is constrained by the special conformal transformation , and its asymptotic behaviour is implied by scaling e.g. $`\psi (\omega )\omega ^{x_\sigma ^1}`$ when $`y_21`$, with $`x_\sigma ^1={\scriptscriptstyle \frac{1}{2}}\eta _{}`$ the surface scaling dimension. Keeping $`w_1`$ fixed inside the square, the two-point correlation function becomes

$`\stackrel{}{\varphi }_{w_1}\stackrel{}{\varphi }_w_f[\kappa (w)]^{x_\sigma }\psi (\omega )`$
$`\kappa (w)=\mathrm{Im}\left[\mathrm{sn}{\displaystyle \frac{2Kw}{L}}\right]\times \left|\left(1\mathrm{sn}^2{\displaystyle \frac{2Kw}{L}}\right)\left(1k^2\mathrm{sn}^2{\displaystyle \frac{2Kw}{L}}\right)\right|^{1/2},`$ (8)
where $`f`$ specifies that the boundary conditions (BCs) are chosen free along the edges of the square. This expression is correct up to a constant amplitude determined by $`\kappa (w_1)`$ which is kept fixed, but the function $`\psi (\omega )`$ is still varying with the location of the second point, $`w`$. This makes this expression not so useful in practice (see e.g. reference  for an application). In order to cancel the role of the unknown scaling function, it is more convenient to work with a density profile $`m(w)`$ in the presence of symmetry breaking surface fields $`\stackrel{}{h}_\mathrm{\Lambda }`$ on some parts of the boundary $`\mathrm{\Lambda }`$ of the lattice $`\mathrm{\Lambda }`$. In the case of fixed-free BCs denoted by $`+f`$ (the spins located on the $`x>0`$ half axis are kept fixed while those of the $`x<0`$ are free), Burkhardt and Xue have shown in the case of Ising and Potts models that in the half-plane the order parameter profile obeys the following expression :

$$m_{+f}(z)=y^{\eta /2}(\mathrm{cos}{\scriptscriptstyle \frac{1}{2}}\theta )^{\eta _{}/2}.$$
(9)
In the square geometry, the fixed-free BCs correspond to keeping the spins fixed for $`u>0`$, and the corresponding profile is expected to obey the following ansatz:

$`m_{+f}(w)\left({\displaystyle \frac{L}{2\mathrm{K}}}\right)^{x_\sigma }[\kappa (w)]^{x_\sigma }[\mu (w)]^{x_\sigma ^1},`$
$`\mu (w)={\displaystyle \frac{1}{\sqrt{2}}}\left(1+\mathrm{Re}\left[\mathrm{sn}{\displaystyle \frac{2Kw}{L}}\right]\times \left|\mathrm{sn}{\displaystyle \frac{2Kw}{L}}\right|^1\right)^{1/2}`$ (10)
where $`\kappa (w)`$ was defined in equation (8).

In this paper we will essentially apply equation (10) in order to determine the critical exponents $`\eta `$ and $`\eta _{}`$ at the KT point. In Ref. , the bulk exponent was deduced from a similar analysis of the order parameter profiles in a square-shaped system with fully fixed BCs.

## 3 Monte Carlo simulations

The application of the simple power-law in equation (10) requires a relatively precise numerical determination of the order parameter profile of the $`2d`$ $`XY`$ model confined inside a square with mixed fixed-free boundary conditions. Fixing the spins along some part $`\mathrm{\Lambda }`$ of the boundary, e.g. $`\stackrel{}{\varphi }_w=(1,0),w\mathrm{\Lambda }`$, plays the role of ordering surface fields $`\stackrel{}{h}_{\mathrm{\Lambda }(w)}`$. The resort to cluster update algorithms (Wolff algorithm here) is necessary in order to prevent the critical slowing down: the central idea is the identification of clusters of sites using bond variables connected to the spin configuration. The spins of the clusters are then updated and the algorithm is particularly efficient if the percolation threshold of the bond process occurs at the transition temperature of the spin model. This ensures the updating of clusters of all sizes in a single MC sweep. For the $`XY`$ model, Ising variables defined in the Wolff algorithm by the sign of the projection of the spin variables along some random direction. The bonds are then introduced through the Kasteleyn-Fortuin random graph representation. The percolation threshold for these bonds coincides with the Kosterlitz-Thouless point , which guarantees the efficiency of the Wolff cluster updating scheme  at $`T_{\mathrm{KT}}`$. When one uses particular boundary conditions, with fixed spins along some surface as required here, the Wolff algorithm should become less efficient, since close to criticality the unique cluster will often reach the boundary and no update is made in this case. To prevent this, we use the symmetry properties of the Hamiltonian (1): even when the cluster reaches the fixed boundaries, it is updated, and the order parameter profile is then measured with respect to the new direction of the boundary spins, $`m(w)=\stackrel{}{\varphi }_w\stackrel{}{\varphi }_\mathrm{\Lambda }_{\mathrm{sq}.}`$. The new configuration reached has the right statistical weight.

Using this procedure, we discard $`10^6`$ sweeps for thermalization, and the measurements are performed on $`10^6`$ production sweeps. The simulations at $`T_{\mathrm{KT}}`$ on a typical system of size $`100\times 100`$ takes a few hours on a standard PC with 733 MHz processor, to get a smooth profile as shown in figure 1.

Figure 1: MC simulations of the $`2d`$ $`XY`$ model inside a square $`48\times 48`$ spins with mixed fixed-free boundary conditions at the KT transition temperature $`k_BT_{\mathrm{KT}}/J=0.893`$ (average over $`10^6`$ MCS/spin after cancellation of $`10^6`$ for thermalization).

## 4 Scaling of the order parameter profile $`m_{+f}(w)`$

In Ref. , we have shown that the logarithmic corrections at the KT point are negligible for the order parameter profile with fixed BCs. Assuming that there will be no substantial difference for fixed-free BCs, we perform a fit of $`\mathrm{ln}m_{+f}(w)`$ against the two-dimensional linear expression, $`\mathrm{const}x_\sigma \mathrm{ln}\kappa (w)+x_\sigma ^1\mathrm{ln}\mu (w)`$, as expected from equation (10). The value obtained for the bulk correlation function exponent $`\eta `$ will be a test for the absence of logarithmic corrections. The simulations are performed at the value $`k_BT_{\mathrm{KT}}/J0.893`$ found in the literature, e.g. in Ref. . The fit leads to accurate determinations of both exponents for different system sizes ranging between $`L=20`$ and $`L=200`$. Although it would easily be possible to produce simulations for larger systems, there is no need here, since the results are already quite stable as shown in table 1 and in figure 2.

Table 1:
Bulk and surface scaling dimensions of the order parameter at the Kosterlitz-Thouless transition. The values that we quote are those which correspond to a fit over all the lattice points (except those exactly at the boundaries) ignoring possible lattice effects.    $`L`$ $`x_\sigma `$ $`x_\sigma ^1`$ $`\chi ^2/\mathrm{d}.\mathrm{o}.\mathrm{f}`$    20 0.149(33) 0.290(15) $`0.691\times 10^3`$  24 0.147(34) 0.287(15) $`0.719\times 10^3`$  28 0.134(22) 0.280(12) $`1.210\times 10^3`$  36 0.133(35) 0.278(12) $`1.010\times 10^3`$  48 0.133(35) 0.279(10) $`0.743\times 10^3`$  64 0.131(35) 0.278(10) $`0.664\times 10^3`$  100 0.126(17) 0.275(14) $`0.544\times 10^3`$  200 0.125(23) 0.271(11) $`0.392\times 10^3`$

This is fully coherent with the result of Kosterlitz and Thouless RG analysis $`\eta =2x_\sigma =1/4`$. We can thus also consider as reliable the value of the surface exponent $`\eta _{}=2x_\sigma ^10.54`$.

Figure 2: Size dependence of the bulk and surface scaling dimensions at the KT transition plotted as a function of the inverse linear size of the lattice.

## 5 Scaling of the correlation function $`\stackrel{}{\varphi }_{w_1}\stackrel{}{\varphi }_w_f`$ with free boundary conditions

Equation (8) provides a test in order to control the quality of the surface exponent obtained above. Indeed, if one plots $`\stackrel{}{\varphi }_{w_1}\stackrel{}{\varphi }_w_{\mathrm{sq}.}\times [\kappa (w)]^{x_\sigma }`$ vs $`\omega `$, one should obtain the universal scaling function $`\psi (\omega )`$ whose asymptotic large distance behaviour should lead to a consistent surface exponent through the expected power-law behaviour $`\omega ^{x_\sigma ^1}`$.

Figure 3: Universal scaling function of the surface-volume correlation function for $`L=100`$. The power-law behaviour $`\psi (\omega )\omega ^{x_\sigma ^1}`$ shown in dotted line agrees with a value $`x_\sigma ^10.277`$. The scaling function $`\psi (\omega )`$ is shown in figure 3. It exhibits a clear power-law behaviour consistent with the previous value of the surface scaling dimension.

We mention here that the data were obtained for both quantities, $`\stackrel{}{\varphi }_{w_1}\stackrel{}{\varphi }_w_f`$ and $`m_{+f}(w)`$, with completely independent runs, since the boundary conditions are different in equations (8) and (10) and the quantities themselves are different (2-point or 1-point correlators). The data collapse is quite good and supports the validity of equation (8).

## 6 Conclusion

In this paper we have performed standard MC simulations of the $`2d`$ $`XY`$ model at the Kosterlitz-Thouless topological transition. Combined to the use of a conformal mapping, they provide a simple numerical determination of bulk and surface critical exponents associated to the algebraic decay of the correlation function. The well known result of Kosterlitz and Thouless for the bulk exponent is recovered and as far as the author knows this is the first determination of the surface exponent $`\eta _{}0.54`$. On the contrary to the bulk case, $`\eta _{}`$ is different from the Ising value (which is equal to $`\eta _{}^{Ising}=1`$.