text:
# Diffusion in supersonic, turbulent, compressible flows

## I Introduction

Laboratory and terrestrial gases and liquids are usually well described by incompressible flows (e.g. Lesieur ). In contrast, the dynamical behavior of typical astrophysical gases, are characterized by poorly understood highly compressible supersonic turbulent motion (see e.g. VS00 ). For example, the large observed linewidths in large molecular clouds show direct evidence for the presence of chaotically oriented velocity fields with magnitudes in excess of the sound speed. This random motion carries enough kinetic energy to counterbalance and sometimes overcompensate the effects of self-gravity of these clouds WBM00 . The intricate interplay between supersonic turbulence and self-gravity determines the overall dynamical evolution of these clouds and their observable features such as their density structure, the star formation rate within them, and their lifetimes. Thus, it is importance for the description of many astrophysical systems to understand in detail the momentum and heat transfer properties of compressible turbulent gases.

Some important clues on the nature and efficiency of mixing associated with the clouds’ supersonic turbulence can be constrained by the observed metallicity distribution of the stars formed within them. In the Pleiades cluster, stars which emerged from the same molecular cloud have nearly identical metal abundance W02 . This astronomical context therefore imposes a strong motivation for a general analysis of the transport and mixing processes in compressible supersonically turbulent media.

Analytical and numerical studies of diffusion processes are typically restricted to certain families of statistical processes, like random walk MK00  or remapping models or certain Hamiltonian systems I92 . The direct numerical modeling of turbulent physical flows mostly concentrates on incompressible media (e.g. DS97 ; MKM99 ; OL01 ), but some studies have been extended into the weakly compressible regime CKM95 ; HCB95 ; PPW92 ; PPW94 ; PPSW99 . Although highly compressible supersonic turbulent flows have been studied in several specific astrophysical contexts PW00 ; SPWHW00 ; SOG98 ; OGS99 ; OSG01 ; M99 ; KHM00 ; SMH00 ; HMK01 ; HZMLN01 ; BM02 ; PVP95 ; PV98 ; VPP95 ; VS92 ; BVS99 ; BP99 ; BCP01 ; PN99 ; PZN00 ; BNP02 ; GPPL01 , the diffusion properties of such flows have not been investigated in detail.

It is the goal of this paper to analyze transport phenomena in supersonic compressible turbulent flows and to demonstrate that – analogous to the incompressible case – a simple mixing length description can be found even for strongly supersonic and highly compressible turbulence. We first briefly recapitulate in Section II the Taylor formalism for describing the efficiency of turbulent diffusion in subsonic flows. In Section III we describe the numerical method which we use to integrate the Navier-Stokes equation. In Section IV we report the diffusion coefficient obtained in our numerical models, and in Section V we introduce an extension of the well known mixing length approach to diffusion into the supersonic compressible regime. Finally, in Section VI we summarize our results.

## II A statistical description of turbulent diffusion

Transport properties in fluids and gases can be characterized by studying the time evolution of the second central moment of some representative fluid-elements’ displacement in the medium,

$$\xi _\stackrel{}{r}^2(tt^{})=[\stackrel{}{r}_i(t)\stackrel{}{r}_i(t^{})]^2_i,$$
(1)
where the average $`_i`$ is taken over an ensemble of passively advected tracer particles $`i`$ (e.g. dye in a fluid, or smoke in air) that are placed in the medium at a time $`t^{}`$ at positions $`\stackrel{}{r}_i(t^{})`$; or where the average is taken over the fluid molecules themselves (or equivalently, over sufficiently small and distinguishable fluid elements). The dispersion in one spatial direction, say along the $`x`$-coordinate, is $`\xi _x^2(tt^{})=[x_i(t)x_i(t^{})]^2_i`$. For isotropic turbulence it follows that $`\xi _x^2=\xi _y^2=\xi _z^2=1/3\xi _\stackrel{}{r}^2`$. For fully-developed stationary turbulence, the initial time $`t^{}`$ can be chosen at random and for simplicity is set to zero in what follows.

The quantity $`\xi _\stackrel{}{r}(t)`$ can be associated with the diffusion coefficient $`D`$ as derived for the classical diffusion equation,

$$\frac{n}{t}=D\stackrel{}{}^2n,$$
(2)
where $`n(\stackrel{}{r}_i,t)`$ is the probability distribution function (pdf) for finding a particle $`i`$ at position $`\stackrel{}{r}_i(t)`$ at time $`t`$ when it initially was at a location $`\stackrel{}{r}_i(0)`$. This holds if the particle position is a random variable with a Gaussian distribution B49 . In the classical sense, $`n(\stackrel{}{r},t)`$ may correspond to the contaminant density in the medium. Equation 2 holds for normal diffusion processes and for time scales larger than the typical particles’ correlation time scale $`\tau `$.

In general, however, the Lagrangian diffusion coefficient is time dependent and can be defined as

$$D(t)=\frac{d\xi _\stackrel{}{r}^2(t)}{dt}=2\stackrel{}{r}_i(t)\stackrel{}{v}_i(t)_i,$$
(3)
where $`\stackrel{}{v}_i(t)=d\stackrel{}{r}_i(t)/dt`$ is the Lagrangian velocity of the particle. The diffusion coefficient along one spatial direction, say along the $`x`$-coordinate, follows accordingly as $`D_x=d\xi _x^2(t)/dt=2x_i(t)v_{x.i}(t)_i`$. Equation 3 holds for homogeneous turbulence with zero mean velocity. From $`\stackrel{}{r}_i(t)=\stackrel{}{r}_i(0)+_0^t\stackrel{}{v}_i(t^{})𝑑t^{}`$ it follows that

$$D(t)=2\left[\stackrel{}{r}_i(0)+_0^t\stackrel{}{v}_i(t^{})𝑑t^{}\right]\stackrel{}{v}_i(t)_i=2_0^t\stackrel{}{v}_i(t^{})\stackrel{}{v}_i(t)_i𝑑t^{}.$$
(4)
The above expression allows us to related $`D(t)`$ to the trace of the Lagrangian velocity autocorrelation tensor $`\mathrm{tr}𝒞(tt^{})=\stackrel{}{v}_i(t^{})\stackrel{}{v}_i(t)_i`$ as

$$D(t)=2_0^t\mathrm{tr}𝒞(tt^{})𝑑t^{}=2_0^t\mathrm{tr}𝒞(t^{})𝑑t^{},$$
(5)
a result which was derived by Taylor already 1921 T21 . This formulation has the advantage that it is fully general and that it allows us to study anomalous diffusion processes. Note, that strictly speaking any transport process with $`\xi _\stackrel{}{r}(t)`$ not growing linearly in time is called anomalous diffusion. This is always the case for time intervals shorter than the correlation time $`\tau `$, but sometimes anomalous diffusion can also occur for $`t\tau `$. If $`\xi _\stackrel{}{r}(t)t^\alpha `$ and if $`\alpha <1`$ transport processes are called subdiffusive, if $`\alpha >1`$ they are called superdiffusive Lesieur ; I92 ; CMMV99 ; LM00 . Studying transport processes directly in terms of the particle displacement, i.e. Equation 1, is useful when attempting to find simple approximations to the diffusion coefficient $`D(t)`$ for example in a mixing length approach.

## III Numerical Method

In order to utilize the above formalism, we carry out a series of numerical simulation of supersonic turbulent flows. A variety of numerical schemes can be used to describe the time evolution of gases and fluids. By far the most widely-used and thoroughly-studied class of methods is based on the finite difference representations of the equations of hydrodynamics (e.g. P77 ). In the most simple implementation, the fluid properties are calculated on equidistant spatially fixed grid points in a Cartesian coordinate system. Finite difference schemes have well defined mathematical convergence properties, and can be generalized to very complex, time varying, non-equidistant meshes with arbitrary geometrical properties. However, it is very difficult to obtain a Lagrangian description, which is essential when dealing with compressible supersonic turbulence with a high degree of vorticity. Methods that do not rely on any kind of mesh representation at all are therefore highly desirable.

For the current investigation we use smoothed particle hydrodynamics (SPH), which is a fully Lagrangian, particle-based method to solve the equations of hydrodynamics. The fluid is represented by an ensemble of particles, where flow properties and thermodynamic observables are obtained as local averages from a kernel smoothing procedure (typically based on cubic spline functions) B90 ; M92 . Each particle $`i`$ is characterized by mass $`m_i`$, velocity $`\stackrel{}{v}_i`$ and position $`\stackrel{}{r}_i`$ and carries in addition density $`\rho _i`$, internal energy $`ϵ_i`$ or temperature $`T_i`$, and pressure $`p_i`$. The SPH method is commonly used in the astrophysics community because it can resolve large density contrasts simply by increasing the particle concentration in regions where it is needed. This versatility is important for handling compressible turbulent flows where density fluctuations will occur at random places and random times. The same scheme that allows for high spatial resolution in high-density regions, however, delivers only limited spatial resolution in low-density regions. There, the number density of SPH particles is small and thus the volume necessary to obtain a meaningful local average tends to be large. Furthermore, SPH requires the introduction of a von Neumann Richtmyer artificial viscosity to prevent interparticle penetration, shock fronts are thus smeared out over two to three local smoothing lengths. Altogether, the performance and convergence properties of the method are well understood and tested against analytic models and other numerical schemes, for example in the context of turbulent supersonic astrophysical flows PRL ; KB00 ; KB01 ; KHM00 , and its intrinsic diffusivity is sufficiently low to allow for the current investigation of turbulent diffusion phenomena L99 .

To simplify the analysis we assume the medium is infinite and isotropic on large scales, and consider a cubic volume which is subject to periodic boundary conditions. The medium is described as an ideal gas with an isothermal equation of state, i.e. pressure $`p`$ relates to the density $`\rho `$ as $`p=c_\mathrm{s}^2\rho `$ with $`c_\mathrm{s}`$ being the speed of sound. Throughout this paper we adopt normalized units, where all physical constants (like the gas constant), total mass $`M`$, mean density $`\rho `$, and the linear size $`L`$ of the cube all are set to unity. The speed of sound is $`c_\mathrm{s}=0.05`$, hence, the sound crossing time through the cube follows as $`t_{\mathrm{sound}}=20`$. In all models discussed here, the fluid is represented by an ensemble of $`\mathrm{205\hspace{0.17em}379}`$ SPH particles which gives sufficient resolution for the purpose of the current analysis.

Supersonic turbulence is known to decay rapidly PRL ; SOG98 ; PN99 ; BM99 ; MB00 ; BM00 . Stationary turbulence in the interstellar medium therefore requires a continuous energy input. To generate and maintain the turbulent flow we introduce random Gaussian forcing fields in a narrow range of wavenumbers such that the total kinetic energy contained in the system remains approximately constant. We generate the forcing field for each direction separately and simply add up the three contributions. Thus, we excite both, solenoidal as well as compressible modes at the same time. The typical ratio between the solenoidal and compressible energy component is between 2:1 and 3:1 in the resulting turbulent flow (see e.g. Figure 8 in KHM00 .). We keep the forcing field fixed in space, but adjust its amplitude in order to maintain a constant energy input rate into the system compensating for the energy loss due to dissipation (for further details on the method see M99 ; KHM00 ). This non-local driving scheme allows us to exactly control the (spatial) scale which carries the peak of the turbulent kinetic energy. It is this property that motivated our choice of random Gaussien fields as driving source. In reality the forcing of turbulence in the interstellar medium is likely to be a multi-scale phenomenon with appreciable contributions from differential rotation (i.e. shear) in the Galctic disk and energy input from supernovae explosions ending the lives of massive stars MK03 . Comparable to the values observed in interstellar gas, we study flows with Mach numbers in the range 0.5 to 10, where we define the Mach number from the one-dimensional rms velocity dispersion $`\sigma _v`$ as $`=\sigma _v/c_\mathrm{s}`$. For each value of the Mach number we consider three different cases, one case where turbulence is driven on large scales only (i.e. with wavenumbers $`k`$ in the interval $`1k2`$), intermediate-wavelength turbulence ($`3k4`$), and small-scale turbulence ($`7k8`$), as summarized in Table 1. Note that our models are not subject to global shear because of the adopted periodic boundary conditions. We call turbulence ”large scale” when the Fourier decomposition of the velocity field is dominated by the largest scales possible for the considered volume $`L^3`$, i.e. the system becomes isotropic and homogeneous only on scales larger than $`L`$. On scales below $`L`$ it may exhibit a considerable degree of anisotropy. This is most noticeable in the case $`1k2`$, because wavenumber space is very poorly sampled and variance effects become significant. The system is dominated by one or two large shock fronts that cross through the medium. In the interval $`7k8`$ the number of Fourier modes contribution to the velocity field is large, and the system appears more isotropic and homogeneous already on distances smaller than $`L`$. This trend is clearly visible in Figure 1.

Similar to any other numerical calculations, the models discussed here fall short of describing real gases in comprehensive details as they cannot include all physical processes that may act on the medium. In interstellar gas clouds, transport properties and chemical mixing will not only be determined by the compressible turbulence alone, but the density and velocity structure is also influenced by magnetic fields, chemical reactions, and radiation transfer processes. Furthermore, all numerical models are resolution limited. The turbulent inertial range in our large-scale turbulence simulations spans over about 1.5 decades in wavenumber. This range is considerably less than what is observed in interstellar gas clouds. The same limitation holds for the Reynolds numbers achieved in the models, they fall short of the values in real gas clouds by several orders of magnitude. Nevertheless, despite these obvious shortcomings, the results derived here do characterize global transport properties in interstellar gas clouds and in other supersonically turbulent compressible flows.

## IV Transport properties

Supersonic turbulence in compressible media establishes a complex network of interacting shocks. Converging shock fronts locally generate large density enhancements, diverging flows create rarefied voids of low gas density. The fluctuations in turbulent velocity fields are highly transient, as the random flow that creates local density enhancements can disperse them again. The life time of individual shock-generated clumps corresponds to the time interval for two successive shocks to pass through the same location in space, which in turn depends on the length scale of turbulence and on the Mach number of the flow.

The velocity field of turbulence that is driven at large wavelengths is found to be dominated by large-scale shocks which are very efficient in sweeping up material, thus creating massive coherent density structures. The shock passing time is rather long, and shock-generated clumps can travel quite some distance before begin disrupted. On the contrary, when energy is inserted mainly on small scales, the network of interacting shocks is very tightly knit. Clumps have low masses and the time interval between two shock fronts passing through the same location is small, hence, swept-up gas cannot travel far before being dispersed again.

The density and velocity structure of three models with large-, intermediate-, and small-wavelength turbulence is visualized in Figure 1. It shows cuts through the centers of the simulated volume. As turbulence is stationary, all times are equivalent, and the snapshot in the upper panel is taken at some arbitrary time. The lower panel depicts the system some time interval later corresponding to $`1/4`$ shock crossing time through the cube. One clearly notices markable differences in the density and velocity field between the three models.

### IV.1 Transport properties in an absolute reference frame

In order to drive supersonic turbulence and to maintain a given rms Mach number in the flow, we use a random Gaussian velocity field with zero mean to ‘agitate’ the fluid elements at each timestep. However, despite the fact that the driving scheme has zero mean, the system is likely to experience a net acceleration and develop an appreciable drift velocity, because of the compressibility of the medium. This evolutionary trend is well illustrated in Figure 2 which plots the time evolution of the three components of the mean velocity for models 2$`\mathrm{}`$, 2i, and 2s, with rms Mach numbers $`3.1`$, where turbulence is driven on (a) large (i.e. with small wavenumbers $`1k2`$), (b) intermediate ($`3k4`$), and (c) small scales (with $`7k8`$). The net acceleration is most pronounced when turbulent energy is inserted on the global scales, as in this case larger and more coherent velocity gradients can build up across the volume compared to small-scale turbulence.

The tendency for the zero-mean Gaussian driving mechanisms to induce significant center-of-mass drift velocities in highly compressible media can be understood as follows. Suppose the gas is perturbed by one single mode in form of a sine wave. If the medium is homogeneous and incompressible, equal amounts of mass would be accelerated in the forward as well as in the backward direction. But, if the medium is inhomogeneous, there would be an imbalance between the two directions and the result would be a net acceleration of the system. If the density distribution remains fixed, this acceleration would be compensated by an equal amount of deceleration after half a period, and the center of mass would simply oscillate. However, if the system is highly compressible and the driving field is a superposition of plane waves, the density distribution would change continuously (and randomly). Any net acceleration at one instance in time would not be completely compensated after some finite time interval later. This will only occur for $`t\mathrm{}`$ assuming ergodicity of the flow. Subsequently, the system is explected to develop a net flow velocity in some random direction for $`t<\mathrm{}`$. This effect is most clearly noticeable for long-wavelength turbulence, where density and velocity structure is dominated by the coherent large-scale structure. But the effect is small for turbulence that is excited on small scales, because in this limit, there is a large number of accelerated ‘cells’ which in turn compensate for another’s acceleration.

The property that the compressible turbulent flows are likely to pick up average drift velocities, even when driven by Gaussian fields with zero mean, has implications for the transport coefficients. Figure 3 shows the time evolution of the absolute (Eulerian) diffusion coefficients $`D_x`$, $`D_y`$, and $`D_z`$ in each spatial direction computed from Equation 1. Note that for stationary turbulence, only time differences are relevant and one is free to chose the initial time. In order to improve the statistical significance of the analysis, we obtain $`D(t)`$ and $`\xi _\stackrel{}{r}(t)`$ by further averaging over all time intervals $`t`$ that ‘fit into’ the full timespan of the simulation.

Due to the (continuous) net acceleration experienced by the system, the quantity $`\xi _\stackrel{}{r}^2(t)`$ grows faster than linearly with time, even for intervals much larger than the correlation time $`\tau `$, i.e. for $`\tau t<\mathrm{}`$. The system resides in a superdiffusive regime, where $`D(t)`$ does not saturate. Instead, $`D(t)`$ grows continuously with time, which is most evident in model 2$`\mathrm{}`$ of large-scale turbulence. The ever increasing drift velocity $`\stackrel{}{v}_i(t)_i`$ causes strong velocity correlations leading continuous growth of the velocity autocorrelation tensor $`_0^t\mathrm{tr}𝒞(t^{})𝑑t^{}`$. This net motion, however, can be corrected for, allowing us to study the dispersion of particles in a reference frame that moves along with the average flow velocity of the system.

### IV.2 Transport properties in flow coordinates

In order to gain further insight into the transport properties of compressible supersonic turbulent flows, we study the time evolution of the relative (Lagrangian) diffusion coefficient. In this prescription,

$$D^{}(t)=\frac{d\xi _\stackrel{}{r}^2(t)}{dt}$$

is obtained relative to a frame of reference which comoves with the mean motion of the system $`\stackrel{}{v}_i(t)_i`$ on the trajectory $`\stackrel{}{r}_i(t)_i=_0^t\stackrel{}{v}_i(t^{})_i𝑑t^{}`$. Then, $`\xi _\stackrel{}{r}^2(tt^{})=[(\stackrel{}{r}_i(t)\stackrel{}{r}_i(t)_i)(\stackrel{}{r}_i(t^{})\stackrel{}{r}_i(t^{})_i)]^2_i`$ (see Equation 1).

In Figure 4, we show the evolution of $`D^{}(t)`$ for each coordinate direction for the complete suite of models. The rms Mach numbers range from about 0.5 to 10, each considered for three cases where turbulence is driven on large, intermediate, and small scales, respectively. The plots are rescaled such that the time-averaged one-dimensional rms velocity dispersion $`\overline{\sigma }_v`$ is normalized to unity (for each direction separately). We also rescale the time $`t`$ with respect to the average shock crossing time scale through the computational volume, $`t_{\mathrm{cross}}=L/\overline{\sigma }_v`$. Recall that $`L=1`$, and note that $`\overline{\sigma }_v`$ usually differs between the three spatial directions because of the variance effects, especially in models of large-scale turbulence.

In Figure 4, we demonstrate that the magnitude of $`D^{}(t)`$ saturates for large time intervals $`t>\tau `$ in all directions. In a reference frame that follows the mean motion of the flow, diffusion in compressible supersonically turbulent media indeed behaves in a normal manner. For small time intervals $`t<\tau `$, however, the system still exhibits an anomalous diffusion even with the mean-motion correction. In this regime $`D^{}(t)`$ grows roughly linearly with time. For $`t>\tau `$ the diffusion coefficient $`D^{}(t)`$ reaches an asymptotic limit. This result holds for the entire range of Mach numbers studied and for turbulence that is maintained by energy input on very different spatial scales.

From Figure 4, we find that diffusion in compressible supersonic turbulent flows follows a universal law. It can be obtained by using the rms Mach number (together with the sound speed $`c_\mathrm{s}`$) as characterizing parameter for rescaling the velocity dispersion $`\sigma _v`$, and the rms shock crossing time scale through the volume $`t_{\mathrm{cross}}=L/(c_\mathrm{s})`$ for rescaling the time. The normalized diffusion coefficient $`D^{}(t)`$ exhibits a universal slope of two at times $`t<\tau `$ (i.e. in the superdiffusive regime), and approaches a constant value that depends only on the length scale but not on the strength (i.e. the resulting Mach number) of the mechanism that drives the turbulence. Even for highly compressible supersonic turbulent flows it is possible to find simple scaling relations to characterize the transport properties — analogous to the mixing length description of diffusive processes in incompressible subsonic turbulent flows.

## V A mixing length description

Incompressible turbulence is often described in terms of a hierarchy of turbulent eddies, where each eddy contains multiple eddies of smaller size on the lower levels of the hierarchy, while itself being part of turbulent eddy at larger scales R22 ; K41 ; O41 . At each level of the hierarchy, an eddy is characterized by a typical lengthscale $`\stackrel{~}{\mathrm{}}`$ and a typical velocity $`\stackrel{~}{v}`$. The typical lifetime of an eddy is its ‘turn-over’ time $`\tau =\stackrel{~}{\mathrm{}}/\stackrel{~}{v}`$. This mixing length prescription is an attempt to characterize the flow properties in terms of the typical scales $`\stackrel{~}{\mathrm{}}`$ and $`\stackrel{~}{v}`$. For example, this classical picture defines an effective ‘eddy’ viscosity $`\mu =\rho \stackrel{~}{\mathrm{}}\stackrel{~}{v}`$, where $`\rho `$ is the density. The mixing length $`\stackrel{~}{\mathrm{}}`$ is interpreted to be the turbulent analogue of the mean free path of molecules in the kinetic theory of gases, with $`\stackrel{~}{v}`$ being the characteristic velocity of the turbulent fluctuation.

In such a model, the velocities of gas molecules within an eddy are strongly correlated within a time interval $`t<\tau `$. They all follow the eddy rotation; the diffusion process is coherent. However, for $`t\tau `$ the velocities of gas molecules become uncorrelated, as the eddy has long been destroyed and dispersed. Hence, the velocity autocorrelation function vanishes for large time intervals, $`𝒞(t)0`$ for $`t\mathrm{}`$. Diffusion becomes incoherent as in Brownian motion or the random walk. The diffusion coefficient in the mixing length approach simply is $`D(t)2\stackrel{~}{v}^2t`$ in the regime $`t<\tau `$, which follows from replacing $`\stackrel{}{r}(t)`$ by $`\stackrel{~}{v}t`$ and $`\stackrel{}{v}(t)`$ by $`\stackrel{~}{v}`$ in Equation 3. As the largest correlation length is the eddy size, $`\stackrel{}{r}(t)`$ is substituted by $`\stackrel{~}{\mathrm{}}=\stackrel{~}{v}\tau `$ for times $`t\tau `$, and the classical mixing length theory yields $`D(t)2\stackrel{~}{\mathrm{}}\stackrel{~}{v}=2\stackrel{~}{v}^2\tau =\mathrm{constant}`$.

Compressible, supersonic, turbulent flows rapidly build up a network of interacting shocks with highly transient density and velocity structure. Density fluctuations are generated by locally converging flows, and their lifetimes are determined by the time $`\tau `$ between two successive shock passages. This time interval is determined by the typical shock velocity, which is roughly the rms velocity of the flow, i.e. the Mach number times the sound speed, $`\sigma _v=c_\mathrm{s}`$. It also depends on the length scale at which energy is inserted into the system to maintain the turbulence, which in our case is $`L/k`$ with $`k`$ being the driving wavenumber and $`L`$ being the size of the considered region (recall that in our models $`L`$ is unity). This length scale is also the typical traveling distance before two shocks interact with each other. As basic ingredients for a supersonic compressible mixing length description we can thus identify:

$`\mathrm{shock}\mathrm{travel}\mathrm{length}:`$ $`\stackrel{~}{\mathrm{}}L/k,`$ (6)
$`\mathrm{rms}\mathrm{velocity}:`$ $`\stackrel{~}{v}\sigma _v=c_\mathrm{s}.`$ (7)
The Lagrangian velocity correlation time scale, $`\tau `$, is analogeous to the time interval during which shock-generated density fluctuation remains unperturbed and moves coherently before it is being dispersed by the interaction with a new shock front. This time interval is equivalent to the time scale a shock travels along its ‘mean free path’ $`\stackrel{~}{\mathrm{}}`$ with an rms velocity $`\stackrel{~}{v}`$. This crossing time is $`\tau =\stackrel{~}{\mathrm{}}\stackrel{~}{v}\sigma _vL/k`$. For $`t<\tau `$ gas molecules can travel coherently within individual shock generated density fluctuations, and the diffusion coefficient in the mixing length prescription follows as

$$D^{}(t)2\stackrel{~}{v}^2t2\sigma _v^2t.$$
(8)
$`D^{}(t)`$ grows linearly with time with slope $`2\sigma _v^2`$. For large times, $`t\tau `$, $`D^{}(t)`$ approaches a constant value,

$$D^{}(t)2\stackrel{~}{v}^2\tau 2\sigma _vL/k.$$
(9)
This mixing length approach (Equations 8 and 9) suggests a unique scaling dependence of the diffusion coefficients in supersonic compressible flows on the Mach number $``$ and on the length scale $`\stackrel{~}{\mathrm{}}`$ of the most energy containing modes with respect to the total size $`L`$ of the system considered.

We can use $``$ (together with the given value of the sound speed) to normalize the rms velocity: $`\sigma _v=c_\mathrm{s}\sigma _v^{}=1`$. And we can also rescale the time with respect to the rms shock crossing time scale through the total volume, which is $`t_{\mathrm{cross}}=L/\sigma _v=L/(c_\mathrm{s})=t_{\mathrm{sound}}/`$ with $`t_{\mathrm{sound}}=L/c_\mathrm{s}`$ being the sound crossing time, so that $`tt^{}=t/t_{\mathrm{cross}}`$. From this normalization procedure, we get $`D^{}(t)D^{\prime \prime }(t^{})=D^{}(t)c_\mathrm{s}L`$ and obtain the following universal profile for the diffusion coefficient,

$`D^{\prime \prime }(t^{})=2t^{}`$ $`\mathrm{for}`$ $`t^{}1/k\mathrm{and}`$ (10)
$`D^{\prime \prime }(t^{})=2/k`$ $`\mathrm{for}`$ $`t^{}1/k.`$ (11)
Note that this result holds for each velocity component separately, as the results in Figure 4 indicate. In this case $`\sigma _v`$ stands for $`\sigma _x`$, $`\sigma _y`$, or $`\sigma _z`$ in Equations 8 and 9, and it hold for the total diffusion coefficient, when using $`\sigma _v=(\sigma _x^2+\sigma _y^2+\sigma _z^2)^{1/2}`$ instead.

The validity of the mixing length approximation is quantified in Figure 5 which plots the mixing length predictions against the values obtained from the numerical models. For large-scale and intermediate scale turbulence, the mixing length approach gives very satisfying results, only for small-scale turbulence it underestimates the diffusion strength. This disparity probably has to do with the numerical resolution of the code, in the sense that driving wavenumbers of $`k8`$ come close to the dissipation scale of the method and hence the inertial range of turbulence is limited KHM00 . That limitation leads to an effective driving for the turbulent motion on somewhat larger scales than $`1/8`$. Consequently, it leads to a stronger diffusion, i.e. somewhat larger diffusion coefficients than those predicted by Equations 10 and 11. The same numerical effects also account for the slightly shallower slope of $`D^{}(t)`$ for $`t\tau `$ for models $`7k8`$.

Figures 4 and 5 indicate that the classical mixing length theory can be extended from incompressible (subsonic) turbulence into the regime of supersonic turbulence of highly compressible media. In this case, driving length $`\stackrel{~}{\mathrm{}}`$ and rms velocity dispersion $`\sigma _v=c_\mathrm{s}`$ act as characteristic length and velocity scales in the mixing length approach. Note, that this only applies to mean-motion corrected transport. In general (i.e. in an absolute reference frame), supersonic turbulence in compressible media leads to superdiffusion as visualized in Figure 3.

## VI Summary

We studied diffusion processes in supersonically turbulent, compressible media. To drive turbulence and maintain the desired rms Mach number in the flow, we insert energy into the system at a pre-specified rate and over a given spatial scale using random Gaussian fields. In our numerical models, the adopted magnitude of the rms Mach numbers range from $`=0.5`$ to $`=10`$, and turbulence was driven on large, intermediate, and small scales, respectively.

Supersonic turbulence in compressible media establishes a complex network of interacting shocks. Converging shock fronts locally generate large density enhancements, diverging flows create voids of low gas density. The fluctuations in turbulent velocity fields are highly transient, as the random flow that creates local density enhancements can disperse them again.

Due to compressibility, supersonically turbulent flows will usually develop noticeable drift velocities, especially when turbulence is driven on large scales, even when it is excited with Gaussian fields with zero mean. This tendency has consequences for the transport properties in an absolute reference frame. The flow exhibits super-diffusive behavior (see also B01 ). However, when the diffusion process is analyzed in a comoving coordinate system, i.e. when the induced bulk motion is being corrected, the system exhibits normal behavior. The diffusion coefficient $`D(t)`$ saturates for large time intervals, $`t\mathrm{}`$.

By extending classical mixing length theory into the supersonic regime we propose a simple description for the diffusion coefficient based on the rms velocity $`\stackrel{~}{v}`$ of the flow and the typical shock travel distance $`\stackrel{~}{\mathrm{}}`$,

$`D^{}(t)=2\stackrel{~}{v}^2t`$ for $`t\stackrel{~}{\mathrm{}}/\stackrel{~}{v},`$
$`D^{}(t)=2\stackrel{~}{v}\stackrel{~}{\mathrm{}}`$ for $`t\stackrel{~}{\mathrm{}}/\stackrel{~}{v}.`$
This functional form may be used in those numerical models where knowledge of the mixing properties of turbulent supersonic flows is required, but where these flows cannot be adequately resolved. This is the case, for example, in astrophysical simulations of galaxy formation and evolution, where the chemical enrichment of the interstellar gas and the distribution and spreading of heavy elements produced from massive stars throughout galactic disks needs to be treated without being able to follow the turbulent motion of interstellar gas on small enough scales relevant to star formation R91 ; B98 . Our results furthermore are directly relevant for understanding the properties of individual star-forming interstellar gas clouds within the disk of our Milky Way. These are dominated by supersonic turbulent motions which can provide support against gravitational collapse on global scales, while at the same time produce localized density enhancements that allow for collapse, and thus stellar birth, on small scales. The efficiency and timescale of star formation in galactic gas clouds depend on the intricate interplay between their internal gravitational attraction and their turbulent energy contentMK03 . The same is true for the statistical properties of the resulting star clusters. For example, the element abundances in young stellar clusters are found to be very homogeneous W02 , implying that the gas out of which these stars formed must be have been chemically well mixed initially. On the basis of the results discussed here, this observation can be used to constrain astrophysical models of interstellar turbulence in star-forming regions. Understanding transport processes and element mixing in supersonic turbulent flows thus is a prerequisite for gaining deeper insight into the star formation phenomenon in our Galaxy.

###### Acknowledgements.

We thank Javier Ballesteros-Paredes, Peter Bodenheimer, Mordecai-Mark Mac Low, and Enrique Vázquez-Semadeni for many stimulating discussions. RSK acknowledges support by the Emmy Noether Program of the Deutsche Forschungsgemeinschaft (DFG, KL1358/1) and subsidies through a NASA astrophysics theory program at the joint Center for Star Formation Studies at NASA-Ames Research Center, UC Berkeley, and UC Santa Cruz.


format: text

text:
# AXPS AND X-RAY DIM NEUTRON STARS: RECENT XMM-NEWTON AND CHANDRA RESULTS

## INTRODUCTION

ROSAT observations revealed the existence of a group of soft X-ray sources with well defined and unique characteristics. To date we know seven or eight such objects which are summarized in Table 1. Their energy spectra are well modeled by blackbody emission with kTBB between 50 and 120 eV and no significant flux changes are observed on time scales of years. No association with a supernova remnant was found and upper limits derived from radio surveys indicate that the radio luminosity is well below that of ordinary radio pulsars (for recent reviews see Treves et al. 2000, Motch 2001). X-ray pulsations were discovered from two sources in ROSAT data and recent Chandra and XMM-Newton observations identified two further stars as pulsars. Low hydrogen column densities derived from the X-ray spectra indicate that we deal with close and therefore low-luminosity objects (X-ray dim isolated neutron stars, XDINs).

Although it is now well established that the ROSAT sources are isolated neutron stars, the origin of their soft X-ray emission is not clear. The temperatures in the range of (0.6 – 1.4)$`\times 10^6`$ K inferred from their blackbody-like spectra are typical for young cooling neutron stars with ages less than 106 years. However these young rotation powered pulsars show radio emission and neutron star spin periods well below 1 s, both properties in contrast to that of the new ROSAT sources. The long periods would require that these neutron stars were born with either an unusually long period or an unusually strong magnetic field to decelerate the rotation of the neutron star to the current value within its life time (Heyl and Hernquist, 1998). A decaying strong magnetic field could provide a significant source for the X-ray luminosity which, at least in the case of RX J0720.4–3125, is too high to be explained by rotational energy losses. Heating by magnetic field decay could also keep the high B field pulsars detectable in X-rays over longer times than ordinary pulsars (Heyl and Kulkarni, 1998).

Alternatively, old ($`>10^9`$ years) neutron stars could be re-heated by accretion from the interstellar medium, as originally proposed by Ostriker, Rees and Silk (1970). Population synthesis studies by Treves and Colpi (1991), Blaes and Madau (1993) and Madau and Blaes (1994) estimated the number of old, isolated neutron stars detectable in the ROSAT all-sky survey to several thousands. Despite intensive searches however, it became clear that the large galactic population of old neutron stars does not accrete from the interstellar medium at the originally predicted rates (Neuhäuser and Trümper, 1999; Motch 2001). Two possible explanations for the elusiveness of old isolated neutron stars were elaborated by Treves et al. (2000). The velocity distribution may peak at $``$200–400 km s-1, strongly suppressing spherical accretion, or magnetic field decay on timescales of 108-9 years could prevent a large fraction of neutron stars from entering the accretor stage. The first evidence for a high space velocity of 185 km s-1 was found for RX J1856.5–3754, the brightest and possibly nearest (117 $`\pm `$ 12 pc derived from parallax measurements) of the new isolated neutron stars (Walter and Lattimer, 2002).

A crucial neutron star parameter for both models is the magnetic field strength and its temporal evolution. The study of spin histories and measuring spin-down rates can provide estimates for the magnetic field strength, when interpreted in terms of magnetic dipole braking. First $`\dot{P}`$ measurements for RX J0720.4–3125 yield a magnetic field strength of the order of 1013 G (Kaplan et al., 2002a; Zane et al., 2002). Although, this value is below those typically quoted for magnetars (1014-15 G) it is still extreme and close to the critical value Bc = 4.4$`\times `$1013 G at which quantum effects become important. Anomalous X-ray pulsars (AXPs, for a recent review see Mereghetti et al., 2002) and soft $`\gamma `$-ray repeaters (SGRs, for a review see Hurley, 2000), which share some common characteristics (Kaspi, 2003), may belong to this class of neutron stars with the strongest known magnetic field strengths. It is remarkable that the spin periods of AXPs, SGRs and XDINS cluster within a narrow range. To further understand the properties of XDINs and their stellar evolution it is therefore crucial to observe them with new generation instruments available on board of Chandra and XMM-Newton. Recent results obtained from high resolution spectroscopy and timing measurements are reviewed in the following.

## MAIN PROPERTIES OF X-RAY DIM ISOLATED NEUTRON STARS

### Optical Identifications

Optical counterparts were identified for the two brightest XDINs RX J1856.5–3754 (Walter and Matthews, 1997) and RX J0720.4–3125 (Motch and Haberl, 1998; Kulkarni and van Kerkwijk, 1998). Their X-ray to optical flux ratios are log(fx/fopt) = 4.8 and 5.3, respectively. For 1RXS J130848.6+212708 Kaplan et al. (2002b) found a likely optical counterpart with log(fx/fopt) = 4.9 which would be even larger if the faint star is not the counterpart. For all other XDINS empty X-ray error circles yield lower limits for log(fx/fopt) in the range 2.5 for RX J1836.2+5925 to 4.0 in the case of RX J1605.3+3249 (see the discovery papers of each individual source given as first reference in the comment column of Table 1). These limits practically rule out any other possibility than an isolated neutron star.

XMM-Newton and in particular Chandra observations will allow to strongly reduce the uncertainties in the X-ray positions, which is highly demanded to optical identification programs of such faint stars. An example is the likely identification of 1RXS J130848.6+212708 from above which is based on a 1 arcsec error radius obtained by Chandra. A first accurate X-ray position for RX J0806.4$``$4123 was recently reported by Haberl et al. (2002).

### X-ray Spectra

The X-ray spectra obtained from most of the XDINs by the ROSAT PSPC were consistent with little absorbed blackbody emission. However the low energy resolution of the PSPC did not allow to look for absorption features which may be created in a stellar atmosphere containing heavy chemical elements. However first high resolution spectra obtained by the RGS instruments on XMM-Newton of RX J0720.4$``$3125 (Paerels et al., 2001) and by the LETGS on Chandra of RX J1856.5$``$3754 (Burwitz et al., 2001, 2003) did not reveal any absorption features. The absence of significant absorption features in the RGS spectra of RX J0720.4$``$3125 in the 0.35–1.25 keV range would appear to exclude magnetic fields of B$``$(0.3–2.0)$`\times `$1011 and (0.5–2.0)$`\times `$1014 G which would cause electron or proton cyclotron resonances, respectively, in the observed energy band (cf. Pavlov et al., 1995; Zane et al., 2000). Similarly, the LETGS spectrum (0.15–0.82 keV) would exclude B$``$(1.3–7.0)$`\times `$1011 and (0.2–1.3)$`\times `$1014 G for RX J1856.5$``$3754. Remarkably the X-ray spectra from both XDINs are best fit with a simple blackbody model.

For some of the fainter XDINs medium resolution CCD spectra are available. RX J0806.4$``$4123 was observed by the EPIC instruments on XMM-Newton (Haberl and Zavlin, 2002) and Chandra ACIS-S spectra were published for 1RXS J130848.6+212708 by Hambaryan et al. (2002). For both sources the spectra are well modeled by blackbody spectra. Temperatures included in Table 1 are those obtained from Chandra and XMM-Newton when available (for RX J0720.4$``$3125, RX J0806.4$``$4123, 1RXS J130848.6+212708 and RX J1856.5$``$3754).

### Spectral energy distribution

The optical fluxes measured from the counterparts of XDINs are higher than one would expect from the extrapolation of the X-ray spectrum. Factors around 5 are derived (which however depend on the parameters used for the X-ray spectrum which are not always consistent between different instruments) for RX J1856.5$``$3754 (Pons et al., 2002; Burwitz et al., 2003) and RX J0720.4$``$3125 (Motch and Haberl, 1998). A similar value is found for the candidate counterpart of 1RXS J130848.6+212708 (Kaplan et al., 2002b).

Pavlov et al. (1996) and Pons et al. (2002) have demonstrated that the light element (hydrogen and helium) non-magnetic neutron star atmosphere models can not reproduce the observed energy distribution of RX J1856.5$``$3754 from the optical to the X-ray band, as they over-predict the optical flux by two orders of magnitude. Heavy element (iron and standard solar-mixture) atmosphere models (Zavlin and Pavlov, 2002) for reasonable values for the gravitational redshift do not fit the X-ray spectrum of RX J1856.5$``$3754 either (Burwitz et al., 2001). This let Pons et al. (2002) to propose a two component model to explain the energy distribution of RX J1856.5$``$3754: a blackbody component with kT = 55.3 eV and small area (R = 8.2 km for a distance of 117 pc) to mainly account for the X-ray emission and another blackbody component with kT = 20 eV and larger area (R = 15.6 km) to match the observed optical flux.

A problem with this hot polar cap interpretation is however (unless the spin axis is either parallel to the magnetic axis or the observers line of sight) the non-detection of X-ray pulsations from RX J1856.5$``$3754. From the analysis of 505 ks Chandra data Ransom et al. (2002) and Drake et al. (2002) derived upper limits of 4.5% (10 ms – 103 s) and 2.7% (10 ms – 104 s), respectively. This limit was further reduced to 1.3% (20 ms – 103 s) by Burwitz et al. (2003) from a 57 ks XMM-Newton observation. In an attempt to explain this small variability, these authors discuss a one-component model with kT = 63 eV and R = 12.0 km (distance 117 pc). To match optical and X-ray data, this model requires a low radiative efficiency at X-ray energies, which may be expected if the neutron star has a condensed matter surface (Zane 2003).

### X-ray pulsations

X-ray pulsations were discovered from two XDINs in ROSAT data: a 8.39 s period for RX J0720.4–3125 (Haberl et al. 1997) and a possible period of 22.7 s (4$`\sigma `$) from RX J0420.0$``$5022 (see Fig. 1). Hambaryan et al. (2002) reported the discovery of a 5.16 s modulation in the X-ray flux of 1RXS J130848.6+212708 observed by Chandra ACIS-S. XMM-Newton observed this XDIN in Dec. 2001 for 20 ks (Hambaryan et al. in preparation). The derived power spectrum is shown in Fig. 2. Besides the strongest peak at the previously determined period of 5.16 s another peak at half the frequency suggests that the true spin period of the neutron star is 10.31 s. The light curve folded at the period of 10.31 s is shown in Fig. 1. Observed pulsed fractions are between 6% (RX J0806.4$``$4123) and 43$`\pm `$14% (RX J0420.0$``$5022). In the case of XDINs being cooling neutron stars the strongest spin modulations are difficult to explain with temperature variations as function of magnetic colatitude which would be caused by e.g. a thermal conductivity dependence on the magnetic field vector (see Cropper et al., 2001 and references therein).

### Spectral Variations

Spectral variations with pulse phase were discovered from the 8.4 s pulsar RX J0720.4$``$3125 (Cropper et al., 2001). The (0.1–0.4 keV)/(0.4–0.8 keV) hardness ratio varies sinusoidal similar to the pulsations but is softest slightly before flux maximum. Spectral fits to phase resolved spectra suggest that the majority of the variation is caused by a change in absorption rather than a change in temperature (Cropper et al., 2001). While direct changes in photo-electric absorption might be conceivable in terms of material co-rotating in the magnetosphere, the hardness ratio modulation is more likely explained by energy dependent beaming effects. If such an energy dependence is similar to that of photo-electric absorption fitting an absorbed blackbody model would necessarily lead to variations in absorption. E.g. Burwitz et al. (2003) modeled the LETG spectrum of RX J1856.5$``$3754 with a blackbody modified by an energy dependent absorption factor of the form $`\alpha `$ = Eβ in addition to photo-electric absorption. They found a best fit with $`\beta `$=1.28$`\pm `$0.30 and a reduced column density. It could well be that in the case of RX J0720.4$``$3125 the spectral variation is caused by a phase dependent absorption factor, as it may be expected for a condensed matter surface, with varying viewing angles relative to the magnetic field. Such models would not require large temperature variations on the surface of the neutron stars and greatly reduce the related problems in explaining large pulsed fractions.

### Spin Histories

To measure spin period changes with time are of great importance to independently derive estimates for the magnetic field strength of the neutron star. In particular in the absence of absorption features in the high resolution spectra of XDINs it may be the only way to obtain this information. Measuring pulse period changes should also allow to discriminate between models of cooling and accreting XDINs, with more erratic changes expected for the latter. The best studied case is RX J0720.4$``$3125, the brightest XDIN showing X-ray pulsations. Kaplan et al. (2002a) derived the first useful upper limit for the period derivative of 3.6$`\times 10^{13}`$ s s-1 (3$`\sigma `$) from the analysis of ROSAT, SAX and Chandra data. Including XMM-Newton data Zane et al. (2002) derived possible values of 3–6$`\times 10^{14}`$ s s-1 which gives a magnetic field strength of 2$`\times 10^{13}`$ G when interpreted in terms of dipolar losses. This yields a dynamic age (P/(2$`\dot{P}`$)) of $``$2.5$`\times `$106 years consistent with the age expected from theoretical cooling curves for the measured temperature. Magnetic field decay may therefore not be required to explain the observed properties of RX J0720.4$``$3125.

A first value for $`\dot{P}`$ of (0.7–2.0) $`\times 10^{11}`$ s s-1 was measured for 1RXS J130848.6+212708 (Hambaryan et al. 2002) relying on one period measurement from Chandra and one uncertain (due to possible aliases) ROSAT value. The new period obtained from XMM-Newton data is however shorter than the Chandra value and a preliminary linear fit to the three available periods (using values twice as long as the originally reported value of 5.16 s) yields only an upper limit of $`\dot{P}<6\times 10^{12}`$ s s-1 leaving an indefinite situation for this XDIN. A combined analysis of the data to explore the P and $`\dot{P}`$ parameter space is currently in progress (Hambaryan et al., in preparation).

## DISCUSSION

The new results obtained from optical and X-ray observations more and more support the model of young cooling neutron stars for XDINs. The large proper motion of RX J1856.5$``$3754 yields a small cross section for accretion from the interstellar matter and makes the X-ray luminosity difficult to explain. Also the steady and low spin down rate observed from RX J0720.4$``$3125 and the infered high magnetic field strength is incompatible with accretion which would require B $`<`$ 1010 G (Haberl et al., 1997).

The properties of AXPs are reviewed in a recent article by Mereghetti et al. (2002). The pulsating XDINs show similar pulse periods and also pulse profiles as AXPs. Figure 3 shows the period distribution for both classes of objects. In the light of the new period of 10.31 s for 1RXS J130848.6+212708 there may be some evidence for somewhat longer periods observed from XDINs on average. However the number of known objects, in particular pulsating XDINS, is still too low to draw definite conclusions. The 10.31 s period of 1RXS J130848.6+212708 implies a double-peaked pulse profile which is similar to that of the AXPs 1E 2259+586 and 4U 0142+61. However the latter changes the relative intensity of the two pulses at energies above 4 keV, a property which may be related to the hard power-law component in the spectra of AXPs which is probably not existent in XDINs. It is worth noting here, that like many other classes of X-ray pulsars also XDINs apparently can exhibit double-peaked profiles, indicating asymmetries between the two emitting areas.

The X-ray spectra of AXPs consist of a blackbody component with typical temperatures around 400 to 600 eV, much hotter than those seen in XDINs of 50 to 120 eV. Further, no hard emission component similar to the power-law seen from AXPs has been detected yet from any of the XDINs. This may indicate that AXPs are younger neutron stars (which is also supported by associations of some of them with supernova remnants) with a larger energy reservoir. The currently favored model for AXPs is that of magnetars, neutron stars with very strong magnetic fields of 1014-15 G and the decay of the field providing the energy. The question is now how do such magnetars look like after their magnetic field has decayed to $`10^{13}`$ G, the value observed from RX J0720.4$``$3125. In other words are XDINs descendents of AXPs?

XDINs are far more numerous as we can detect them only in a small volume around the Sun (of a few hundred pc radius). The X-ray luminosities of 1034-36 erg s-1 for AXPs allow their detection within the whole Galaxy. From volume considerations maybe 103-4 more XDINs than AXPs currently exist in the Galaxy, a number however highly uncertain because it is not clear how representative the solar neighborhood is in terms of volume density of XDINs. A factor of about 100 (assuming evolutionary time scales of 106 and 104 years for XDINs and AXPs, respectively, during which we can detect these sources) could be explained by evolution of AXPs into XDINs. Hence, only 1-10% of XDINs may be descendents of AXPs, the majority being born with weaker (1013 G) magnetic fields. One such case may be RX J0720.4$``$3125, for which no significant magnetic field decay is required to explain its observed spin period history and temperature. The existence of such objects is naturally expected if one accepts the idea of magnetars with magnetic field strengths above $`10^{14}`$ G and a field strength distribution at neutron star birth down to the 1012 G typically observed from X-ray binary pulsars.

If AXPs cool down during their evolution, loose their X-ray activity and become XDINs, an important question is also why we do not see any isolated neutron star with intermediate temperature of 100–400 eV. MS 0317.7$``$6647 was the first object proposed for an old accreting neutron star on the basis of its soft X-ray spectrum: Stocke et al. (1995) presented a ROSAT PSPC spectrum and derived a blackbody temperature of 180 eV. The source is located in the field of the nearby galaxy NGC 1313 (and is also called NGC 1313 X-2) and an interpretation as ultra-luminous X-ray source in NGC 1313 could not be excluded. This latter possibility is now strongly supported by XMM-Newton observations of NGC 1313 which reveal an additional hard component in the spectrum of X-2 and its similarity to X-1, another ultra-luminous X-ray source in NGC 1313 (Miller et al., 2003). Therefore, the ROSAT discovered XDINs are probably the only XDINs we now today. However, it should be remarked here, that their low temperatures might partly be due to a selection effect. Searches in ROSAT all-sky survey data were performed among the very softest sources (Haberl et al., 1998) and the simply larger number of harder ROSAT sources makes the identification as XDIN more difficult. Evolutionary scenarios remain therefore speculation until more information about magnetic field strength of XDINs is obtained and more of these enigmatic objects are discovered.


text:
# Selective addressing of high-rank atomic polarization moments

## Abstract

We describe a method of selective generation and study of polarization moments of up to the highest rank $`\kappa =2F`$ possible for a quantum state with total angular momentum $`F`$. The technique is based on nonlinear magneto-optical rotation with frequency-modulated light. Various polarization moments are distinguished by the periodicity of light-polarization rotation induced by the atoms during Larmor precession and exhibit distinct light-intensity and frequency dependences. We apply the method to study polarization moments of 87Rb atoms contained in a vapor cell with antirelaxation coating. Distinct ultra-narrow (1-Hz wide) resonances, corresponding to different multipoles, appear in the magnetic-field dependence of the optical rotation. The use of the highest-multipole resonances has important applications in quantum and nonlinear optics and in magnetometry.

High-rank polarization moments (PM) and associated high-order coherences have recently drawn attention (see Łobodziński and Gawlik (1996, 1997); Suter et al. (1993); Suter and Marty (1994); Xu et al. (1997a, b); Budker et al. (2002a); Matsko et al. (2003) and references therein) because they may enhance nonlinear optical effects important in applications such as electromagnetically induced transparency Harris (1997), creation of nonclassical atomic and photonic states Parkins et al. (1993), realization of photon blockades Imamoglu et al. (1997), quantum gates Turchette et al. (1995), photonic switches Harris and Yamamoto (1998), and atomic magnetometry Alexandrov et al. (1997); Okunevich (2001).

While signatures of high-order PM were detected in several experiments Suter et al. (1993); Suter and Marty (1994); Xu et al. (1997a, b); Matsko et al. (2003), the methods used in these investigations are not sufficiently selective and/or do not allow real-time manipulation of particular multipoles. Here we describe a method, based on nonlinear optical rotation with frequency-modulated light (FM NMOR) Budker et al. (2002b), by which one can selectively induce, control, and study any possible multipole moment. Applying the method to 87Rb atoms in a paraffin-coated cell Alexandrov et al. (1996, 2002), we have verified the expected power and spectral dependences of the resonant signals and obtained a quantitative comparison of relaxation rates for the even-rank moments.

The density matrix in the $`M,M^{}`$ representation for a state with total angular momentum $`F`$ can be decomposed into PM of rank $`\kappa =0\mathrm{}2F`$, uncoupled under rotations, with components $`q=\kappa \mathrm{}\kappa `$:

$$\rho _q^{(\kappa )}=\underset{M,M^{}=F}{\overset{F}{}}(1)^{FM^{}}F,M,F,M^{}|\kappa ,q\rho _{M,M^{}},$$

where $`\mathrm{}|\mathrm{}`$ indicate the Clebsch-Gordan coefficients (see, for example, Varshalovich et al. (1988)). For a given choice of quantization axis, $`|q|0`$ components of the PM are related to coherences between Zeeman sublevels for which $`\mathrm{\Delta }M=MM^{}=q`$, while $`q=0`$ components depend on sublevel populations. Thus, coherences $`\rho _{M,M^{}}`$ contribute to all PM with $`\kappa |\mathrm{\Delta }M|`$ and a $`|\mathrm{\Delta }M|=2F`$ coherence is uniquely associated with the highest PM for a given state, e.g., the quadrupole moment ($`\kappa =2`$) for $`F=1`$, or the hexadecapole ($`\kappa =4`$) for $`F=2`$. The method introduced here exploits the different axial symmetries of the PM (2-fold and 4-fold for the quadrupole and hexadecapole, respectively; Fig. 1) to selectively create and detect them (see also Suter et al. (1993); Suter and Marty (1994); Xu et al. (1997a, b)).

While multipole moments of rank $`\kappa 2`$ can be easily generated and detected with weak light (since a photon has spin one), higher-rank moments require multi-photon interactions for both production *and detection*. In the present method, we use a single laser beam (which is still of sub-mW power) for the nonlinear interactions required to pump and probe the high-multipole moments.

Under conditions of our experiment, FM NMOR (Fig. 2) can be understood as a three-stage (pump, precession, probe) process: atoms are polarized in an interaction with the laser beam (whose diameter is much smaller than the vapor cell dimensions), then leave the beam and bounce around the cell while undergoing Larmor precession, and finally return into the laser beam region and undergo the “probe” interaction. The laser light is frequency modulated, causing the optical pumping and probing to acquire a periodic time dependence. When the pumping rate is synchronized with the precession of atomic polarization, a resonance occurs and the atomic medium is pumped into a polarized state which rotates around the direction of the magnetic field at the Larmor frequency $`\mathrm{\Omega }_L`$. The optical properties of the medium are modulated at the frequency $`\kappa \mathrm{\Omega }_L`$, due to the symmetry of atomic polarization with rank $`\kappa `$. For example, for the quadrupole moment the modulation is at $`2\mathrm{\Omega }_L`$, and for the hexadecapole it is at $`4\mathrm{\Omega }_L`$ \[Fig. 1(a,d)\]. This periodic change of the optical properties of the atomic vapor modulates the angle of the light polarization, leading to the FM NMOR resonances. If the time-dependent optical rotation is measured at the first harmonic of $`\mathrm{\Omega }_m`$, a resonance is seen when $`\mathrm{\Omega }_m`$ coincides with $`\kappa \mathrm{\Omega }_L`$ (Fig. 3) 111Additional resonances can be observed at higher harmonics Budker et al. (2002b)..

At the resonance for a PM of rank $`\kappa `$ (which should be absent for states with $`2F\kappa `$) we expect the signal amplitude to go as the $`\kappa `$-th power of light intensity at low intensities, as discussed below. Next, we show that these predictions are verified in this experiment.

In the present experiment (Fig. 2), we used the FM NMOR technique Budker et al. (2002b) with 87Rb atoms. The central laser frequency was tuned near various hfs components of the D1 line. The typical light power was a few hundred $`\mu `$W and the laser beam diameter was $``$$`3`$mm. The laser frequency was modulated at $`\mathrm{\Omega }_m/(2\pi )`$ from $`50`$Hz to 1 kHz, and the frequency modulation amplitude was approximately $`40`$MHz (peak-to-peak). The vapor cell, with isotopically enriched 87Rb, is 10 cm in diameter and has an antirelaxation coating and no buffer gas. The cell is surrounded with four layers of magnetic shielding. A system of coils inside the innermost shield is used to compensate the residual fields (at a level $`0.1\mu `$G) and first-order gradients, and to apply a well-controlled, arbitrarily directed magnetic field to the atoms. This allows observation of FM NMOR resonances with magnetic-field widths of about 1$`\mu `$G in the low-light-intensity limit.

Figure 3 shows the magnetic-field dependence of the observed FM NMOR signals. The central laser frequency was tuned to the low-frequency slope of the $`F=2F^{}=1`$ absorption line as shown in Fig. 4. At relatively low light power \[Fig. 3(a)\], there are three prominent resonances: one at $`B=0`$, and two corresponding to $`2\mathrm{\Omega }_L=\mathrm{\Omega }_m`$. Much smaller signals, whose relative amplitudes rapidly grow with light power, are seen at $`4\mathrm{\Omega }_L=\mathrm{\Omega }_m`$, the expected positions of the hexadecapole resonances.

The spectral dependences for both types of resonance signals with fixed magnetic field and modulation frequency are shown in Fig. 4. While the quadrupole resonance signals \[$`\mathrm{\Omega }_m=2\mathrm{\Omega }_L`$; Fig. 4(a)\] are observed for both ground-state hyperfine components, no signals are observed for the hexadecapole resonances \[$`\mathrm{\Omega }_m=4\mathrm{\Omega }_L`$; Fig. 4(b)\] near the lines involving the $`F=1`$ ground state, which can not support a hexadecapole moment 222The spectral dependences shown in Fig. 4(a,b) were obtained by subtracting the average of spectra of the quadrature component taken approximately $`15\mu `$G above and below the resonance from the spectrum recorded with magnetic field set at the resonance value. This procedure renders a measurement free from the small residual contribution from the zero-field resonance arising from the “transit” effect Budker et al. (2002b).. The different shapes of the quadrupole and hexadecapole spectra near the $`F=2`$ transition group require further analysis.

Light-intensity dependences of the resonance amplitudes are shown in Fig. 5. The observed low-intensity ($`I`$) asymptotics of these curves (which show saturation at higher intensities) scale approximately as $`I^2`$ and $`I^4`$. These power dependences and many other salient features of the observed resonances can be understood with the help of a simple model of an $`F=2F^{}=1`$ transition with independent pumping, evolution in a magnetic field, and probing, in which to lowest-order, the quadrupole and hexadecapole moments are first ($`\rho ^{(2)}I_{pump}`$) and second order ($`\rho ^{(4)}I_{pump}^2`$) in light intensity, respectively. The component of the signal due to the induced optical rotation in the probe beam proportional to the quadrupole moment goes as $`I_{probe}\rho ^{(2)}`$, while that proportional to the hexadecapole moment goes as $`I_{probe}^2\rho ^{(4)}`$. The intensity dependences predicted by this model, as well as the periodicity of the signals with respect to Larmor precession and the relative widths of the resonances, match the observations from the experiment, in which a single beam serves as pump and probe.

Figure 6 shows the dependence of the resonance widths on power at a fixed magnetic field. While both the quadrupole and hexadecapole resonances exhibit power broadening, it is much less pronounced in the latter case. This is important for applications such as magnetometry because it allows operation at higher light powers with better statistical sensitivity Budker et al. (2000a). In the zero-power limit, we find that the resonance widths (Fig. 6) for the two types of resonances tend to values near $`1\mu `$G with ratio $`\mathrm{\Delta }B_q/\mathrm{\Delta }B_h=0.94(4)`$. The expected ratio of these widths is twice the ratio of the light-independent relaxation rates for the PM with $`\kappa =2`$ and 4 333It can be shown that the width of an FM NMOR resonance corresponding to a PM of rank $`\kappa `$ goes as $`\gamma _\kappa /(\kappa g\mu )`$, where $`\gamma _\kappa `$ is the rate of light-independent relaxation for a PM of rank $`\kappa `$, $`g`$ is the gyromagnetic ratio, and $`\mu `$ is the Bohr magneton.. Thus, we find that the light-independent relaxation rate for the hexadecapole is approximately twice that of the quadrupole. Relaxation of the PM in our experimental conditions is dominated by the residual relaxation on the paraffin-coated cell walls and spin-exchange collisions between Rb atoms. The electron-randomization collision model (see, for example, Knize et al. (1988)), predicts that the quadrupole moment relaxes at a rate $`3/8`$ that of the hexadecapole moment, which relaxes at the electron-randomization rate. Thus, the observed ratio of the widths is close to the expected.

In conclusion, we have developed and applied a new technique for the study of high atomic-polarization moments, allowing their selective creation and detection via corresponding well-resolved resonances. The study of the highest-rank PM signals in FM NMOR is important for the applications of this technique in Earth-field magnetometry Budker et al. (2002b). This is because energy separation for Zeeman sublevels with $`\mathrm{\Delta }M=2F`$ of the $`F=I+1/2`$ states of the alkali atoms is linear in the magnetic field, while separations between other sublevels are generally nonlinear due to magnetic-field mixing of states of different $`F`$. Therefore, it is advantageous to use of the FM NMOR resonances due to the $`\mathrm{\Delta }M=2F`$ coherence, as will be discussed in detail elsewhere.

Future work will explore the extension of the technique to separated pump and probe light beams and application to higher angular-momentum states. This will allow further optimization of the selective control of multipole moments and will elucidate the exact mechanisms responsible for the nonlinear light-atom interactions involved at the pump and the probe stages. Specifically, we will explore the role of conversion of the high-rank moments into alignment and orientation under the combined action of the magnetic field and the light shifts, which can be of crucial importance in nonlinear optical rotation Budker et al. (2000b); Rochester and Budker (2001); Rochester et al. (2001); Budker et al. (2002a).

We thank M. Zolotorev for useful discussions and A. Vaynberg for crafting parts of the apparatus. This work has been supported by the Office of Naval Research (grant N00014-97-1-0214); by NRC US-Poland Twinning grant 015369, and by the Polish Committee for Scientific Research (grant PBZ/KBN/043/PO3/2001) and the National Laboratory of AMO Physics in Toruń, Poland; by a US-Armenian bilateral Grant CRDF AP2-3213/NFSAT PH 071-02; by NSF; and by the Director, Office of Science, Nuclear Science Division, of the U.S. Department of Energy under contract DE-AC03-76SF00098. D.B. also acknowledges the support of the Miller Institute for Basic Research in Science.