text:
# Radiation Pressure Acceleration by X-rays in Active Galactic Nuclei

## 1 Introduction

Highly ionized gas (HIG) is common in both types of active galactic nuclei (AGN). It is seen as strong absorption features in many type-I (Seyfert 1) galaxies (Reynolds 1997; George et al. 1998, Kaspi et al. 2000 and references therein) and as strong, large equivalent width emission lines in many type-II (Seyfert 2) AGN (e.g. Turner et al. 1997; Netzer, Turner & George 1998). The situation regarding the high luminosity AGN (the bright quasars, see George et al. 2000 for review and references) is still unclear, mostly because of the limited signal-to-noise ASCA observations used in such studies.

The strongest observed HIG features in type-I AGN are several absorption edges, at around 1 keV, mostly due to O vii (0.74keV) & O viii (0.87keV). The deduced hydrogen column density is in the range of $`10^{21}10^{24}`$ cm<sup>-2</sup>, the covering fraction is close to unity and the ionization parameter about two orders of magnitude larger than the broad line region (BLR) ionization parameter (George et al. 1998). Recent Chandra observations reveal the presence of numerous absorption and emission lines (e.g., Kaspi et al. 2000). The absorption lines are narrow ($`150\mathrm{km}\mathrm{s}^1`$)and exhibit outflow velocities of $`<1000\mathrm{km}\mathrm{s}^1`$. Currently, there is no clear indication regarding the location and hence the mass of this component. The very few observations of variable absorption features (e.g., Guainazzi et al. 1996; and references therein) suggest a dimension not too different from that of the BLR. If correct, this would mean that the HIG observed in type-II AGN is of different origin since its distance from the centre must be larger than $`10`$ pc.

Many type-I AGN show also narrow UV absorption features due to C iv$`\lambda 1548`$ and other resonance lines. In most cases the lines are blueshifted with respect to the systemic velocity, implying outflow with typical velocities of $`1000\mathrm{km}\mathrm{s}^1`$. In several well studied cases the lines are resolved, showing intrinsic width of $`100300\mathrm{km}\mathrm{s}^1`$ (Crenshaw et al. 1999; Crenshaw & Kraemer 1999; Srianand 2000; and references therein). There have been several suggestions linking the origin and dynamics of the UV and X-ray absorption features (Mathur et al. 1994; Mathur et al. 1995; Murray & Chiang 1995). This issue is still under discussion (e.g. the new observations of 3C288.1 by Hamann, Netzer & Shields 2000).

Radiation pressure driven flows are thought to be common in astrophysical environments. Stellar winds can be accelerated to their observed velocities ($`>1000\mathrm{km}\mathrm{s}^1`$) by radiation pressure force (e.g., Castor, Abbott & Klein 1975). Broad absorption line (BAL) flows in AGN (BALQSO) can be driven to sub-relativistic velocities ($`30,000\mathrm{km}\mathrm{s}^1`$) by the same physical mechanism (e.g., Arav, Li & Begelman 1994; Murray et al. 1995; Proga, Stone & Kallman 2000) although other mechanisms have also been proposed (Begelman, de Kool & Sikora 1991). In both cases, the absorption line widths are of the order of the outflow velocity which lead to a continuous flow geometry. In this case, the physics can be described by the Sobolev approximation which is valid provided the velocity changes on length scales over which the other gas properties (e.g., ionization state, temperature, density) do not change considerably. Obviously not all types of flow can be driven by radiation pressure forces, among these are jets and supernovae envelopes.

A second, almost orthogonal approach to the study of flows is to consider moving “clouds” or condensations (hereafter, the “cloud model”). Such clouds are confined by means of external pressure and are likely to be in hydrostatic equilibrium. These models have been studied thoroughly in the context of broad line region (BLR) dynamics (e.g., Blumenthal & Mathews 1979; Mathews 1986). The clouds are assumed to be accelerated ballistically, with no appreciable inner velocity gradients and their spectra consists of narrow, blueshifted absorption lines (e.g., Mathews 1975). Cloud confinement and stability have been studied, extensively, in many works (e.g., Krolik 1979; Mathews 1986; Mathews 1992). Possible confinement mechanisms include thermal confinement (e.g., Krolik, McKee & Tarter 1981) with its well known problems (e.g., Mathews & Ferland 1987), and nonthermal confinement, e.g., by magnetic fields (Rees 1987; Bottorff, Korista & Shlosman 2000). Despite the many unresolved issues, the cloud model has had great success in accounting for the observed features of the BLR clouds (e.g. Netzer 1990), NLR clouds (e.g., Kraemer, Ruiz & Crenshaw 1998), and BAL flows (e.g. Arav et al. 1994).

While acceleration by UV radiation, in AGN, has been studied extensively, there are only very few works discussing the acceleration by X-ray radiation (e.g., Mathews & Capriotti 1985; Mathews & Veilleux 1989, none in great detail). In particular, there is no systematic study of the acceleration of gas with properties similar to those commonly observed in the X-ray spectrum of many AGN (for a heuristic treatment of the problem see Reynolds & Fabian 1995, for a magneto-hydrodynamical model with comparison to observations see Bottorff, Korista & Shlosman 2000). As shown below, considerations of the relevant parameter space within the framework of the cloud model, result in very different dynamics for the HIG clouds compared with previous works.

This paper investigates the acceleration of the HIG clouds. The purpose is to study the micro-physics of the gas and the radiation pressure force, and to derive the gas velocity under various conditions prevailing in type I AGN. This will enable to explain observed spectral features in new, high resolution X-ray spectra.

The paper is organized as follows: in §2 we outline the main ingredients of the model. Section 3 describes the calculation procedure, and present the results concerning the cloud structure, the radiation pressure force and the gas motion, with a following discussion on the model limitations. The main conclusions are summarized in §4.

## 2 The model

### 2.1 Basic assumptions

Our purpose is to solve for the structure and the dynamics of HIG clouds in AGN. Such gas has been detected by ASCA and earlier missions but the evidence used here is found mainly in the recent Chandra observations such as the one reported in Kaspi et al. (2000) for NGC 3783 and in Kaastra et al (2000) for NGC 5548. In both cases there is a rich absorption spectrum indicating high level of ionization, large column density ($`10^{22}\mathrm{cm}^2`$) and outflow motion of $`<1000\mathrm{km}\mathrm{s}^1`$. The lines are barely resolved (Kaspi et al. 2001) and the typical Doppler width is ($`300\mathrm{km}\mathrm{s}^1`$). The observed properties are, therefore, different from those found in stellar winds and in BAL quasars. They can be interpreted as either due to cloud motion or continuous wind. This paper addresses the first possibility and the wind solution will be discussed in a future publication.

The main assumptions regarding the cloud model are (see also §3.5):

1. The general AGN environment can support various type of clouds. Known examples are the BLR and the NLR clouds (e.g. Netzer 1990; see a recent study by Kaiser et al. 2000). Such entities with higher level of ionization cannot be ruled out by any known AGN property and have, in fact, been addressed, by several recent publications (e.g., Contini & Viegas 1999; Bottorff et al. 2000). The HIG region is modeled as a spherical cloud system of large covering factor. The justification is that such regions are very common in low luminosity AGN (e.g. George et al. 1998). We concentrate on thin spherical shells although this is not a requirement of the model (see §3.5.)
2. We assume that the clouds are confined by external pressure $`P_{\mathrm{ext}}`$. This is required since clouds which are not self-gravitating evaporate on sound crossing time scales (which is typically much shorter than dynamical time scales). While the origin of this pressure is not the topic of the paper, we note that confinement by a hot gas (e.g. Krolik, McKee & Tarter 1981) is problematic (Mathews & Ferland 1987) and magnetic pressure is the most likely explanation (see Rees (1987) for a general formulation, de Kool & Begelman (1995) for a discussion of radiatively driven, magnetically confined flows and Sivron & Tsuruta (1993) for a discussion of the magnetically confined cloud model). The pressure is assumed to be time-independent and its radial dependence is given by

$$P_{\mathrm{ext}}R^\alpha ,$$
(1)
Typically, $`10/3\alpha 3/2`$ (e.g., Netzer 1990). The global spherical symmetry is not critical to our discussion and is relaxed in §3.5.
3. We restrict the discussion to geometrically thin HIG clouds, $`\mathrm{\Delta }R/R1`$ where $`R`$ is the location of the illuminated surface (see for example Reynolds & Fabian 1995). This is based on variability studies which set rough lower limits on the density ($`10^6\mathrm{cm}^3`$, see Otani et al. 1996; Guainazzi et al. 1996) and thus, through the known ionization parameter (see below) on the location of the clouds. The external pressure difference $`\mathrm{\Delta }P_{\mathrm{ext}}P_{\mathrm{ext}}(R)P_{\mathrm{ext}}(R+\mathrm{\Delta }R)`$ satisfies $`\mathrm{\Delta }P_{\mathrm{ext}}/P_{\mathrm{ext}}=\alpha \mathrm{\Delta }R/R1`$. This assumption is checked for self-consistency in §3.5.
4. The clouds are assumed to be in photo-ionization equilibrium. This is justified as long as the recombination time is the shortest time scale. This is true in our parameter range provided the ionizing continuum is not rapidly varying (Krolik & Kriss 1995).
5. The dynamical problem is treated by assuming radiation pressure force, gas pressure force, and gravity. The local gas acceleration, $`a(r)`$, at position $`r`$ inside a cloud whose illuminated face is at a distance $`R`$ from the centre, is given by

$$a(r)=a_{\mathrm{rad}}(r)g\frac{1}{\rho }\frac{dP_{\mathrm{gas}}}{dr}$$
(2)
where $`a_{\mathrm{rad}}`$ is radiation pressure acceleration pressure, $`g`$ the gravitational acceleration, $`\rho `$ the gas density, and $`P_{\mathrm{gas}}`$ the gas pressure. Note that the dependence of $`g`$ and $`a_{\mathrm{rad}}`$ on $`R`$ has been omitted since $`R`$ is constant across the geometrically thin cloud. Furthermore, we assume quasi-hydrostatic cloud structure in the cloud’s rest frame. As shown by Blumenthal & Mathews (1979), and others, this is an adequate description of BLR and similar clouds. The clouds considered here are illuminated by the same radiation field and are roughly at the same location, hence the same general arguments apply. The clouds are therefore regarded as moving with a certain, distance dependent, bulk velocity and acceleration ($`a`$, i.e., $`a(r)=a`$), with no significant velocity gradients. Given this, we do not have to treat the full hydrodynamic problem but assume that the clouds remain stable and do not disperse as they are accelerated. The quasi-hydrostatic assumption is checked for self-consistency in §3.5.

### 2.2 Ionization structure and radiation pressure force

The ionization structure of the cloud is determined by the local radiation field. This is usually described in terms of the illuminated face ($`r=0`$) ionization parameter, $`U`$,

$$U=\frac{_{\nu _0}^{\nu _1}\left(L_\nu /h\nu \right)𝑑\nu }{4\pi R^2n_Hc}$$
(3)
where $`L_\nu `$ is the monochromatic continuum luminosity ($`\mathrm{erg}\mathrm{s}^1\mathrm{Hz}^1`$), $`n_H`$ the hydrogen number density at $`r=0`$ and $`c`$ the speed of light. Our aim is to study the acceleration of highly ionized clouds. We therefore use the X-ray ionization parameter, $`U_x`$ (Netzer 1996), defined by $`h\nu _0=0.1`$ keV and $`h\nu _1=10`$ keV. The ionization and thermal structure are calculated by ion99, the 1999 version of the photo-ionization code ion (Netzer 1996, and references therein). The code is very detailed regarding X-ray related processes and solves for the temperature profile and the radiation pressure force distribution inside the cloud. All our calculations apply to fully exposed clouds (i.e., no obscuration or shielding).

The absorption of continuum radiation by bound-bound, bound-free, Compton scattering and free-free processes results in a net momentum transfer to the clouds. The relative importance of the various processes depends on the atomic cross sections and the spectral energy distribution (SED). The largest cross sections are for bound-bound transitions. These occur at energies lower than the ionization energy of the relevant ion where, for the assumed AGN-type SED (see below) the value of $`L_\nu `$ is larger than its value at the ionization threshold. However, large cross sections imply large optical depths and self shielding. Thus, bound-bound absorption contributes the most at the illuminated surface and decays quickly into the cloud (Fig. 1). Bound-free processes are the next largest contributers. The Compton cross-section is relatively small, and the optical depth is negligible in all cases considered here (column density $`<1.5\times 10^{24}\mathrm{cm}^2`$). Thus, Compton scattering contributes equally (per particle) across the cloud. In almost all cases (see however §3.1), free-free absorption contributes very little to the radiation pressure force ($`0.01`$ Compton). We neglect internal line radiation pressure (see Elitzur & Ferland 1986) which do not contribute significantly to the total pressure because of the low excitation temperature of the strongest X-ray transitions.

We use the force multiplier formalism (e.g., Mihalas & Mihalas 1999). The force multiplier is defined as the ratio between the radiation pressure force due to all processes relative to the force due to Compton scattering,

$$M(r)\frac{_0^{\mathrm{}}\chi _\nu (r)L_\nu e^{\tau _\nu (r)}𝑑\nu }{n_e(r)\sigma _TL_{\mathrm{tot}}},$$
(4)
where $`\chi _\nu (r)`$ is the frequency dependent total absorption and scattering coefficient, $`n_e`$ is the number density of free electrons, $`L_{\mathrm{tot}}`$ is the total, unattenuated ($`r=0`$) luminosity, and $`\tau _\nu (r)`$ is the position dependent optical depth. This gives

$$a_{\mathrm{rad}}(r)=\frac{n_e(r)\sigma _TL_{\mathrm{tot}}}{4\pi R^2\rho (r)c}M(r).$$
(5)
### 2.3 Model parameters

Several parameters influence the solution. The SED detrmines the ionization structure and the radiation pressure force over the various energy bands. We have experimented with several SEDs typical of type-I AGN and found only little differences regarding the structure and dynamics of the clouds. We therefore take as our standard and only case an SED characterized by a “UV bump”,$`\alpha _{ox}=1.4`$, $`\alpha _{(0.51keV)}=1.06`$, and $`\alpha _{(110keV)}=0.9`$.

The internal line widths determine the line optical depth and thus the bound-bound contributions to $`M`$. We have experimented with three possibilities: a) Thermal line width, $`v_T`$. b) Gaussian-shaped profiles with line widths for all ions equal to the thermal hydrogen line width (denoted by $`v_{\mathrm{sound}}`$). c) Gaussian-shaped, line profiles, characterized by their FWHM (typically several hundred $`\mathrm{km}\mathrm{s}^1`$) and motivated by the Crenshaw et al. (1999) and Kaspi et al. (2001) observations. Mechanisms responsible for such broadening have been discussed, extensively, in the literature (see Bottorff & Ferland (2000) for discussion of MHD turbulences, including many older references, as well as many papers discussing microturbulences and expansion in stellar atmospheres).

The range of ionization parameter and column density is determined by the observations (Reynolds & Fabian 1995; Reynolds 1997; George et al. 1998; Kaspi et al. 2001). We chose $`0.02<U_x<4`$ and $`10^{20}N_H10^{23}\mathrm{cm}^2`$. Many of the results pertain to a “standard model” defined by $`U_x=10^{0.8},n_H=10^8\mathrm{cm}^3,N_H=10^{22}\mathrm{cm}^2`$, with line width given by $`v_T`$. We consider only the following solar-like composition, $`\mathrm{H}:\mathrm{He}:\mathrm{C}:\mathrm{N}:\mathrm{O}:\mathrm{Ne}:`$ $`\mathrm{Mg}:\mathrm{Si}:\mathrm{S}:\mathrm{Ar}:\mathrm{Fe}=10^4:10^3:3.7:1.1:8:1.1:0.35:0.37:0.16:0.03:0.4`$

### 2.4 Method of solution

We follow the general method described in Weymann (1976), Mathews & Blumenthal (1977) and Blumenthal & Mathews (1979) where the cloud structure obeys the hydrostatic equation

$$\frac{dP_{\mathrm{gas}}(r)}{dr}=\rho (r)\left[a_{\mathrm{rad}}(r)ga\right],$$
(6)
where $`a`$ is the average (bulk) acceleration. This is satisfied for every distance $`R`$. The radiation pressure force is changing across the cloud due to optical depth effects. This results in a non-uniform gas pressure profile which must obey the boundary condition $`P_{\mathrm{gas}}(r=0)=P_{\mathrm{ext}}(R)P_{\mathrm{ext}}(R+\mathrm{\Delta }R)=P_{\mathrm{gas}}(r=\mathrm{\Delta }R)`$. The exact distribution depends on the ionization structure, the level population of the various ions, and the temperature at every point inside the cloud.

The cloud accelerates ballistically and we make no use of the continuity equation. Nevertheless, the cloud column density can change due to the external conditions since it can accrete mass, lose mass, or expand along the way. In this work we treat the 1D problem and consider the following $`N_H(R)`$: 1) constant mass gas shells which subtend a constant solid angle, $`\mathrm{\Omega }4\pi `$, i.e. $`\mathrm{\Omega }R^2N_H=\mathrm{const}`$, 2) constant column density clouds. The latter model is included as an example for extreme cases where the column density decreases with distance much slower than $`R^2`$. For example, spherical constant mass clouds imply $`N_HR^{2/3\alpha }`$. This allows us to study the sensitivity of our results to different column densities. In §3.5 we consider several other implications of these assumptions.

To calculate the acceleration $`a`$, we integrate equation (6) with the given boundary conditions (i.e., pressure) at both cloud surfaces. Given our parameter space and set of assumptions, the contribution of the small external pressure gradients to the acceleration via $`\mathrm{\Delta }P_{\mathrm{ext}}/(\rho \mathrm{\Delta }R)`$ is negligible ($`10^4g`$) and the bulk acceleration is given by

$$a=\frac{_0^{\mathrm{\Delta }R}\rho (r)a_{\mathrm{rad}}(r)𝑑r}{_0^{\mathrm{\Delta }R}\rho (r)𝑑r}g,$$
(7)
which depends on $`R`$ via the boundary conditions, the value of $`N_H`$ and via $`a_{\mathrm{rad}}(r)`$.

Keeping with the definition of the average acceleration, we define the average force multiplier,

$$M=\frac{_0^{\mathrm{\Delta }R}n_e(r)M(r)𝑑r}{_0^{\mathrm{\Delta }R}n_e(r)𝑑r},$$
(8)
where we have made use of the fact that $`U_x`$ is large and hydrogen and helium are fully ionized, thus $`\rho (r)/n_e(r)1.18m_H`$. Thus we obtain the hydrostatic equation in the force multiplier representation,

$$\frac{dP_{\mathrm{gas}}(r)}{dr}=\frac{n_e(r)\sigma _TL_{\mathrm{tot}}}{4\pi R^2c}\left[M(r)M\right].$$
(9)
The cloud dynamics is therefore reduced to the ballistic case where the cloud is accelerated coherently with the following equation of motion

$$v\frac{dv}{dR}=a(R)=\frac{1}{1.18m_H}\frac{\sigma _TL_{\mathrm{tot}}}{4\pi cR^2}M(R)\frac{GM_{\mathrm{BH}}}{R^2}$$
(10)
where $`v`$ is the velocity of the cloud and $`M_{\mathrm{BH}}`$ is the mass of the central black hole. Thus, $`M`$ is the only parameter which depends on the inner structure of the cloud.

Defining a new variable, $`l`$, which is proportional to $`1/M_{\mathrm{BH}}`$, by

$$l=\frac{1}{1.18m_H}\frac{\sigma _TL_{\mathrm{tot}}/c}{4\pi GM_{\mathrm{BH}}}=\frac{L_{\mathrm{tot}}}{1.18L_{\mathrm{edd}}},$$
(11)
simplifies the equation of motion to

$$v\frac{dv}{dR}=\frac{1}{1.18m_H}\frac{\sigma _TL_{\mathrm{tot}}}{4\pi cR^2}\left(M(R)\frac{1}{l}\right).$$
(12)
In the above formalism, the dynamical problem is decoupled from the cloud structure problem (i.e. the structure is independent of the cloud velocity).

The actual calculations start by assuming a uniform density profile and by solving for the temperature and radiation pressure force profiles. We calculate $`M`$, solve the hydrostatic equation and obtain the first estimate of $`P_{\mathrm{gas}}(r)`$. Given the first temperature profile, we obtain a new density profile which is used, as an input, for the next iteration. This procedure is repeated several times until convergence. Tests show that the final results depend very little on the initial density profile. The iteration stops when the density profile does not change by more than 1% between one iteration to the next.

Given the pressure and ionization profile, the solution of the equation of motion is straight forward. The initial conditions include the distance from the centre at the base of the flow, $`R_0`$, $`P_{\mathrm{ext}}(R_0)`$, $`N_H(R_0)`$, the column density as a function of $`R`$ (i.e., constant column clouds or constant mass shells), and the AGN parameters ($`\alpha `$, $`L_{\mathrm{tot}}`$ and $`l`$). This allows us to solve for the cloud structure and for $`v(R)`$ at each $`R`$.

## 3 Results and Discussion

This section describes the results of the calculations regarding the cloud’s internal structure, the average force multiplier $`M`$, and the cloud motion. Unless otherwise specified the standard model parameters are assumed.

### 3.1 Cloud structure

We allow the ionization parameter to vary in the range $`10^{1.8}10^{0.6}`$ assuming all other parameters remain constant. The contributions of the various absorption and scattering processes to $`M`$ are computed from the ionization structure. As seen in Fig. 2, the force multiplier decreases with $`U_x`$. The reasons are the fewer bound-bound transitions at large $`U_x`$ and the lower available driving flux at high energies.

A common property of all models is the increasing density into the cloud (Fig. 2), due to the temperature decrease as a function of $`r`$ and the constant external pressure. The most extreme case shown here is for $`U_x=10^{0.8}`$ where the density increases by a factor $`2`$ above its initial value. This is also the ionization parameter where the ionization and thermal structure are dominated by H-like and He-like oxygen.

Fig. 2 also shows the internal pressure profile. The rise near $`r=0`$ is caused by the pressure gradient that must balance the radiation pressure force which is larger at this location due to the contribution of bound-bound processes. The peak pressure depends on the difference $`M(r=0)M`$ and the location $`r`$ where $`M(r)=M`$. For very large values of $`U_x`$, the peak is lower since the gas is highly ionized and bound-bound transitions are less important.

Finally, since the value of the temperature at the illuminated surface, $`T`$, is independent of the column density, we can approximate it by $`T10^{5.5}U_x^{0.5}`$ which is a good approximation for the abundances, SED and the range of $`U_x`$ considered here. The same functional form, with different coefficients, provides a good fit for other SEDs (see, e.g., Komossa & Meerschweinchen 2000).

Changes in column density affect the pressure and temperature profile because of the increased opacity. This is shown in Fig. 3. For low column densities, and the standard $`U_x`$, the clouds have almost uniform ionization structure. The ionization structure and the radiation pressure force distribution change considerably for $`N_H10^{22}\mathrm{cm}^2`$ since this column density corresponds to an optical depth of $`1`$ in the O vii and O viii bound-free transitions.

A common feature for all column densities is the fast decline of the bound-bound force multiplier, $`M_{\mathrm{bb}}`$, beyond $`N_H10^{17}\mathrm{cm}^2`$, due to the large line opacity. Deeper layers contribute much more to the bound-free force multiplier, $`M_{\mathrm{bf}}`$. This behavior is column density dependent since, for large columns, the increased bound-free opacity causes lower energy bands of the continuum to become more significant in producing the local radiation pressure (for the typical SED considered here, there are many more photons in these bands). The increase in $`M_{\mathrm{bf}}`$ is the reason for the peculiar pressure profile in the $`N_H>10^{22}\mathrm{cm}^2`$ models.

We study clouds with a large range of densities ($`1<n_H<10^{14}\mathrm{cm}^3`$) at the illuminated surface. As shown in Fig. 4, the cloud structure is not very sensitive to $`n_H`$ for $`n_H(r=0)10^{11}\mathrm{cm}^3`$ since the micro-physics remains essentially unchanged. All changes can therefore be reduced to a scale change in $`\mathrm{\Delta }R`$. Above $`n_H=10^{11}\mathrm{cm}^3`$, the micro-physics changes considerably due to collisional excitation and free-free heating. This results in more uniform temperature and density profiles and hence changes in $`M_{\mathrm{bb}}`$. In general, the very high density clouds are characterized by higher ionization, higher temperature and smaller $`M`$.

The internal velocity spread is dynamically important. Wider profiles decrease the line opacity and increase the thickness of the zone where $`M(r)`$ is dominated by bound-bound transitions. This shifts the peak of the pressure profile deeper into the cloud (Fig. 5). A large FWHM does little to affect the local ionization and temperature balance, as expected from purely radiative considerations. However, the increased importance of $`M_{\mathrm{bb}}`$ results, in low column density clouds, in a more uniform pressure profile.

### 3.2 The dependence of $`M`$ on $`U_x`$, $`N_H`$, $`n_H`$ and the line width

The average force multiplier is the primary property determining the cloud velocity. It is therefore important to investigate its behavior as a function of the cloud properties.

The dependence of $`M`$ on the ionization parameter, $`U_x`$ is shown in Fig. 6. Since the clouds are more ionized for larger $`U_x`$, $`M`$ is more sensitive to higher frequency bands, where the incident flux is lower. The result is that $`M`$ is a decreasing function of $`U_x`$ for all column densities. For low column densities ($`N_H<10^{22}\mathrm{cm}^2`$) the change in $`M`$ is roughly linear with $`U_x`$. For larger columns, the structure is more complicated (e.g., Fig. 3) and so is the dependence on $`U_x`$. A crude approximation of the form $`\mathrm{log}M(U_x)=C_1+C_2\mathrm{log}(U_x)`$, for several column densities and for line widths given by $`v_{\mathrm{sound}}`$, is

$$\mathrm{log}\left(\begin{array}{c}M_{20}\\ M_{21}\\ M_{22}\\ M_{23}\end{array}\right)=\left(\begin{array}{c}0.7\\ 0.6\\ 0.5\\ 0.7\end{array}\right)\mathrm{log}(U_x)+\left(\begin{array}{c}0.6\\ 0.5\\ 0.4\\ 0.4\end{array}\right),$$
(13)
where $`M_n`$ is the average force multiplier for $`N_H=10^n\mathrm{cm}^2`$. This approximates $`M`$ within a factor of $`2`$.

Fig. 6 also shows the dependence of $`M`$ on the column density. This is characterized by a decrease followed by an increase with a minimum at $`N_H10^{22}\mathrm{cm}^2`$. The minimum is mostly due to the fact that at this column density, the optical depth due to absorption by oxygen is close to unity. For $`N_H10^{22}\mathrm{cm}^2`$, the clouds are ionized throughout (Fig. 3) and the decline in $`M`$ reflects a decrease in the available X-ray flux. For $`N_H>10^{22}\mathrm{cm}^2`$, the ionization structure is such that the leading edge is neutral enough to absorb a large fraction of the far UV radiation. Since the photon flux at those bands is much larger, it results in an increasing $`M`$. The changes in $`M`$ are more pronounced for low values of $`U_x`$ because of the SED shape.

The gas density has only a small influence on the cloud structure, and on $`M`$, for $`n_H(r=0)<10^{11}\mathrm{cm}^3`$ (Fig. 7). For higher densities, the micro-physics changes (see §3.1) and the changes in $`M`$ can amount to a factor $`2`$ between $`n_H(r=0)10^8\mathrm{cm}^3`$ and $`n_H(r=0)10^{14}\mathrm{cm}^3`$. As shown below, for the $`P_{\mathrm{ext}}R^2`$, $`N_H=\mathrm{const}.`$ case, this corresponds to a change of $`40\%`$ in terminal velocity. Such changes can be critical for marginally accelerating clouds.

The dependence of $`M`$ on the FWHM of oxygen lines is shown in Fig. 8 where we consider FWHM$`<16,000\mathrm{km}\mathrm{s}^1`$, similar to BAL flow velocities. We do not show cases with $`v_{\mathrm{sound}}`$ which are roughly equivalent to FWHM$`70\mathrm{km}\mathrm{s}^1`$.

As expected, broader lines produce larger force multipliers ($`6`$ times larger for FWHM$`3000\mathrm{km}\mathrm{s}^1`$ compared with thermal case). The effect is more pronounced for small $`U_x`$ where there are more lines available to absorb the larger continuum flux. For $`U_x=10^{0.2}`$ there is a decrease in $`M`$ for FWHM$`>3000\mathrm{km}\mathrm{s}^1`$ due to changes in ionization structure which mask the relatively small increase in line contributions. For FWHM$`>16,000\mathrm{km}\mathrm{s}^1`$, $`M`$ is relatively constant since the lines are optically thin relative to the bound-free processes.

While the cloud structure and, therefore, the radiation pressure force are relatively insensitive to the exact shape of the SED (i.e., $`\alpha _x`$, $`\alpha _{ox}`$) for the same $`U_x`$, the exact value of $`M`$ depends roughly linearly on the X-ray (0.1-1keV) to total luminosity ratio (equation 4) which in our case is $`0.04`$.

### 3.3 Approximations for $`M`$

Fig. 9a illustrates the most important contributions to $`M`$ with respect to the cloud dynamics. While the exact value depends on the cloud’s structure, the bound-free force multiplier, $`M_{\mathrm{bf}}`$, is, in general, a good approximation to $`M`$. It is possible to estimate $`M_{\mathrm{bf}}`$ by considering oxygen, the most abundant metal. Assuming cross sections of the form $`\sigma =\sigma _0(E/E_0)^s`$, where $`E_0`$ is the ionization energy, $`\sigma _0`$ the edge cross section, $`s2.53`$, and an X-ray spectrum with an energy power-law index $`\alpha _x`$. For an ion $`X`$ with a number density $`n_x`$ we find

$$M_{\mathrm{bf}}=\frac{1}{\alpha _x+s1}\left(\frac{n_X}{n_e}\right)\left(\frac{\sigma _0}{\sigma _T}\right)\left(\frac{L_E(E_0)E_0}{L_{\mathrm{tot}}}\right),$$
(14)
where $`L_E(E_0)`$ is the monochromatic flux at the ionization edge. For the range of interest, $`10^{1.8}U_x10^{0.2}`$, and the chosen SED, the main contributions to the force multiplier are O vii, O viii, and He ii. Applying the relevant constants, and assuming $`X_{\mathrm{O}\mathrm{vii}}+X_{\mathrm{O}\mathrm{viii}}1`$, we obtain the following estimate:

$$M_{\mathrm{bf}}1.6+1.3X_{\mathrm{O}\mathrm{viii}}+2\left(\frac{X_{\mathrm{He}\mathrm{ii}}}{10^4}\right),$$
(15)
where $`10^6<X_{\mathrm{He}\mathrm{ii}}<10^4`$. This $`M_{\mathrm{bf}}`$ is within a factor $`3`$ of $`M`$ for line width given by $`v_T`$. This result as well as the general (all elements) $`M_{\mathrm{bf}}`$ is significantly different from the force multiplier calculated by Mathews & Veilleux (1989) who considered gas with much lower ionization stages and calculated only the contribution due to the carbon edge. For the typical $`U_x`$ and SED considered here this gives a value of $`M`$ which is a factor of 50(!) smaller than the one calculated by us.

### 3.4 Cloud dynamics

Solving the equation of motion (equation 12) requires the value of $`M(R)`$ for every $`R`$. The formal solution, assuming $`v(R_0)=0`$, is

$$v(\xi )=v_c\left[_1^\xi \left(M(\xi ^{})\frac{1}{l}\right)\frac{d\xi ^{}}{\xi ^2}\right]^{1/2},$$
(16)
where $`\xi =R/R_0`$. The characteristic velocity $`v_c`$ is defined as

$$v_c=\sqrt{\frac{1}{1.18m_H}\frac{\sigma _TL_{\mathrm{tot}}}{2\pi cR_0}}1350\left(\frac{L_{45}}{R_{17}}\right)^{1/2}\mathrm{km}\mathrm{s}^1,$$
(17)
where $`L_{45}`$ and $`R_{17}`$ are the luminosity in units of $`10^{45}\mathrm{erg}\mathrm{s}^1`$ and the initial distance ($`R_0`$) in units of $`10^{17}\mathrm{cm}`$. The dynamical problem is thus reduced to the calculation of a single quantity, $`M(R)`$, at each position along the cloud trajectory.

We have studied velocity profiles for two generic cases: shell-like clouds with constant mass ($`\mathrm{\Omega }R^2N_H=\mathrm{const}`$) and constant column density clouds ($`N_H(R)=\mathrm{const}`$). The exact numerical solutions are shown in Fig. 10. The diagram illustrates the lower velocities that are obtained for the constant column density clouds since these are more massive (note that a comparison between the two is not straight forward since $`M`$ is a non-monotonic function of $`N_H`$, e.g., §3.3). For the standard model, the velocities are of the same order as the escape velocity at the base of the flow. In particular, velocities of 1000–3000 $`\mathrm{km}\mathrm{s}^1`$ are typical of HIG clouds originating in or just outside the BLR. Much larger velocities are achieved for higher density clouds that originate inside the BLR (related, perhaps, to the inner accretion disc). For example, the standard model with $`n_H=10^{12}\mathrm{cm}^3`$ and constant column density results in a terminal velocity of $`20,000\mathrm{km}\mathrm{s}^1`$ for line width given by $`v_{\mathrm{sound}}`$. Thus, acceleration by X-ray radiation pressure is capable of producing BAL-type flows in AGN. This conclusion applies only to the ballistic motion of hydrostatically stable clouds. Wind-type flows may be different and will be discussed in future works.

Three cases of external pressure powerlaws were examined for constant column density clouds and for constant mass shells, all with $`n_H<10^{11}\mathrm{cm}^3`$ (see Fig. 10). The main results are explained using the constant column density cloud case.

The first case $`U_x(R)=\mathrm{const}`$ ($`\alpha =2`$). Here $`M(R)`$ is a constant of motion and equation (16) results in

$$v(R)=v_{\mathrm{final}}\sqrt{1\frac{R_0}{R}}.$$
(18)
where the final (asymptotic) velocity is

$$v_{\mathrm{final}}=v_c\left[M(U_x)\frac{1}{l}\right]^{1/2}.$$
(19)
The final velocity of marginal outflows is therefore very sensitive to the value of $`l`$ since $`1/lM`$ (Fig. 9b) The acceleration time-scale, defined as the time it takes the cloud to reach 90% of its asymptotic velocity, is

$$t_{\mathrm{accel}}150\frac{v_c}{v_{\mathrm{final}}}R_{17}\mathrm{years}.$$
(20)
In the two other cases, $`U_x`$ decreases with $`R`$ ($`\alpha >2`$) and $`U_x`$ increases with $`R`$ ($`\alpha <2`$) (see Fig. 10). In the first case, the acceleration due to the radiation pressure force decreases as $`1/\xi `$ to some (positive) power for our range of $`\alpha `$. For $`\xi 1`$, the clouds become asymptotically neutral and the radiation acceleration drops to zero.

In the second case, $`U_x`$ increases with $`R`$, the clouds become more ionized as they move out, bound-bound and bound-free processes become less efficient and at large $`\xi `$, the only driving force is due to Compton scattering.

The condition for asymptotic outflow ($`\xi 1`$) is that the cloud velocity exceeds the escape velocity at some distance,

$$v(R)>v_{\mathrm{esc}}(R)=\sqrt{\frac{1}{1.18m_H}\frac{\sigma _TL_{\mathrm{tot}}}{2\pi cRl}}=v_c\sqrt{\frac{1}{l}}.$$
(21)
Contrary to other wind solutions, in which the final velocity is of the order of $`v_{\mathrm{esc}}(R_0)`$, the final velocity in our model depends on the value of $`\left[M(R)l1\right]`$ (equation 16) which can be much larger than unity.

Fig. 10 also shows that a fair fraction of the final velocity is obtained at small $`\xi `$. This is the reason for the similar terminal velocity in the $`\alpha =2`$ and the $`\alpha =3/2`$ cases. A different behaviour is seen for $`\alpha =10/3`$ where the gravitational force exceeds the radiation pressure force at $`\xi 3`$, resulting in a marginal outflow.

The dynamics of constant mass gas shells follows the same trend for the various external pressure profiles, but with higher terminal velocities. This results from the decrease in column density with distance and the fact that for $`N_H<10^{22}\mathrm{cm}^2`$, $`M`$ is a monotonically decreasing function of $`N_H`$ (see Fig. 6).

So far we have considered radial motion only. Other AGN components, such as the BLR and the accretion disc, are thought to be rotating around the central mass, and a plausible assumption is to associate a rotational motion to the HIG. Introducing Keplerian velocities (see also §3.5) at $`R_0`$, while conserving angular momentum along its trajectory, will decrease the effective gravity (in equation 12, the term $`1/l`$ will be replaced by $`1/l(11/\xi ^{})`$) thus increasing the final velocities (see Fig. 9b).

The obtained dynamical results differ significantly from the few previous calculations involving acceleration of clouds by X-ray radiation pressure force. It is an order of magnitude difference from the Mathews & Capriotti (1985) results (see their table 1) and this is not surprising since previous works did not consider the relevant parameter space for the HIG and, therefore, did not include the relevant radiation absorption mechanisms.

### 3.5 Limitations and extensions of the model

We have made two critical assumptions: thin shell clouds (which is always true, by construction, at $`R_0`$) and hydrostatic equilibrium at all locations. Below we derive general expressions for a critical, normalized distance, $`\xi _c`$, beyond which one or both of these assumptions fail.

The critical distance for the thin-shell approximation is

$$\xi _c^{(\mathrm{thin})}=\{\begin{array}{cc}\left(\frac{R_0}{\mathrm{\Delta }R_0}\right)^{\frac{1}{\beta \alpha 1}}\hfill & \beta \alpha >1\hfill \\ \mathrm{}\hfill & \mathrm{otherwise}\hfill \end{array}$$
(22)
where we have used $`n_HP_{\mathrm{ext}}`$. Here, $`\beta =0,2`$ for constant column density clouds and constant mass shells, respectively. Since the initial conditions are such that $`\mathrm{\Delta }R_0/R_01`$, we see that in general $`\xi _c^{(\mathrm{thin})}1`$.

A quasi-hydrostatic cloud structure can be obtained provided the cloud’s sound crossing time is shorter than the time it takes the external conditions (i.e., external pressure) to change considerably (e.g., Blumenthal & Mathews 1975) as a result of the cloud’s motion. This time is roughly $`R_0/v_{\mathrm{final}}`$. The quasi-hydrostatic approximation is therefore valid provided $`\xi \xi _c^{(\mathrm{hyd})}`$ where

$$\xi _c^{(\mathrm{hyd})}\left(\frac{v_{\mathrm{sound}}}{v_{\mathrm{final}}}\frac{R_0}{\mathrm{\Delta }R_0}\right)^{\frac{1}{\beta \alpha }}.$$
(23)
In cases where $`\beta \alpha 0`$, the cloud structure and dynamics can be treated by our model for all $`\xi 1`$ provided $`R_0v_{\mathrm{sound}}/(\mathrm{\Delta }R_0v_{\mathrm{final}})>1`$. In conclusion, cases for which $`\xi _c=\mathrm{min}(\xi _c^{(\mathrm{thin})},\xi _c^{(\mathrm{hyd})})<1`$ cannot be treated by our model.

Considering specific examples, we note that for the standard model, with $`L_{45}=1,R_{17}=1,l=0.3,U_x(R_0)=10^{0.8}`$, and $`\alpha =2`$, numerical calculations give $`\xi _c20`$. At $`\xi =4`$ the cloud is already moving at 90% of its final velocity. For two other cases of $`\alpha =2`$ ($`N_HR^2`$) and $`\alpha =1.5`$, the approximations hold throughout the region. The case of $`\alpha =10/3`$ is different due to the rapid expansion of the cloud. In this case, for $`N_H=\mathrm{const}`$, our model approximations break at small $`\xi `$ ($`\xi _c4`$).

Throughout this paper we assumed an external pressure that depends only on $`R`$. This is a natural assumption for the constant mass (thin spherical shell) case. For the other case, where the column density decrease with distance is slower than $`N_HR^2`$, we can imagine a situation where the cloud subtends a small solid angle ($`\mathrm{\Omega }<4\pi `$) and the external pressure profile is angle dependent. Clouds with finite lateral extents are prone to evaporation through their rims due to pressure gradients inside the cloud (e.g., Mathews & Veilleux 1989, equation 1). This sets a lower limit on the cloud’s lateral extents and conversly on the spatial angle it subtends of

$$\frac{\mathrm{\Omega }}{4\pi }>10^2\left(\frac{v_{\mathrm{sound}}}{v_{\mathrm{final}}}\right)^2.$$
(24)
Filamentary structure, and different dynamics in different directions, can be introduced into our formalism. The treatment of the general two dimensional problem is beyond the scope of this paper (for simple pancake shaped clouds see Blumenthal & Mathews 1979).

In this work we have considered general confinement, having in mind external magnetic pressure. Confinement by hot thermal gas is another possibility that has been discussed, extensively, in the literature. The presence of highly ionized, dilute gas between the clouds cannot be ruled out. Relative velocities between such gas and the cloud will result in drag forces that tend to decrease the final velocities. The cloud-gas boundary may be subjected to various instabilities which may affect the cloud structure considerably and ultimatly lead to its destruction. This complex situation cannot be solved using the methods adopted here.

## 4 Conclusions

We have made detailed numerical calculations of the dynamics of highly ionized gas clouds that are in pressure equilibrium with external magnetic pressure and are ionized by a typical AGN continuum. The self-consistent calculations include the ionization and thermal equilibrium of the cloud, the radiative transfer and the one-dimensional hydrostatical equilibrium solution. The principal conclusion is that such highly ionized clouds can be accelerated to high velocities by means of radiation pressure acceleration. The dynamical problem can be reduced to the calculation of a single parameter, the average force multiplier, $`M`$. The terminal velocity of the flow scales like the escape velocity at its origin but can exceed this velocity by a large factor ($`\sqrt{M}`$). The cloud velocity is sensitive to the confining pressure profile ($`\alpha `$), especially for marginal flows. In particular, flows that originate just outside the BLR in objects with $`L_{\mathrm{tot}}10^{45}\mathrm{erg}\mathrm{s}^1`$, will reach velocities of $`10003000\mathrm{km}\mathrm{s}^1`$, similar to the velocities measured for the X-ray and UV absorption lines. Highly ionized AGN clouds are driven mainly by bound-free absorption and bound-bound processes are less important unless significant non-thermal line broadening or very low column densities ($`<10^{20}\mathrm{cm}^2`$) are involved.

We thank Nahum Arav for useful discussions, and the referee for some valuable comments. We acknowledge financial support by the Israel Science Foundation and the Jack Adler Chair of Extragalactic Astronomy.


format: text

text:
# Heat capacity and pairing transition in nuclei

## I Introduction

Nuclear structure at low excitation energy depends critically on the presence of Cooper pairs. Thermal and rotational breaking of these $`J=0`$ nucleon pairs gives abrupt structural changes, such as increased level density and rotational-spin alignments. These critical phenomena were early addressed in several theoretical papers .

A very exciting feature is the gradual reduction of pair correlations as function of temperature. Recently, Schiller et al.  reported the experimental critical temperature $`T_c`$ for the pairing transition. The findings were based on using the canonical heat capacity as thermometer. An S-shaped heat capacity as function of temperature was observed in the 161,162Dy and 171,172Yb isotopes. Around $`T_c0.5`$ MeV, a local maximum in the heat capacity signals the breaking of Cooper pairs and quenching of the pair correlations. This property has also been observed in the 166,167Er nuclei .

Similar fingerprints have been obtained in various calculations. Finite-temperature Hartree Fock Bogoliubov calculations  for 164Er give almost identical S-shape as observed for 166Er. In relativistic Hartree Fock-BCS calculations  the proton and neutron pairing gaps are seen to vanish around $`T0.40.5`$ MeV for 166,170Er. Furthermore, in Shell Model Monte Carlo simulations (SMMC) , the heat capacities for iron isotopes show a pairing transition around temperatures of $`0.7`$ MeV.

The thermal breaking of a Cooper pair results in a tenfold increase in number of available energy levels. In this process particles are thermally scattered on available single particle states, giving rise to increased entropy. Recently, it was shown that each thermal particle carries an entropy of $`1.7`$, a feature which is valid for midshell nuclei with mass number $`A>40`$.

The present work aims to present a simple model for hot nuclei that includes the main features found experimentally. In Sec. II the model is described within the canonical-ensemble theory, and in Sec. III some relevant model properties are discussed. The results are compared with recent theoretical and experimental data on various deformed nuclei in Sec. IV. Concluding remarks are given in Sec. V.

## II Model

The total partition function $`Z`$ is described in the canonical ensemble, where thermal particle excitations, rotations, and vibrations are treated adiabatically according to

$$Z(T)Z_{\mathrm{par}}Z_{\mathrm{rot}}Z_{\mathrm{vib}}.$$
(1)
Thermodynamical quantities such as the entropy $`S`$, the average excitation energy $`E`$ and the heat capacity $`C_V`$ can then be calculated from the Helmholtz free energy:

$$F(T)=T\mathrm{ln}Z(T),$$
(2)
by

$`S(T)`$ $`=`$ $`\left({\displaystyle \frac{F}{T}}\right)_V,`$ (3)
$`E(T)`$ $`=`$ $`F+TS,`$ (4)
$`C_V(T)`$ $`=`$ $`\left({\displaystyle \frac{E}{T}}\right)_V,`$ (5)
where the Boltzmann constant is set to unity ($`k_B=1`$) and $`T`$ is measured in units of MeV.

The theoretical basis for the particle partition function $`Z_{\mathrm{par}}`$ was earlier presented in Ref. , however, for completeness we make a short outline here. The particle excitations are described by spin $`1/2`$ fermions scattered into doubly degenerated single-particle levels with equal energy spacing $`ϵ`$. There are infinitely many single-particle levels with one set of levels for protons and one set for neutrons. The partition function for one spin-up fermion is given by

$$z_1^{}=\underset{i=0}{\overset{\mathrm{}}{}}\mathrm{exp}(iϵ/T)=\frac{1}{1\mathrm{exp}(ϵ/T)}.$$
(6)
When two identical spin-up fermions are placed into the same single-particle level scheme, the summing is restricted due to the Pauli principle by

$`z_2^{}`$ $`=`$ $`{\displaystyle \underset{i=0}{\overset{\mathrm{}}{}}}\left[\mathrm{exp}(iϵ/T){\displaystyle \underset{j=i+1}{\overset{\mathrm{}}{}}}\mathrm{exp}(jϵ/T)\right]`$ (7)
$`=`$ $`{\displaystyle \frac{1}{1\mathrm{exp}(2ϵ/T)}}{\displaystyle \frac{\mathrm{exp}(ϵ/T)}{1\mathrm{exp}(ϵ/T)}}`$ (8)
$`=`$ $`z_1^{}{\displaystyle \frac{\mathrm{exp}(ϵ/T)}{1\mathrm{exp}(2ϵ/T)}}.`$ (9)
Generally, for $`n`$ identical spin-up fermions the partition function is

$$z_n^{}=z_{n1}^{}\frac{\mathrm{exp}((n1)ϵ/T)}{1\mathrm{exp}(nϵ/T)}.$$
(10)
We now introduce spin up and down for the fermions, and evaluate the corresponding partition function $`z_n`$. For one fermion with spin up or down we obtain a degeneration of two, since $`z_1^{}=z_1^{}`$:

$$z_1=z_1^{}+z_1^{}=2z_1^{}.$$
(11)
If two fermions occupy the same single-particle level scheme, the partition function becomes more complicated. With the spins of the two fermions antiparallel ($`m=m_1+m_2=0`$), no Pauli blocking is present and

$$z_2(m=0)=z_1^{}z_1^{}=\left(z_1^{}\right)^2.$$
(12)
For the contribution with parallel spin, the Pauli principle gives

$$z_2(m=\pm 1)=z_2^{}+z_2^{}=2z_2^{}.$$
(13)
It is now straightforward to evaluate $`z_n`$ for $`n`$ fermions, allowing spin up and down. The five lowest partition functions are:

$`z_1`$ $`=`$ $`2z_1^{},`$ (14)
$`z_2`$ $`=`$ $`2z_2^{}+\left(z_1^{}\right)^2,`$ (15)
$`z_3`$ $`=`$ $`2z_3^{}+2z_2^{}z_1^{},`$ (16)
$`z_4`$ $`=`$ $`2z_4^{}+2z_3^{}z_1^{}+\left(z_2^{}\right)^2,`$ (17)
$`z_5`$ $`=`$ $`2z_5^{}+2z_4^{}z_1^{}+2z_3^{}z_2^{}.`$ (18)
The partition function for nuclei with protons and neutrons may now be constructed from $`z_n`$. In, e.g., the three-particle case of an odd-proton system the $`\pi \nu ^2`$ and the $`\pi ^3`$ partition functions have to be included. For an even-even two-particle system a $`\pi `$ or a $`\nu `$ pair can be broken, and so on. The lowest partition functions for odd-mass ($`Z_{1,3,5}`$) and even-even systems ($`Z_{2,4}`$) become:

$`Z_1`$ $`=`$ $`z_1,`$ (19)
$`Z_2`$ $`=`$ $`2z_2,`$ (20)
$`Z_3`$ $`=`$ $`z_3+z_1z_2,`$ (21)
$`Z_4`$ $`=`$ $`2z_4+(z_2)^2,`$ (22)
$`Z_5`$ $`=`$ $`z_5+z_3z_2+z_1z_4.`$ (23)
The partition functions for the odd-odd system can be constructed in a similar manner. The two-particle partition function is $`(z_1)^2`$, giving an entropy of twice the single particle entropy of the odd-mass system. Higher numbers of particles yield the following partition functions for the odd-odd case:

$`\stackrel{~}{Z}_2`$ $`=`$ $`(z_1)^2,`$ (24)
$`\stackrel{~}{Z}_4`$ $`=`$ $`2z_3z_1,`$ (25)
$`\stackrel{~}{Z}_6`$ $`=`$ $`2z_5z_1+(z_3)^2.`$ (26)
The creation of particles costs energy. In order to break one Cooper pair, the energy $`2\mathrm{\Delta }`$ is necessary. The pairing gap parameter is given by the empirical formula

$$\mathrm{\Delta }=12A^{1/2}\mathrm{MeV}.$$
(27)
The total particle partition functions of even-even (ee), odd (odd) and odd-odd (oo) nuclei can then be expressed as

$`Z_{\mathrm{par}}^{\mathrm{ee}}`$ $`=`$ $`1+Z_2e^{2\mathrm{\Delta }/T}+Z_4e^{4\mathrm{\Delta }/T}+\mathrm{},`$ (28)
$`Z_{\mathrm{par}}^{\mathrm{odd}}`$ $`=`$ $`Z_1+Z_3e^{2\mathrm{\Delta }/T}+Z_5e^{4\mathrm{\Delta }/T}+\mathrm{},`$ (29)
$`Z_{\mathrm{par}}^{\mathrm{oo}}`$ $`=`$ $`\stackrel{~}{Z}_2+\stackrel{~}{Z}_4e^{2\mathrm{\Delta }/T}+\stackrel{~}{Z}_6e^{4\mathrm{\Delta }/T}+\mathrm{}`$ (30)
The number of active particles $`n`$ (i.e. particles not bound in Cooper pairs) can be evaluated by

$`n^{\mathrm{ee}}`$ $`=`$ $`(2Z_2e^{2\mathrm{\Delta }/T}+4Z_4e^{4\mathrm{\Delta }/T}+\mathrm{})/Z_{\mathrm{par}}^{\mathrm{ee}},`$ (31)
$`n^{\mathrm{odd}}`$ $`=`$ $`(Z_1+3Z_3e^{2\mathrm{\Delta }/T}+5Z_5e^{4\mathrm{\Delta }/T}+\mathrm{})/Z_{\mathrm{par}}^{\mathrm{odd}},`$ (32)
$`n^{\mathrm{oo}}`$ $`=`$ $`(2\stackrel{~}{Z}_2+4\stackrel{~}{Z}_4e^{2\mathrm{\Delta }/T}+6\stackrel{~}{Z}_6e^{4\mathrm{\Delta }/T}+\mathrm{})/Z_{\mathrm{par}}^{\mathrm{oo}}.`$ (33)
The rotational partition function is given by

$$Z_{\mathrm{rot}}=\underset{I=0,2,4\mathrm{}}{}\mathrm{exp}\left[A_{\mathrm{rig}}I(I+1)/T\right],$$
(34)
where the rotational parameter

$$A_{\mathrm{rig}}=\mathrm{}^2/2\theta _{\mathrm{rig}}$$
(35)
is expressed by the rigid moment of inertia  $`\theta _{\mathrm{rig}}=2/5MR^20.0137A^{5/3}`$ $`\mathrm{}^2`$MeV-1, where $`M`$ and $`R`$ are the nuclear mass and radius, respectively. The vibrational partition function includes zero and one-phonon states:

$$Z_{\mathrm{vib}}=\underset{v=0,1}{}W_v\mathrm{exp}(v\mathrm{}\omega _{\mathrm{vib}}/T),$$
(36)
where the multiplicity for the zero-phonon state is $`W_0=1`$. For one-phonon states, we take into accounts three vibrations (e.g. $`\beta `$, $`\gamma `$ and octupole vibrations), giving multiplicity $`W_1=3`$. Higher order phonon states are neglected, and the phonons are assumed to carry the same energy quantum $`\mathrm{}\omega _{\mathrm{vib}}`$. In both partition functions we have omitted the spin degeneracy ($`2I+1`$), since we are dealing with levels and not states.

## III Model properties

Figure 1 shows the Helmholtz free energy $`F`$, the average excitation energy $`E`$, the entropy $`S`$ and the heat capacity $`C_V`$ as functions of temperature. The model parameters are taken from 162Dy (see Table I) as a typical set applicable for rare earth nuclei. In order to calculate the thermodynamical quantities up to $`T1`$ MeV, at least 10 broken nucleon pairs have to be incorporated. The free energy $`F`$ and the average excitation energy $`E`$ behave smoothly as functions of temperature. Around $`T0.65`$ MeV the nuclei are excited to energies comparable to their respective neutron binding energies. The even-even, oddWe use the abbreviation odd for odd-even and even-odd systems, since these systems are equivalent in our model. and odd-odd systems have different excitation energies at one and the same temperature, where the even-even system requires the highest $`E`$ value.

The entropy $`S`$ and heat capacity $`C_V`$ represent first and second derivatives of $`F`$, and are thus more sensible to thermal changes. For the lowest temperatures the entropy difference is $`2`$ between the three mass systems. However, the entropy curves coincide for $`T>0.6`$ MeV, displaying almost identical behavior. It is interesting to test if our model reproduces a Fermi gas with the absence of pairing at these temperatures§§§A simplified Fermi gas has $`S2aT+\mathrm{const}.`$. The level density parameter $`a`$ for a Fermi gas with a level spacing $`ϵ`$ is given by

$$a=\frac{\pi ^2}{6}(g_p+g_n)\frac{\pi ^2}{3}g=\frac{\pi ^2}{3ϵ}.$$
(37)
Here, the single-particle level-density parameters for protons and neutrons ($`g_p`$ and $`g_n`$) are assumed to be approximately equal. Inserting $`ϵ=0.13`$ MeV for 162Dy, we obtain $`a=25`$ MeV-1. This is in exact agreement with the slope ($`a=S/2T`$) of the entropy curves for $`T>0.6`$ MeV.

The lower right panel of Fig. 1 shows the typical S-shape of the heat capacity. As this shape is an important fingerprint for pairing transitions in nuclei , we will in the following focus on its origin.

The contribution to $`C_V`$ from collective excitations is negligible, and has no influence on the S shape. This is shown in Fig. 2 for the even-even system, where the components of particle, rotational, and vibrational degrees of freedom are displayed. The collective contribution to the entropy $`S`$ is small and fairly constant with increasing temperature. However, one should note that $`S_{\mathrm{rot}}`$ is in fact the main component at the lowest temperatures with $`T<0.3`$ MeV.

For the total entropy and heat capacity, also the effect of the number of pairs is shown in Fig. 2. The $`C_V`$ curves are identical for 6 and 10 pairs up to $`T0.7`$ MeV, but from there on the 6 pair system is exhausted and not able to absorb energy at the same rate as when more Cooper pairs are present. The figure shows that the S shape can be evaluated rather accurately without taking very many pairs into account. On the other hand, if too few pairs are included, one may easily misinterpret the shape of the $`C_V`$ curve. For the 162Dy mass region, 6 Cooper pairs are sufficient to determine the S shape.

Figure 3 shows the heat capacity (upper panels) and the corresponding number of unpaired particles (lower panels) for 10 pairs of nucleons. In the left panels a realistic pairing gap parameter of $`\mathrm{\Delta }=0.94`$ MeV has been chosen. We identify three temperature regions of interest: (i) The $`T0.20.8`$ MeV region with different S shapes for even-even, odd and odd-odd-systems, (ii) the Fermi gas regime with $`T0.81.0`$ MeV, exhibiting the same linear heat capacity for the three systems and (iii) the dramatical change for $`T>1`$ MeV due to the final number of particles in the systems.

The pronounced $`C_V`$ maximum at $`T1.01.2`$ MeV is due the low number of available Cooper pairs; the lower left panel of Fig. 3 shows that $`15`$ nucleons are already unpaired at this temperature. As $`T\mathrm{}`$, we find $`C_V`$ 20, 21 and 22 for the even-even, odd and odd-odd systems, respectively. This corresponds to the situation where all pairs are broken, giving a $`C_V`$-value equal to the total number of nucleons in the systemFor the respective systems, we have 10 Cooper pairs plus 0, 1 or 2 particles.. Since a system with one particle in an infinite harmonic oscillator gives $`C_V1`$ for $`T\mathrm{}`$ and the effective Pauli blocking between the nucleons is negligible, we obtain $`C_V=n`$.

In the right panels of Fig. 3, the pairing gap is reduced to $`\mathrm{\Delta }=ϵ=0.13`$ MeV. Now the S shape of the $`C_V`$ curve can be seen to vanish in the $`\mathrm{\Delta }0`$ limit, telling that the S shape is connected with the pairing strength. Furthermore, we see that the heat capacity reaches a maximum level of $`C_V35`$ at lower temperature $`T0.7`$ MeV. At this point the remaining number of particles bound in Cooper pairs is low, and the depairing mechanism again looses its capability to create more heat capacity. It might be surprising though that the number of depaired particles does not increase much faster with temperature. The answer to this is that in our model, not only the energy $`2\mathrm{\Delta }`$ is required to break up a pair, but also the unpaired particles have to be placed into the single-particle level scheme at some finite energy $`E`$, since the Pauli blocking prevents them to populate the lowest levels in the single-particle level scheme, which are all occupied by other unpaired nucleons.

In Fig. 4 an interpretation of the S shape is given for the even-even system. In the left panel, a ”background” is subtracted from $`C_V`$, using a straight line through 0.3 and 0.9 MeV of temperature. This line is intended to mimic the underlying heat capacity originating from a Fermi gas type of system. However, the gas properties are strongly connected to the depairing process, and the subtracted peak in the lower left panel should be taken with care.

The resulting peak shows some resemblance with the Schottky anomaly, which describes a particle placed in a two level system . However, the maximum heat capacity of such a model is $`C_V0.45`$ at a temperature of 40% of the energy gap between the levels. For atomic nuclei, we may define a similar, but extended Schottky model: Either no pairs are broken at excitation energy $`E=0`$, or one pair is broken at $`E=2\mathrm{\Delta }`$, or two pairs are broken at $`E=4\mathrm{\Delta }`$ and so on. This picture is an extreme simplification, and we include here only zero, one and two pairs. Thus, our Schottky-like partition function reads

$$Z_{\mathrm{Sch}}=1+W_{2\mathrm{\Delta }}e^{2\mathrm{\Delta }/T}+W_{4\mathrm{\Delta }}e^{4\mathrm{\Delta }/T}.$$
(38)
The multiplicities in front of the Boltzmann factors represent the number of levels at $`E=2\mathrm{\Delta }`$ and $`E=4\mathrm{\Delta }`$ and are estimated from experimental data using $`W=\rho \delta E`$, where $`\rho `$ is the level density and $`\delta E`$ is the energy window considered.

In the right panels of Fig. 4 $`C_V`$ is evaluated by the use of Eq. (38) for 162Dy. The $`\rho `$ values are estimated from Table I of Ref. . With an energy window of $`\delta E=0.2`$ MeV, we obtain $`W_{2\mathrm{\Delta }}=22`$ and $`W_{4\mathrm{\Delta }}=556`$; numbers which are also consistent with the numbers of seniority $`𝒮=`$ 2 and 4 states in the model of Ref. . The high multiplicities give a strength, width, and position of the peak that compare qualitatively well with the peak in the lower left panel. From this comparison, we interpret the local maximum of the S-shaped $`C_V`$ curve at $`T0.5`$ MeV as the point at which the depairing process is at the strongest. Hence, we define this point as the critical temperature $`T_c`$ for the pairing transition.

To close this section, we comment on the weakening of the S shape when going from even-even to odd and odd-odd systems; an effect which is clearly seen in the lower right panel of Fig. 1. The reason for the weakening can be understood from the entropy plot displayed in the right upper panel of Fig. 1. There, we saw that the entropy of each valence nucleonInterpreted as the entropy gap between the various S-curves . is reduced from a value of $`S2`$ at $`T0.3`$ MeV to zero at $`T0.6`$ MeV. Since the corresponding heat capacity relates to $`TS/T`$, a negative contribution to $`C_V`$ appears and a quenching of the amplitude of the S-curve is apparent.

## IV Comparison with data

Our model is described within the canonical ensemble, while experimental data refer to the microcanonical ensemble. However, there are two ways to compare our model with experiments. With a known experimental level density $`\rho `$, the partition function can be constructed from

$$Z(T)=\underset{i}{}\delta E_i\rho (E_i)e^{E_i/T},$$
(39)
where $`E_i`$ is the excitation energy and $`\delta E_i`$ are the energy bins used. From this partition function all thermodynamical quantities can be deduced, see e.g. Eqs. (2-5). The drawback of this method, is that the level density function has to be known up to high excitation energies, typically $`E40`$ MeV. The other way is to evaluate the microcanonical level density from our canonical model. This can be performed by an inverse Laplace transformation of $`Z`$. Using the saddle-point approximation (Fawler-Darwin method), we obtain

$$\rho (E)=\frac{e^S}{T\sqrt{2\pi C}}.$$
(40)
The level density can then be compared with certain experimental values, and complete knowledge on the level density up to high excitation energies is not necessary. One drawback using the saddle-point approximation is that the excitation energy $`E`$ is a thermal average with a large standard variation of $`\sigma _E=T\sqrt{C}T\sqrt{2aT}`$, with $`a`$ being the level-density parameter. A second drawback is the approximation itself, which we have tested by a ”forward-backward” Laplace transformation. In Fig. 5 we compare the original entropy and heat capacity, with the ones obtained by using Eq. (40) to obtain $`\rho `$, and then using Eq. (39) to obtain a new $`Z`$ and its corresponding $`S`$ and $`C_V`$. The comparison reveals a general smoothing for temperatures $`T<0.50.6`$ MeV. The entropy is seen to be reproduced rather well, however, the heat capacity is more sensitive to the approximation.

The model presented in Sec. II rests on the assumption that the single-particle level scheme can be approximated by equidistant levels. This assumption is never fulfilled in atomic nuclei, but within the canonical ensemble, the various deduced thermodynamical quantities are strongly smoothed with respect to excitation energy. Therefore, even with some non-uniformity, the model might still give realistic results. Probably, heavy and strongly deformed nuclei are the best candidates for our model.

For test cases, we have chosen midshell nuclei around 58Fe, 106Pd, 162Dy and 234U. All four mass regions reveal a rather uniform Nilsson single-particle energy distribution without large energy gaps. The pairing gap $`\mathrm{\Delta }`$ and rotational parameter $`A_{\mathrm{rot}}`$ are calculated according to Eqs. (27) and (35). The vibrational energy quantum $`\mathrm{}\omega _{\mathrm{vib}}`$ is taken as the energy at which the first vibrational states appear in the experimental level schemes of the respective mass region . The last parameter needed is the level-gap parameter $`ϵ`$, which is expected to have a value in between the single-particle spacing $`ϵ_{\mathrm{sp}}`$ and the BCS quasiparticle spacing $`ϵ_{\mathrm{qp}}=\sqrt{(ϵ_{\mathrm{sp}}\lambda )^2+\mathrm{\Delta }^2}\mathrm{\Delta }`$, where $`\lambda `$ is the Fermi level.

In this work, $`ϵ`$ is chosen as a free parameter, determined from a fit to known experimental level densities. For each nucleus we adopt two level density anchor points, as deduced in Ref. . The lower anchor point is based on the counting of known discrete levels. This method is rather accurate, except for the odd-odd nuclei, where the number of levels might be several hundred per MeV, and thus difficult to measure. The other anchor point is based on average neutron-resonance spacing data at the neutron-binding energy.

Figure 6 shows anchor points and level densities\**\**\**The experimental data are taken at a given excitation energies $`E`$, while the calculations give $`\rho `$ as function of average excitation energy $`E`$. calculated using Eq. (40). The parameters are listed in Table I, where $`ϵ`$ is adjusted to obtain the approximate slope of $`\rho `$ in the log-plot. The agreement with the anchor points is good, and also the odd and odd-odd systems fall nicely into the systematics. For the 161,162Dy isotopes, the experimental level densities from Ref.  are shown as well. The calculations reveal good agreement with experiment as function of excitation energy.

With the established parameters, we may now calculate the heat capacity $`C_V`$ from our model and compare with experiments. Due to the imperfection of the saddle-point approximation, we take $`\rho (E)`$ from the saddle-point approximation (which is tuned to experiments in Fig. 6), and generate $`C_V`$ from the corresponding $`Z`$ function of Eq. (39). In the lower part of Fig. 7 the theoretical $`C_V`$ curves are shown for 161,162Dy. The experimental $`C_V`$ curves, shown in the upper part, are based on experimental level density data from Ref. . Since the construction of $`Z`$ requires data to much higher energies than experimentally known ($`E78`$ MeV), the level density has to be extrapolated. Here, we have used the parametrization of Egidy et al. , in accordance with our recent work . The effective level density parameter ($`a17.5`$ MeV-1) has a major impact on the $`C_V`$ curve for $`T>0.5`$ MeV. Figure 7 shows rather good agreement between experiment and model for $`T<0.50.6`$ MeV. In particular the odd-even mass entropy difference is well reproduced. At higher temperatures the comparison is poor due to the arbitrary extrapolation of the experimental level densities.

In Fig. 8 the heat capacity $`C_V`$ from our model, using the parameters of Table I, is displayed for the four mass regions. The curves look similar; the differences are mainly the change in the scaling of the $`C_V`$ and $`T`$ axes for the various mass regions. The figure also includes the calculations on 58,59Fe performed by Liu and Alhassad . These data points show some discrepancies with the work of Rombouts et al. . However, both SMMC calculations obtain a critical temperature around 0.7 MeV, while we obtain a value around 1.2 MeV. Even when changing freely the number of particles, $`ϵ`$ and $`\mathrm{\Delta }`$ in our model, we are not able to reproduce the SMMC results. In particular, the feature that $`C_V`$ becomes constant , or even drops  above $`T0.7`$ MeV is surprising. This effect might originate from shell gaps (which are not included in our model) or too few particles and/or orbitals considered in the SMMC calculations. It would be interesting to test if the SMMC calculations are capable of reproducing the anchor points for Fe in Fig. 6. Doing so, one should remember that our anchor points represent densities of levels and not states, which are usually employed in SMMC calculations.

## V Conclusions

A simple model for hot nuclei has been outlined. The main properties of the model are determined by the pairing-gap parameter $`\mathrm{\Delta }`$ and the energy-gap parameter $`ϵ`$ associated with an infinite single-particle level scheme for protons and neutrons. We have demonstrated how various thermodynamical quantities can be extracted.

The model properties have been discussed with emphasize on the S shape of the heat capacity. This shape is strongly related to the bunch of newly created levels from the pair-breaking process at $`2\mathrm{\Delta }`$ and $`4\mathrm{\Delta }`$ of excitation energy.

The model calculations are compared with experimental data from the $`A`$ 58, 106, 162 and 234 mass regions. Using the saddle-point approximation with only one free parameter, the experimental level densities of even-even, odd and odd-odd systems are reproduced.

The critical temperature of the pair transition for 161,162Dy was calculated to be $`T_c0.5`$ MeV, in agreement with experiments. In the 58,59Fe region we calculate $`T_c1.2`$ MeV, which is significantly higher than the value of $`0.7`$ MeV obtained in Shell Model Monte Carlo simulations. Thus, further theoretical and experimental efforts are needed to understand the thermodynamics of these hot iron isotopes.

We wish to acknowledge the support from the Norwegian Research Council (NFR).


text:
# Anomalous self-similarity in a turbulent rapidly rotating fluid

## Abstract

Our velocity measurements on the quasi-two-dimensional turbulent flow in a rapidly rotating annulus yield an inverse energy cascade with $`E(k)k^2`$ rather than the expected $`E(k)k^{5/3}`$. The probability distribution functions for longitudinal velocity differences, $`\delta v(\mathrm{})=v(x+\mathrm{})v(x)`$, are self-similar (scale independent) but strongly non-Gaussian, which suggests that the coherent vortices play a significant role. The structure functions, $`[\delta v(\mathrm{})]^p\mathrm{}^{\zeta _p}`$, exhibit anomalous scaling: $`\zeta _p=\frac{p}{2}`$ rather than the expected $`\zeta _p=\frac{p}{3}`$.

In large scale flows in the earth’s atmosphere and oceans or in gaseous planets, the Coriolis force dominates other forces and to lowest order is balanced by pressure gradients (geostrophic balance). The dimensionless number that characterizes this regime is the Rossby number, the ratio of magnitudes of the inertial term in the Navier-Stokes equation ($`\stackrel{}{u}\stackrel{}{u}`$, where $`\stackrel{}{u}`$ is the velocity) to the Coriolis term ($`2\stackrel{}{\mathrm{\Omega }}\times \stackrel{}{u}`$, where $`\mathrm{\Omega }`$ is the rotation rate). For planetary flows on large scales the Rossby number is typically in the range 0.05-2, and turbulence in such flows is quite different from three-dimensional (3D) turbulence in an inertial reference frame. Despite the importance of planetary flows, geostrophic turbulence has scarcely been examined in laboratory experiments. Our experiments in a rotating annulus are the first to determine the scaling properties of turbulence in a low Rossby number flow.

One of the most significant effects of rapid rotation on a fluid is the two-dimensionalization of the flow . Recent numerical simulations of a rotating turbulent flow recovered quasi-2-dimensionality even in the case of 3D forcing . Turbulence in a 2D rapidly rotating flow has some similarities to other 2D turbulent flows, but the rapid rotation leads to the formation of large vortices and jets, hence a flow that is more inhomogeneous and anisotropic .

The statistics of velocity differences for 2D turbulence are predicted to be self-similar , i.e., the probability distribution function (PDF) for the difference between velocities measured at two points, $`\delta v(\mathrm{})=v(x+\mathrm{})v(x)`$ (where $`v`$ is along the line connecting the two points), should have a functional form independent of the separation $`\mathrm{}`$. Self-similarity in a turbulent flow can also be determined from the scaling of structure functions, $`S_p(\mathrm{})<[\delta v(\mathrm{})]^p>\mathrm{}^{\zeta _p}`$: for any self-similar flow, $`\zeta _p`$ will vary linearly with $`p`$ . In particular, Kolmogorov’s 1941 theory for homogeneous isotropic turbulence predicts the existence of an inertial range with $`\zeta _p=\frac{p}{3}`$.

Experiments and numerical simulations have shown that 3D turbulence is $`\mathrm{n}\mathrm{o}\mathrm{t}`$ self-similar . The observed deviations from a linear dependence of $`\zeta _p`$ on $`p`$ are attributed to the stretching and folding of vortices. As a vortex stretches, it shrinks in lateral extent until it collapses onto a singular line. Consequently, the energy dissipation rate varies with length scale. Vortex stretching is not allowed in 2D turbulence, and recent numerical simulations of 2D turbulence yielded self-similar behavior  (see also ). A recent experiment yielded self-similar behavior with a Gaussian PDF  while two other experiments showed scale-dependent statistics , as will be discussed later.

Our experiments on turbulent flow in a rapidly rotating annulus yield a self-similar PDF for the velocity differences, but the PDF is not Gaussian and the structure function exponent values $`\zeta _p`$ differ from the expected $`p/3`$. We now describe our experiments and then present the results.

Our apparatus consists of an annular tank filled with water and covered with a solid lid; the inner diameter of the tank is 21.6 cm and the outer diameter is 86.4 cm . The depth of the tank increases from 17.1 cm at the inner radius to 20.3 cm at the outer radius . Flow in the annulus is produced by continuously pumping water in a closed circuit through two concentric rings of 120 holes each in the bottom of the tank; the source ring is at a radius of 18.9 cm and the sink ring is at 35.1 cm. The radially outward flux from the pumping couples with the Coriolis force to generate a counter-rotating jet, which is wide and turbulent over a wide range of parameters. The rapid rotation leads to a flow that is essentially 2D except in the thin Ekman boundary layers at the top and bottom surfaces .

In the present experiments the tank rotates at 11.0 rad/s, sufficiently fast to produce an essentially 2D flow, and the flux is 150 cm3/s, sufficiently large to produce an inverse energy cascade. The azimuthal velocity is measured using two hot film probes which are inserted through the top lid and extend 1 cm into the water on opposite sides of the tank, midway between the inner and outer walls. Each probe was sampled at 150 Hz for periods of two hours, giving $`10^6`$ data points per probe for each run, and the measurements were repeated four times, yielding a total of $`8\times 10^6`$ data points. Using the maximum velocity ($`U_{max}=22`$ cm/s) as the velocity scale and the distance between the forcing rings ($`L=16.2`$ cm) as the integral length scale, we calculate the Reynolds number to be $`35,000`$ and the Rossby number ($`U_{max}/2\mathrm{\Omega }L`$) to be $`0.06`$.

Instantaneous 2D velocity fields were obtained using a Particle Image Velocimetry (PIV) system with a horizontal light sheet at mid-fluid depth and a rotating camera above the tank. For each flow condition, 50 instantaneous velocity fields were obtained, equivalent to approximately $`2\times 10^5`$ velocity values at the radius of the hot film probes. Though this sample size is inadequate for high order statistics, the spatial information provided by the PIV measurements complements the long velocity time series obtained with the hot film probes. We also made PIV measurements at the same rotation rate at a higher pumping rate of 550 cm3/s ($`Re=100,000`$).

Two vorticity profiles obtained with the PIV system are shown in Fig. 1. When the pump is first turned on, the flow consists of rings of cyclonic and anti-cyclonic vortices that form above the outlet and inlet rings, respectively. Vortices of like sign immediately begin merging and growing in size, and a counter-rotating jet forms between the two rings. At long times the inverse energy cascade leads to large vortices which are larger at higher pumping rates, as Fig. 1 illustrates.

We compute energy power spectra from the velocity time series data assuming Taylor’s frozen turbulence hypothesis, which is applicable because the turbulent intensity (ratio of the rms velocity fluctuation to the mean velocity) is less than 10%. The spectra contain a region with $`E(k)k^2`$ (Fig. 2), in contrast with Kraichnan’s prediction of $`E(k)k^{5/3}`$ for the inverse cascade . Energy spectra obtained from PIV measurements of the azimuthal velocity data at $`Re=100,000`$ show a scaling consistent with that obtained from the time series data, as shown in the upper curve of Fig. 2. A prediction of $`E(k)k^2`$ has been made for a turbulent flow made up of Lundgren spiral vortices that were not allowed to stretch in the third dimension . Such scaling has also been predicted for rotating turbulent flows as they become 2D .

At high wavenumbers, our power spectra are consistent with those previously found for the forward (enstrophy) cascade, $`E(k)k^n`$, where $`3n4`$. However, our spectral range is too small to deduce a value for $`n`$.

By plotting the PDF for the velocity differences $`\delta v(\mathrm{})`$ (Fig. 3), we obtain the first indication of self-similarity: data for different separations $`\mathrm{}`$ collapse onto a single curve when normalized. This curve is far from Gaussian, and the enhanced probability in the tails is likely due to the strong velocity differences that arise as coherent vortices pass the probes, as noted by She et al. .

Another indication of anomalous behavior is the scaling of the standard deviation of the velocity differences, $`\delta v_{rms}`$. Our data suggest $`\delta v_{rms}(\mathrm{})\mathrm{}^{1/2}`$ (see insets of Fig. 3) rather than the expected $`\delta v_{rms}(\mathrm{})\mathrm{}^{1/3}`$; the observed scaling of $`\delta v_{rms}`$ is consistent with the scaling of the energy spectrum in the inverse cascade range (Fig. 2).

Now we examine the structure function scaling, plotting $`S_p`$ as a function of $`\mathrm{}`$, as shown in Fig. 4(a). There is a scaling region (labeled A) from about 2 cm to 8 cm, and some indication of another scaling region (labeled B) for $`\mathrm{}<2`$ cm; the transition between the two regions can be seen more clearly in the plot of $`S_{10}/\mathrm{}^{5.5}`$ vs. $`\mathrm{}`$ in the inset. The 2 cm lower limit of scaling of the structure functions is below the 5 cm ($`k=1.2`$ cm-1 in Fig. 2) length where the energy spectrum behavior changes; this difference might be due to systematic error in calibration of the probes (note in Fig. 2 the difference in two probes at high wavenumbers), or perhaps because the scaling regions in real space and in Fourier space do not exactly correspond, as discussed in Ref.  (p. 62).

One test for the existence of an inverse energy cascade in 2D turbulence is the sign of $`S_3(\mathrm{})`$: Kolmogorov’s four-fifths law can be written for anisotropic flows as $`\epsilon =\frac{1}{4}_{\mathrm{}}|\delta 𝐯(\mathrm{})|^2\delta 𝐯(\mathrm{})`$, where $`𝐯`$ is the full 3D velocity vector and $`\epsilon `$ is the mean rate of energy transfer (Ref. , p. 88). For the 3D forward cascade $`\epsilon `$ is positive; a negative value of $`\epsilon `$ corresponds to an inverse cascade. If the anisotropy is not too strong, the longitudinal $`S_3`$ dominates the transverse $`S_3`$, and one can obtain the sign of $`\epsilon `$ from measurements of the longitudinal structure function alone. While our flows are anisotropic, we find that PDFs for velocity differences from PIV measurements in the radial and azimuthal direction are in agreement, indicating that the anisotropy does not strongly affect the velocity differences. The hot film data yield $`S_3>0`$ for $`\mathrm{}<10`$ cm, which suggests that the inverse energy cascade stops at that point. This length corresponds visually to the scale of the largest vortices; thus region A corresponds to the inverse cascade. Scales below about 1 cm (the distance between the holes in the inner ring) presumably correspond to the forward enstrophy cascade.

The now standard way of extracting structure function exponents from data with limited inertial range is the technique introduced by Benzi et al.  called Extended Self-Similarity (ESS), where the $`\zeta _p`$ values are obtained from the slope of log-log plots of $`S_p`$ vs. $`S_3`$. The ESS plots of our data, Fig. 4(b), exhibit power law scaling through both ranges A and B, for lengths $`0.5<\mathrm{}<15`$ cm. Since $`S_3`$ should contain the same trends as the other structure functions, it is not surprising that the switching from region A to B is lost by this method.

The exponent values $`\zeta _p`$ deduced from the plot of $`S_p(\mathrm{})`$ \[region A in Fig. 4(a)\] and from the ESS plot \[Fig. 4(b)\] are compared in Fig. 5. Region A of Fig. 4(a) yields $`\zeta _{p,A}=p/2`$, in contrast with theory, although $`\zeta _{p,ESS}=p/3`$, as it must for a self-similar flow.

Two recent experiments on turbulence in a quasi-2D fluid, a soap film, found departures from linear scaling of $`\zeta _p`$ vs. $`p`$ , indicating that these flows were not self-similar. Another experiment on a quasi-2D flow, a magnetically driven shallow layer of electrolyte, yielded an inverse energy cascade with self-similar behavior as indicated by both the PDFs for $`\delta v`$ and by the scaling of the structure function computed using ESS, which scaled as $`\zeta _p=p/3`$ ; however, the PDFs for $`\delta v`$ were close to Gaussian, in contrast to those yielded by our flow (Fig. 3). A recent numerical simulation of 2D turbulence also yielded self-similar PDFs for $`\delta v`$, where the departure from Gaussianity, though small, contained significant information about the flow .

The $`p/2`$ scaling obtained for $`\zeta _p`$ in our flow may be a consequence of an inverse cascade driven by a radial velocity shear. This process would yield an energy flux $`\epsilon \delta v_s(\delta v_{})^2/\mathrm{}`$, where $`v_s`$ would be a radial velocity due to the interaction of the vortices with the mean shear, $`v_{}`$ the azimuthal velocity and $`\mathrm{}`$ the azimuthal separation. For a scale-invariant $`\epsilon `$, the resulting $`p^{th}`$ order velocity structure function would have the scaling exponent $`p/2`$.

In summary, we find that turbulence in our rapidly rotating flow is self-similar, as indicated by the collapse of the PDFs onto a single curve and by the linear scaling of the structure function exponents. This self similarity extends over both the inverse and forward cascades. In the inverse cascade range, the flow exhibits anomalous scaling of the energy, $`E(k)k^2`$, and the structure functions, $`S_p\mathrm{}^{p/2}`$, but this anomalous scaling is missed in an Extended Self Similarity analysis. We conclude that the rotation and strong shear make this flow different from 2D turbulent flows in non-rotating systems. Our results suggest that the transport and mixing in low Rossby number flows such as the atmosphere and oceans may be different from these processes in isotropic homogeneous 2D turbulent flows.

This research was supported by a grant from the Office of Naval Research. Z-S. S. acknowledges partial support from the Minister of Education in China and the Natural Science Foundation of China.