text:
# Electric-field dependent spin diffusion and spin injection into semiconductors

## Abstract

We derive a drift-diffusion equation for spin polarization in semiconductors by consistently taking into account electric-field effects and nondegenerate electron statistics. We identify a high-field diffusive regime which has no analogue in metals. In this regime there are two distinct spin diffusion lengths. Furthermore, spin injection from a ferromagnetic metal into a semiconductor is enhanced by several orders of magnitude and spins can be transported over distances much greater than the low-field spin diffusion length.

.

Semiconductor devices based on the control and manipulation of electron spin (semiconductor spintronics) have recently attracted considerable attention wolf . Spin transport and injection properties of semiconductors and heterostructures strongly constrain the design of new spintronic devices. In theoretical studies of spin transport and injection in semiconductors schmidt ; rashba ; smith  the spin polarization is usually assumed to obey the same diffusion equation as in metals son ,

$$^2(\mu _{}\mu _{})(\mu _{}\mu _{})/L^2=0,$$
(1)
where $`\mu _{()}`$ is the electrochemical potential of up-spin (down-spin) electrons. In this diffusion equation, the electric field does not play any role, and spin polarization decays away on a length scale of $`L`$ from an injection point. This is reasonable for metals because the electric field E is essentially screened. For semiconductor spintronic devices, however, the semiconductor often is lightly doped and nondegenerate, and moderate electric field can dominate the carrier motion. Equation (1) corresponds to neglecting drift in the more general drift-diffusion equation for the spin polarization,

$$^2(n_{}n_{})+\frac{e𝐄}{k_BT}(n_{}n_{})\frac{(n_{}n_{})}{(L^{(s)})^2}=0,$$
(2)
where $`n_{}n_{}`$ is the difference between up-spin and down-spin electron densities and $`L^{(s)}`$ is the intrinsic spin diffusion length.

If Eq. (1) holds, spin injection from a ferromagnetic metal to a semiconductor without a spin-selective interfacial barrier is virtually impossible due to the “conductivity mismatch”, or more precisely, a mismatch between effective resistances in the metal ($`L^{(f)}/\sigma _f`$) and in the semiconductor ($`L^{(s)}/\sigma _s`$) schmidt ; rashba ; smith . Here $`L^{(f)}`$ and $`L^{(s)}`$ are the spin diffusion lengths for the ferromagnetic metal and the semiconductor, and $`\sigma _f`$ and $`\sigma _s`$ are conductivities for the two materials. Even for spin injection from ferromagnetic semiconductors, $`L^{(f)}/\sigma _fL^{(s)}/\sigma _s`$, and the spin polarization is much less than 99%, so the large spin injection percentages achieved from ZnMnSe molenkamp ; jonker1  and GaMnAs ohno  are difficult to understand via Eq. (1).

Here we clarify the central role of the electric field on spin transport in semiconductors. We obtain the drift-diffusion equation (2) for the spin polarization in a semiconductor. Equation (2) consistently takes into account electric-field effects and nondegenerate electron statistics. We identify a high-field diffusive regime which has no analogue in metals. This regime occurs for field as small as 1 V/cm at low temperatures. Two distinct spin diffusion lengths now characterize spin motion, i.e., up-stream ($`L_u`$) and down-stream ($`L_d`$) spin diffusion lengths, which can differ in orders of magnitude with realistic fields: $`E2.5`$ V/cm at $`T=3`$ K and $`E250`$ V/cm at $`T=300`$ K. These two length scales play distinctive but both favorable roles in spin injection from a ferromagnetic metal to a semiconductor. We find that the effective semiconductor resistance determining the injection efficiency is $`L_u/\sigma _s`$ rather than $`L^{(s)}/\sigma _s`$, which may be comparable to $`L^{(f)}/\sigma _f`$ given that $`L_u`$ can be shorter than $`L^{(s)}`$ by several orders of magnitude in the high-field regime. Moreover, the decay length scale for the spin polarization injected into the semiconductor is $`L_d`$, which would be much longer than $`L^{(s)}`$ in the presence of a strong field. Our results suggest a simple and practical approach to increase spin injection and spin coherence in semiconductors, namely, increasing the electric field, or equivalently, increasing the total injection current in semiconductors. Our results are consistent with the significant current dependence observed for spin injection from Fe to GaAs jonker2 . We further note that strong fields also substantially enhance spin injection in structures with an interfacial barrier.

The semiconductor we consider here is lightly or moderately $`n`$-doped ($`p`$-doped semiconductors can be analyzed similarly), which is typical in spintronic devices. We assume that there is no space charge and the material is homogeneous. The current for up-spin and down-spin can be written as $`𝐣_{()}=\sigma _{()}𝐄+eDn_{()},`$ which consists of the drift current and the diffusion one. Here $`D`$ is the electron diffusion constant, $`\sigma _{()}`$ the up-spin (down-spin) conductivity, and $`n_{()}`$ the up-spin (down-spin) electron density. The spin-dependent conductivity is proportional to the electron density for individual spins, $`\sigma _{()}=n_{()}e\nu _e,`$ where the mobility $`\nu _e`$ is assumed to be independent of field and density. The rate at which spin-up (spin-down) electrons scatter to spin-down (spin-up) electrons is denoted by $`1/\tau _{}`$ ($`1/\tau _{}`$). In steady state, the equations of continuity for individual spins read

$`𝐣_{}`$ $`=`$ $`\sigma _{}𝐄+\sigma _{}𝐄+eD^2n_{}=\left({\displaystyle \frac{n_{}}{\tau _{}}}{\displaystyle \frac{n_{}}{\tau _{}}}\right)e,`$
$`𝐣_{}`$ $`=`$ $`\sigma _{}𝐄+\sigma _{}𝐄+eD^2n_{}=\left({\displaystyle \frac{n_{}}{\tau _{}}}{\displaystyle \frac{n_{}}{\tau _{}}}\right)e.`$
In nondegenerate semiconductors, $`\tau _{}^1=\tau _{}^1\tau ^1/2`$.

For a homogeneous semiconductor without space-charge, local variation of electron density $`\mathrm{\Delta }n`$ should be balanced by a local change of hole concentration $`\mathrm{\Delta }p`$. In doped semiconductors, spin polarization can be created without changing electrons or hole densities $`\mathrm{\Delta }n=\mathrm{\Delta }p=0`$ neutrality , and therefore,

$$\mathrm{\Delta }n_{}+\mathrm{\Delta }n_{}=0.$$
(3)
Here $`\mathrm{\Delta }n_{()}=n_{()}n_0/2`$, and $`n_0`$ is the total electron density in equilibrium. From Poisson’s equation, $`𝐄=(\mathrm{\Delta }n_{}+\mathrm{\Delta }n_{})e/ϵ=0.`$ By using the Einstein’s relation, $`D=k_BT\nu _e/e`$, where $`k_B`$ is the Boltzmann constant and $`T`$ is temperature, we obtain the differential equation (2) for $`n_{}n_{}`$, the measure of the spin polarization in semiconductors, with $`L^{(s)}=\sqrt{D\tau }`$.

Equation (2), together with the local charge neutrality constraint Eq. (3), dramatically alters the spin transport behavior in semiconductors from that expected from Eq. (1). The general form of solution to Eq. (2) (restricting variation to the $`x`$-direction) is

$$n_{}n_{}=A\mathrm{exp}(x/L_1)+B\mathrm{exp}(x/L_2),$$

where $`\lambda _1=1/L_1`$ and $`\lambda _2=1/L_2`$ are the roots of the quadratic equation, $`\lambda ^2\lambda eE/k_BT1/(L^{(s)})^2=0.`$ One of the roots of the above equation must be positive and the other negative. The choice of roots is constrained by the boundary conditions at $`\pm \mathrm{}`$. To understand the physical consequence of the electric field on the spin transport, we suppose that a continuous spin imbalance is injected at $`x=0`$, and the electric field is along the $`x`$ direction. The spin polarization will gradually decay in size as the distance from the point of injection increases. In Fig. 1, we plot the spin polarization as a function of position for different fields. In the absence of the field, as shown in Fig. 1(a), the spin polarization decays symmetrically along $`x`$ and $`+x`$ with a single length scale, $`L^{(s)}`$. When an electric field is applied, the decay of the spin polarization becomes spatially asymmetric. For spin diffusion opposite to the field direction (down-stream for electrons), the decay length of the spin polarization is longer than $`L^{(s)}`$. For spin diffusion along the field direction (up-stream for electrons), the decay length is shorter than $`L^{(s)}`$. As we change the strength of the field, the spatial distribution of the spin polarization can change dramatically.

We define two quantities $`L_d`$, $`L_u`$,

$$L_d=\left[\frac{|eE|}{2k_BT}+\sqrt{\left(\frac{eE}{2k_BT}\right)^2+\frac{1}{(L^{(s)})^2}}\right]^1,$$
(4)
$$L_u=\left[\frac{|eE|}{2k_BT}+\sqrt{\left(\frac{eE}{2k_BT}\right)^2+\frac{1}{(L^{(s)})^2}}\right]^1.$$
(5)
The distribution of the spin polarization in Fig. 1 is then described by $`n_{}n_{}\mathrm{exp}(x/L_d)`$ for $`x>0`$, and $`n_{}n_{}\mathrm{exp}(x/L_u)`$ for $`x<0`$. Thus $`L_d`$ ($`L_d>L^{(s)}`$) and $`L_u`$ ($`L_u<L^{(s)}`$) are the down-stream and up-stream spin diffusion lengths, respectively.

Figure 2 shows $`L_d`$ and $`L_u`$ as a function of the electric field. In the absence of the field, the down-stream and up-stream lengths are equal to the intrinsic diffusion length $`L^{(s)}`$. With increasing field the down-stream diffusion length $`L_d`$ increases, whereas the up-stream diffusion length $`L_u`$ decreases. A high-field regime for spin transport in semiconductors can be defined by $`E>E_c`$, where $`eE_c/k_BT=1/L^{(s)}`$. In this regime, $`L_u`$ and $`L_d`$ deviate from $`L^{(s)}`$ considerably and the spin diffusion behavior is qualitatively different from that in low fields. We emphasize that since $`L^{(s)}`$ is large in semiconductors, this regime is not beyond realistic fields where most spintronic devices operate. For a typical spin diffusion length, $`L^{(s)}=10^4`$ nm length , $`E_c=25`$ V/cm at $`T=300`$ K and $`E_c=0.25`$ V/cm at $`T=3`$ K.

The physics of the field effects on the spin diffusion becomes clearer at the strong-field limit, where $`|eE|/K_BT1/L^{(s)}`$. In this limit, the electrons move with velocity $`|E|\nu _e`$ and so does the spin polarization. $`L_d`$ is simply the distance over which the carriers move within the spin life time $`\tau `$, $`L_d|E|\nu _e\tau =|E|\frac{e}{k_BT}D\tau =(L^{(s)})^2|eE|/k_BT`$. For the up-stream diffusion length $`L_u`$ at this limit, $`L_uk_BT/|eE|`$, which simply corresponds to a Boltzmann distribution of electrons in a retarding field.

A similar field-dependent diffusion phenomenon has been observed and studied in charge transport of minority carriers in doped semiconductors book . In fact, if $`n_{}n_{}`$ is substituted by $`\mathrm{\Delta }p`$ and $`L^{(s)}`$ is regarded as the intrinsic charge diffusion length, Eq. (2) becomes the diffusion equation for the disturbance of minority carrier in $`n`$-doped semiconductors. It is known that the electric field leads to two distinct charge diffusion lengths in this case as well as a modification of carrier injection book .

As an application of our field-dependent spin transport theory, we study how the electric field affects spin injection from a ferromagnetic metal to a semiconductor. We consider a simple one-dimensional spin injection structure to elucidate the underlying physics of electric field and nondegenerate electron statistics effects. This injection structure, as shown in the inset of Fig. 3, comprises a semi-infinite metal ($`x<0`$) and a semi-infinite semiconductor ($`x>0`$). Electrons are injected from the metal to the semiconductor, and therefore, the electric field is antiparallel to the $`x`$-axis. In the ferromagnetic metal the electrochemical potentials for individual spins satisfy the equations hershfield ,

$$\frac{d^2}{dx^2}\left(\begin{array}{c}\mu _{}\\ \mu _{}\end{array}\right)=\left(\begin{array}{cc}(D_{}^f\tau _{})^1& (D_{}^f\tau _{})^1\\ (D_{}^f\tau _{})^1& (D_{}^f\tau _{})^1\end{array}\right)\left(\begin{array}{c}\mu _{}\\ \mu _{}\end{array}\right),$$

where $`D_{()}^f`$ is the up-spin (down-spin) electron diffusion constant. In metals the conductivity and the diffusion constant are related via $`\sigma _{()}^f/D_{()}^f=e^2N_{()}(E_F)`$, and $`N_{()}(E_F)`$ is the up-spin (down-spin) density of states at Fermi energy. It is readily seen that the above equations lead to Eq. (1) if $`L^{(f)}=[(D_{}^f\tau _{})^1+(D_{}^f\tau _{})^1]^{1/2}`$. The general solution can be written as

$$\frac{1}{eJ}\left(\begin{array}{c}\mu _{}\\ \mu _{}\end{array}\right)=\frac{x}{\sigma _{}^f+\sigma _{}^f}\left(\begin{array}{c}1\\ 1\end{array}\right)+C_1e^{\frac{x}{L^{(f)}}}\left(\begin{array}{c}1/\sigma _{}^f\\ 1/\sigma _{}^f\end{array}\right),$$
(6)
where $`J`$ is the total electron current, which is a constant throughout the structure in steady state. In the semiconductor, according to Eqs. (2) and (3),

$$\mathrm{\Delta }n_{}=\mathrm{\Delta }n_{}=C_2\mathrm{exp}(x/L_d),$$
(7)
and $`J=\sigma _sE`$. In order to match boundary conditions at the interface between the metal and the semiconductor, it is desirable to know the electrochemical potentials for up-spin and down-spin electrons in the semiconductor, which are related to the electron density for individual spins via

$$\mu _{()}=k_BT\mathrm{ln}\left(1+\frac{2\mathrm{\Delta }n_{()}}{n_0}\right)+eExC_0.$$
(8)
This relation can be readily derived based on the definition of the electrochemical potential in nondegenerate semiconductors $`n_{()}\mathrm{exp}[(\mu _{()}+e\psi )/k_BT]`$, where $`Ed\psi /dx`$.

The three unknown coefficients $`C_i`$ ($`i=0,1,2`$) in Eqs. (6)-(8) will be determined by the boundary conditions at the interface. For a clean and transparent interface, i.e., no spin-flip scattering at the interface and no interface resistance, both the electrochemical potential and the current for individual spins are continuous, giving rise to three independent equations: (1) $`\mu _{}(0^{})=\mu _{}(0^+)`$, (2) $`\mu _{}(0^{})=\mu _{}(0^+)`$, and (3) $`j_{}(0^{})j_{}(0^{})=j_{}(0^+)j_{}(0^+)`$. The current can be calculated using $`j_{()}=\sigma _{()}\frac{d(\mu _{()}/e)}{dx}`$.

The spin injection in the semiconductor is usually defined via the spin polarization of the current, $`\alpha (x)=[j_{}(x)j_{}(x)]/J`$, which is found to be proportional to the spin polarization of the electron density $`n_{}n_{}`$,

$$\alpha (x)=\frac{n_{}(x)n_{}(x)}{n_0}\left(1\frac{k_BT}{eEL_d}\right).$$
(9)
Thus the solution of $`n_{}n_{}`$ in Eq. (7) indicates $`\alpha (x)=\alpha _0e^{x/L_d}`$, where $`\alpha _0`$ is the spin injection efficiency. We obtain an equation for $`\alpha _0`$, noting $`1k_BT/eEL_d=k_BT/eEL_u`$,

$$\frac{2L^{(f)}(\alpha _0p_f)}{(1p_f^2)\sigma _f}=\frac{k_BT}{eE\sigma _s}\mathrm{ln}\frac{k_BT/eEL_u+\alpha _0}{k_BT/eEL_u\alpha _0},$$
(10)
where $`\sigma _f=\sigma _{}^f+\sigma _{}^f`$, and $`p_f=(\sigma _{}^f\sigma _{}^f)/\sigma _f`$ is the spin polarization in the metal. We solve Eq. (10) and plot the spin injection efficiency $`\alpha _0`$ as a function of the electric field in Fig. 3. We see that the electric field can enhance the spin injection efficiency considerably. When $`\mathrm{\Delta }n_{()}/n_01`$, i.e., small spin polarization in the semiconductor, $`\alpha (x)`$ can be expressed in a compact form,

$$\alpha (x)=\left[\frac{L^{(f)}}{(1p_f^2)\sigma _f}+\frac{L_u}{\sigma _s}\right]^1\frac{p_fL^{(f)}}{(1p_f^2)\sigma _f}e^{\frac{x}{L_d}}.$$
(11)
This remarkable expression shows that the electric-field effects on spin injection can be described in terms of the two field-induced diffusion lengths. Both diffusion lengths affect spin injection favorably but in a different manner. The up-stream length $`L_u`$ controls the relevant resistance in the semiconductor, which determines the spin injection efficiency. With increasing field this effective resistance, $`L_u/\sigma _s`$, becomes smaller, and accordingly the spin injection efficiency is enhanced. The transport distance of the injected spin polarization in the semiconductor, however, is controlled by the down-stream length $`L_d`$. As the field increases, this distance becomes longer.

We now contrast Eq. (11) with that obtained by previous calculations schmidt ; rashba ; smith  based on Eq. (1). The spin injection

$$\alpha (x)=\left[\frac{L^{(f)}}{(1p_f^2)\sigma _f}+\frac{L^{(s)}}{\sigma _s}\right]^1\frac{p_fL^{(f)}}{(1p_f^2)\sigma _f}e^{\frac{x}{L^{(s)}}}$$
(12)
is given by the zero-field result of Eq. (11). As $`L^{(f)}L^{(s)}`$ and $`\sigma _f\sigma _s`$, the effective resistance in the metal, $`L^{(f)}/\sigma _f`$, is much less than its counterpart in the semiconductor, $`L^{(s)}/\sigma _s`$. Thus Eq. (12) suggests that this resistance mismatch makes it virtually impossible to realize an appreciable spin injection from a ferromagnetic metal to a semiconductor. However, the more general description of the spin transport in semiconductors indicates that the effective semiconductor resistance to be compared with $`L^{(f)}/\sigma _f`$ should be $`L_u/\sigma _s`$ rather than $`L^{(s)}/\sigma _s`$. Since $`L_u`$ can be smaller than $`L^{(s)}`$ by orders of magnitude in the high-field regime, this “conductivity mismatch” obstacle may be overcome with the help of strong electric fields, or equivalently, large injection currents darryl . For example, if the parameters of a spin injection device are as follows, $`p_f=0.8`$, $`L^{(f)}=100`$ nm, $`L^{(s)}=10^4`$ nm, and $`\sigma _f=100\sigma _s`$, at zero field the spin injection efficiency is 0.04%, which can be increased to 4.2% at $`|eE|/k_BT=0.02`$ nm-1, which corresponds to $`|E|=50`$ V/cm, or $`|J|=50`$ A/cm2 for a typical semiconductor conductivity $`\sigma _s=1`$ ($`\mathrm{\Omega }`$ cm)-1, at $`T=3`$ K. This may explain the large spin injection percentages from ZnMnSe to ZnSe molenkamp ; jonker1  and from Fe to GaAs ploog ; jonker2 , as well as the dramatic increase in spin injection with current in Ref. jonker2 . Finally, we note that spin injection enhancement from a spin-selective interfacial barrier between the ferromagnetic metal and the semiconductor, which has been identified in the low-field regime rashba ; smith , becomes more pronounced in the high-field regime.

In summary, we have derived the drift-diffusion equation for spin polarization in a semiconductor by consistently taking into account electric-field effects and nondegenerate electron statistics. This equation provides a framework to understand spin transport in semiconductors. We have identified a high-field diffusive regime which has no analogue in metals. In this regime, there are two distinct spin diffusion lengths, i.e., the up-stream and down-stream spin diffusion lengths. The high-field description of the spin transport in semiconductors predicts that the electric field can effectively enhance spin injection from a ferromagnetic metal into a semiconductor and substantially increase the transport distance of the spin polarization in semiconductors. Our results suggest that the “conductivity mismatch” obstacle in spin injection may be overcome with the help of high field injection in the diffusive regime.

We would like to thank N. Samarth for pointing out the relevance of resistance mismatch for spin injection from ferromagnetic semiconductors. This work was supported by DARPA/ARO DAAD19-01-0490.


format: text

text:
# Time-Varying Fine-Structure Constant and Quantum SuperStrings

## Abstract

In this note we observe that Webb and co-workers’ observation of an evolving fine structure constant has since been reconfirmed. We also confirm Kuhne’s contention that this implies a cosmological constant and exhibit a cosmological mechanism for Webb’s result. We also see how this provides a mechanism to reach elementary particle scales from the Quantum Superstring scale.

00footnotetext: E-mail:birlasc@hd1.vsnl.net.in

Kuhne suggested sometime ago that the then preliminary results of Web et al that there was a slow cosmic evolution of the fine struture constant $`\alpha `$, could be explained in terms of a cosmology with a cosmological constant. Since then Web and co-workers have reconfirmed their initial observations, and the acceleration of the universe, implying a cosmological constant has also been repeatedly confirmed. In this note we confirm Kuhne’s conclusion and provide a mechanism for the variation of $`\alpha `$.
Kuhne derives the following equation,

$$\frac{\dot{\alpha }_z}{\alpha _z}=\alpha _z\frac{\dot{H}_z}{H_z},$$
(1)
from a well known equation of Teller, where $`H`$ is the Hubble constant, $`z`$ denotes the value of the parameterr at the cosmological redshift $`z`$ and the dot represents the derivative with respect to time. He then compares (1) with the finding of Web and co-workers, viz.,

$$\frac{\dot{\alpha }_z}{\alpha _z}1\times 10^5H_z.$$
(2)
to exhibit the cosmological constant.
We consider the above results in the light of fluctuational cosmology. In this cosmology, $`\sqrt{N}`$ particles are fluctuationally created, somewhat on the lines suggested by Hayakawa-and this happens within the minimum unphysical time interval, the Compton time $`\tau `$. Such a cosmology is not only consistent with observation, but it also deduces from the theory the otherwise empirical Large Number Relations of Dirac as also the mysterious empirical Weinberg relation between the pion mass and the Hubble Constant. The cosmology also predicted an ever expanding and accelerating universe, as indeed was discovered shortly thereafter and reconfirmed since (Cf. also ). In this cosmological scheme, the cosmological constant $`\mathrm{\Lambda }`$ is given by,

$$\mathrm{\Lambda }\dot{H}0(H^2)$$
(3)
In what follows, equalities are approximate equalities in the order of magnitude sense. We now observe that substitution of (3) in (1) gives

$$\frac{\dot{\alpha }_z}{\alpha _z}=\beta H_z$$
(4)
where $`\beta <\alpha _z<10^2`$.
It can immediately be seen that (4) is compatible with (2).
We finally give the derivation of (2) in the above context wherein, as the number of particles in the universe increases with time, we go from the Planck scale to the Compton scale. It is known that this Compton length, due to Zitterbewegung causes a correction to the electrostatic potential which an orbiting electron experiences, rather like the Darwin term.
Infact we have

$$\delta V=V(\stackrel{}{r}+\delta \stackrel{}{r})V(\stackrel{}{r})$$

$$=\delta r\frac{V}{r}+\frac{1}{2}\underset{ıj}{}\delta r_ı\delta r_j\frac{^2V}{r_ır_j}$$

$$0(1)\delta r^2^2V$$
(5)
Remembering that $`V=e^2/r`$ where $`r10^8cm`$, from (5) it follows that if $`\delta rl`$, the Compton wavelength then

$$\frac{\mathrm{\Delta }\alpha }{\alpha }10^5$$
(6)
where $`\mathrm{\Delta }\alpha `$ is the change in the fine structure constant from the early universe. (6) is an equivalent form of (2) (Cf.ref.), and is the result originally obtained by Webb et al (Cf.refs.).
We now consider briefly Quantum Superstrings.
One of the criticisms levelled against Quantum Superstrings (QSS) Theory is that it is not verifiable, as it deals with Planck scale phenomena. Indeed ’t Hooft has gone so far as to say that it is not even a theory. The question is, are there any effects which can be observed at scales corresponding to elementary particles? For example could we go from QSS to the Compton scale? We will now discuss briefly, exactly such a mechanism.
As is well known, at the scale of QSS we have a non commutative geometry ,

$$[x,y]0(l^2),[x,p_x]=\mathrm{}[1+0(l^2)]etc.$$
(7)
Another way of looking at (7) is that there is a generalised Uncertainity Principle in operation at this scale ( \- )

$$\mathrm{\Delta }x\frac{\mathrm{}}{\mathrm{\Delta }p}+l_p^2\frac{\mathrm{\Delta }p}{\mathrm{}}$$
(8)
$`l_p`$ being the Planck length.
The first term on the right side of equation (8) gives the usual Uncertainity relation, while the second term represents the duality effect - as we go down to the Planck scale we infact are lead to the larger scale represented by the second term (Cf.ref.).
Before proceeding further, it may be remarked that the non commutativee geometry in (7) is a manifestation of the non zero spatial extension of the strings.
We now come to a mechanism by which the Compton scale arises quite naturally in the above context. As was pointed out a long time ago by Hayakawa, the fluctuation in the mass of a typical elementary particle due to the fluctuation of the particle number, which is $`\sqrt{N},N10^{80}`$ being the number of elementary particles in the universe, is given by

$$\frac{G\sqrt{N}m^2}{c^2R}$$

Using now the Uncertainity Principle, we get,

$$(\mathrm{\Delta }mc^2)T=\frac{G\sqrt{N}m^2}{R}T=\frac{G\sqrt{N}m^2}{c}$$
(9)
$`T`$ being the age of the universe and $`R`$ its radius which equals $`cT`$. It can be easily seen that the right side of (9) equals $`\mathrm{}`$! That is, we have,

$$\mathrm{}\frac{G\sqrt{N}m^2}{c}$$
(10)
Equation (10) immediately gives us the Compton wavelength. Interestingly a similar line of reasoning leads to a fluctuational model of cosmology, which explains the many so called Large Number coincidences as also Weinberg’s apparently mysterious empirical relation between the pion mass and the Hubble Constant, from theory, and moreover predicts an ever expanding accelerating universe, as indeed has been repeatedly verified in the past few years.
We now observe that $`n\sqrt{N}10^{40}`$. In the context of (8), let us consider a secondary fluctuation $`\sqrt{n}`$ of Planck scale particles. Then the second term of (8) gives

$$\mathrm{\Delta }x=l_p^2\frac{\sqrt{n}m_pc}{\mathrm{}}=l_p\sqrt{n},$$
(11)
$`m_p`$ being the Planck mass. The right side of (11) is precisely the Compton wavelength $`l_\pi `$ of a typical elementary particle, viz., the pion. Immediately we recover therefrom two other consistent relations for the mass of the pion and its Compton time viz.,

$$m_\pi =\frac{m_p}{\sqrt{n}},t_\pi =t_p\sqrt{n}$$

As mentioned after (8), this is the well known duality effect - as we go down to the Planck scale, we end up at the Compton scale.
Thus the QSS scale and the Compton scale are related by the duality relation.
We finally make the following remark. In the light of the above considerations leading to (9) and (11), it is interesting that we can consider the energy of a Planck mass to be the gravitational energy of the fluctuational $`\sqrt{n}`$ Planck masses within the Compton length $`l_\pi `$. That is, we should have,

$$\frac{G\sqrt{n}m_p^2}{l_\pi }m_pc^2.$$

This is indeed so.


text:
# Can randomness alone tune the fractal dimension?

## I Introduction

The history of describing natural objects by geometry is as old as science itself. Traditionally, this has involved Euclidean geometry that restricts the natural boundary of our imagination to lines, circles, cones, spheres, cuboids and so on. However nature is not restricted to Euclidean shapes only. Instead, most of the natural objects we see around us are so irregular and complex in shape that they can be described as geometrically chaotic, since they are not just merely complex but often contain different degrees of complexity. In 1975, Mandelbrot introduced the concept of fractal geometry to characterize these geometric monsters quantitatively. This extended our horizon of understanding and enables to appreciate that there exists some kind of order even in these seemingly complex and disordered geometric structural patterns. Often fractals are perceived as some brightly-colored computer generated patterns. This is due to the lack of enough analytically solvable theoretical models. The coastline of an island, a river network, the structure of a cabbage or broccoli, or even the networks of nerves and blood vessels in the normal human retina can be best described as fractals. Yet, more than twenty five years after they were first introduced, there is no generally-accepted definition of a fractal. However, it is loosely defined as a shape made of parts similar to the whole, in some sense. It is typically quantified by a non-integer exponent called the fractal dimension that can uniquely characterize the geometric structure. This definition immediately confirms the existence of scale invariance, that is, the objects look the same on different scales of observation. To understand fractals, their physical origin and how they appear in nature we need to be able to model them theoretically. This is the main motivation of our present work.

The simplest way to construct an ideal fractal is by repeating a given operation over and over again deterministically. The classical Cantor set is such a simple text book example which is typically used for learning about the self-similar properties of an ideal fractal. It is created by dividing a line into $`n`$ equal pieces, removing $`(nm)`$ of the parts created and repeating the process with the $`m`$ remaining pieces ad infinitum kn.mandelbrot . However, it is well understood that in our world nothing is stationary or strictly deterministic and all natural complex objects, even intelligence and life, occur through some kind of evolution where randomness is an essential ingredient. That is, nature favors randomness and fractals in nature appear through a continuous kinetic process, whereas the classical Cantor set is discrete in time. The simplest way to incorporate these essential ingredients is by generalizing the Cantor set. For example, instead of dividing an interval into $`n`$ equal pieces, one can divide it randomly into $`n`$ pieces and throw away $`(nm)`$ of these also randomly kn.hassan1 .

The notion of random fractal has been widely used in almost every disciplines of science both theoretically and experimentally since its inception. Yet, the mechanism by which nature creates fractals and the relationship between the degree of order and the fractal dimension is poorly understood. In this article, we introduce a stochastic process that can generate a continuous spectrum of fractals (ranging from random to nonrandom) controlled by some intrinsic tuning parameters. This may be considered as a natural kinetic counterpart of the classical Cantor construction and which can be a potential candidate in order to understand essential governing rules of creating complex geometric objects. We present a generalized version of the classical Cantor set. The rules of the process can be described as follows. The process starts with an initiator of the unit interval  and a generator which divides this interval into two pieces randomly, deleting some parts from either sides of both pieces at each time step. We introduce two intrinsic parameters to determine the degree of randomness and the rate at which a given operation is to be repeated in order to create a fractal.

The construction of this generalized Cantor set is not at all pedagogical. One immediate and potential application of the present model is the kinetics of the irreversible and sequential breakup of particles that occurs in a variety of physical processes and which has important applications in science and technology. These include erosion kn.has0 , grinding and crushing of solidskn.has1 , polymer degradation and fiber length reduction kn.has2 , breakup of liquid droplets kn.has3 , etc. to name just a few. In recent years, there has been an increasing interest in studying fragmentation, allowing variations that increase the flexibility of the theory in matching conditions of real phenomena such as the extension to higher dimension kn.has4 , agglomerate erosion kn.has0 , mass loss kn.has5 , volume change kn.has6  and fragmentation-annihilation kn.has7 . The generalized Cantor construction, as we formulated it, can describe the kinetics of fragmentation with a continuous mass loss kn.hassan2 ; kn.has5 . This is relevant in the kinetics of fragmentation of particles where mass loss might occur due to evaporation, oxidation, sublimation, dissolution, melting etc., or even in the Yule-Furry process of cosmic shower theory where energy loss occurs due to collision kn.cos .

Another potential application of the present work is the first order phase transformation via the kinetics of nucleation and growth processes. The random deposition of cuts in the generalized Cantor set can be thought of as the random nucleation of size-less seeds on the substrate (metastable phase). On the other hand, the continuous removal of parts is in fact equivalent to the decrease of substrate due to the growth of seeds (stable phase). The nucleation of a stable phase usually occurs due to large thermal fluctuations in the metastable phase and the subsequent growth occurs due to the competition between the surface tension and the free energy density difference between the stable and metastable phase. The classical work of such processes has been done independently by Kolmogorov, Johnson, Mehl and Avrami (KJMA model) kn.kolmogorov  and has been extensively used to describe non-equilibrium phenomena during phase transformation kn.krapivsky . In this model the nucleation of a stable phase, within the metastable phase, occurs at random obeying the Poisson stochastic process in space. In other words nucleation centers are homogeneously distributed throughout the entire metastable phase. Once nucleated, the stable phase keeps growing isotropically until it comes into contact with another similarly growing stable phase thus forming an interface.

Interesting questions arise from the present work: (i) What is the role of the fractal dimension during the pattern formation? (ii) Is there any relation between the fractal dimension and the degree of order? (iii) What are the relevant parameters to tune the degree of order and what are their physical meaning? The present work is an attempt to answer these questions.

## II Formulation of the cut and delete process

In order to create a complex geometric objects here we give the general scheme for the formulation of the cut and delete model. To do so, we define the interval size distribution function or the concentration $`c(x,t)`$ so that $`c(x,t)dx`$ describes the number of intervals of size $`x`$ at time $`t`$ in the interval range $`x`$ and $`x+dx`$. The distribution function $`c(x,t)`$ then must have the following form of evolution equation

$$\frac{c(x,t)}{t}=\frac{c(x,t)}{t}|_{\mathrm{cut}}+\frac{c(x,t)}{t}|_{\mathrm{delete}}.$$
(1)
The first term on the right hand side represents the irreversible and sequential cut process which can be described by the following master equation

$`{\displaystyle \frac{c(x,t)}{t}}|_{\mathrm{cut}}`$ $`=`$ $`c(x,t){\displaystyle _0^x}𝑑yF(y,xy)`$ (2)
$`+`$ $`2{\displaystyle _x^{\mathrm{}}}𝑑yc(y,t)F(x,yx).`$
Here $`F(x,y)`$ is the breakup kernel describing the rules and the rate ($`a(x)=_0^x𝑑yF(y,xy)`$) at which an interval of size $`x+y`$ breaks into sizes $`x`$ and $`y`$. The two terms in Eq. (2) describe the destruction and creation of interval of size $`x`$ due to the deposition of cuts on size $`x`$ and larger than $`x`$ respectively. The second term of Eq. (1) on the other hand, represents the deletion of parts of an interval from either ends soon after a cut is deposited. The deletion of parts continues until these parts encounter another cut or the interval becomes a dust like interval; thereby stopping the loss of their masses. This deletion process can be expressed in the following differential form

$$\frac{c(x,t)}{t}|_{\mathrm{delete}}=\frac{\left(q(x)c(x,t)\right)}{x}.$$
(3)
where $`q(x)`$ describes the rate at which an interval of size $`x`$ is lost by the deletion process.

### II.1 Model I: Random scission model

We first consider here the simplest possible model of this kind, namely

$$F(x,y)=1.$$
(4)
This is typically known as the random scission model kn.ziff  as it assumes that every points of the interval are equally likely to break. That is, the present model describes the following recursive cut and delete process. We start with an interval of size $`x+y`$. At each time step, we choose a cut point with a uniform probability density and produce two smaller intervals of size $`x`$ and $`y`$. As soon as the new intervals are born, they start shrinking from either side, although the mechanism whereby they shrink would depend on the choice of $`q(x)`$. In general, q(x) is assumed to follow a power law form with respect to particle size i.e. $`q(x)x^\gamma `$. However, once the choice for the breakup kernel is made, the exponent $`\gamma `$ is strictly restrictive, due to dimensional consistency. Therefore, in the present context it is essential to choose $`q(x)=mx^2`$ where $`m`$ is a positive, real and dimensionless constant. The rate equation for the distribution function $`c(x,t)`$ then becomes

$$\frac{c(x,t)}{t}=xc(x,t)+2_x^{\mathrm{}}c(y,t)𝑑y+m\frac{x^2c(x,t)}{x}.$$
(5)
Instead of trying to solve the rate equation directly, we attempt to solve the rate equation for $`M_n(t)`$, the $`n^{\mathrm{th}}`$ moment of $`c(x,t)`$, defined as $`M_n(t)=_0^{\mathrm{}}x^nc(x,t)𝑑x`$ with $`n>0`$. Incorporating this definition into Eq. (5), we obtain

$$\frac{dM_n(t)}{dt}=\left[\frac{n^2+n(1+1/m)1/m}{(n+1)}\right]M_{n+1}(t).$$
(6)
Note that for $`m=0`$ the total mass or the first moment $`M_1(t)`$ is a conserved quantity. However, for $`m>0`$ the system violates the mass conservation due to the continuous mass removal by the deletion process. The interesting feature of the above equation is that for $`m>0`$, there are infinitely many conserved quantities depending on $`m`$. We can find these conserved quantities by simply inspecting the value of $`n=d_f`$ for which the moment $`M_{d_f}(t)`$ is independent of time which is obtained by setting $`dM_n(t)/dt=0`$. The problem then rests on solving a quadratic equation in $`n`$. The real positive root of the equation then is

$$d_f=\frac{1}{2}\left(1+\frac{1}{m}\right)+\frac{1}{2}\sqrt{\left(1+\frac{1}{m}\right)^2+\frac{4}{m}}$$
(7)
while the other root is $`\alpha =(d_f+1+\frac{1}{m})`$ thus $`\alpha <0`$. We can re-write the rate equation for $`M_n(t)`$ as

$$\frac{dM_n(t)}{dt}=\frac{(nd_f)(n\alpha )}{(n+1)}M_{n+1}(t).$$
(8)
Iterating the above equation enables us to obtain all the derivatives of $`M_n(t)`$ at $`t=0`$. Using this and the monodisperse initial condition into the Taylor series expansion of $`M_n(t)`$ about $`t=0`$ gives the explicit and exact solution for $`M_n(t)`$

$$M_n(t)={}_{2}{}^{}F_{1}^{}(nd_f,n\alpha ;n+1;t)$$
(9)
where $`{}_{2}{}^{}F_{1}^{}(a,b;c;d)`$ is the hypergeometric function. In principle, finding a proper integral transformation of Eq. (9) and comparing it with the definition of $`M_n(t)`$ would enable us to recover the exact solution for the distribution function $`c(x,t)`$ (widely known as Charlesby’s method kn.charlesby ). However, in the present work we are more interested in $`M_n(t)`$ and its asymptotic temporal behavior than in $`c(x,t)`$. The asymptotic expansion of the hypergeometric function kn.luke  reveals that the exact time dependence of $`M_n(t)`$ is $`M_n(t)t^{(nd_f)}`$. The asymptotic decay of $`M_n(t)`$ is thus linear in $`n`$ which implies that the system obeys a simple scaling. This will be shown vigorously in the following section.

#### II.1.1 Scaling theory

A closer inspection of Eq. (5) immediately reveals that there are two governing parameters $`x`$ and $`t`$ that completely describe the governed parameter $`c(x,t)`$. This is due to the fact that the dimension of the two variables $`x`$ and $`t`$ are interlocked by the choice of the breakup kernel and through the rate equation itself. That is, $`x`$, the co-factor of $`c(x,t)`$ of the first term on the right hand side of Eq. (5), must bear the dimension inverse of time $`t`$. This enables us to choose any of the two as an independent parameter while the other can be expressed in terms of this. That is, we can define $`\xi =tx`$ as a dimensionless quantity. So, if we assume $`t`$ to be the independent quantity then $`c(x,t)`$ can be expressed in terms of $`t`$ alone i.e.

$$c(x,t)t^\theta \varphi (\xi ).$$
(10)
Here $`\theta `$ is yet an undetermined exponent and must take the value for which $`t^\theta `$ bears the dimension of concentration or the distribution function $`c`$. This means that all plots of $`\frac{c}{t^\theta }`$ vs $`xt`$ for any initial condition must fall onto each other which is the hallmark for the existence of scaling typically known as data collapse formalism kn.stanley . Substituting the scaling ansatz (Eq. (10)) into the definition of $`M_n(t)`$ and demanding that $`M_{d_f}`$ be a conserved quantity, one immediately obtains $`\theta =1+d_f`$.

#### II.1.2 Fractal analysis

Associating broken objects with disorder is perhaps the most immediate and natural thought that cross our mind. Indeed, an assembly of fragmented intervals can be considered as one of the simplest yet most common example of a disordered system. The existence of scaling is already a proof that an assembly of such sequentially broken intervals has some degree of regularity. Owing to the random nature of the present model and due to the mass removal term, it is clear that when the process evolves ad infinitum, it will not create the whole set that describes a line. Rather, it will produce a subset of the full set (line). In order to measure the size of the set created in the long time limit, we define the average interval size as

$$\delta (t)=\frac{_0^{\mathrm{}}xc(x,t)𝑑x}{_0^{\mathrm{}}c(x,t)𝑑x}t^1.$$
(11)
This shows that the mean or typical interval size decreases in time algebraically. We can use this as a yard stick and count the number of the segments we need to cover the resulting set. We find that the number of yard stick $`N(\delta )`$ scales as

$$\underset{\delta 0}{lim}N(\delta )\delta ^{d_f}$$
(12)
where $`N=_0^{\mathrm{}}c(x,t)𝑑x`$. The exponent $`d_f`$ is known as the fractal dimension or the Hausdorff-Besicovitch dimension of the resulting set created by the cut and delete model. The distance between the points in space is the key to the definition of the Hausdorff-Besicovitch dimension kn.feder . Note that $`d_f`$ is the root of the quadratic equation in $`n`$ obtained by setting $`\frac{dM_n(t)}{dt}=0`$ (a condition for stationary solution). The expression for $`d_f`$, Eq. (7), reveals that as the value of $`m`$ increases, the fractal dimension decreases very sharply and in the limit $`m\mathrm{}`$, $`d_f0`$. Fig. (1) does indeed show that as $`m`$ increases the strength of the mass deletion term increases with respect to the terms responsible for describing the cut processes. This means that as $`m`$ increases, the size of the arising set decreases sharply due to the fast disappearance of its member. While as $`m0`$, $`d_f1`$, that is, we recover the full set from which the subset is derived.

### II.2 Model II: Gaussian cut and delete model

Next we treat the following quadratic type of breakup kernel

$$F(x,y)=(xy)^{b1}.$$
(13)
This is a generalized version of the previous model and it describes the following recursive cut and delete process. We start with an interval $`(x+y)`$ and at each time event we choose a cut point with a Gaussian probability density if $`b>1`$. The spread of the density of the cut around the center of the interval depends on the values of $`b`$ as we will show later. Of course we recover the previous random scission model when $`b=1`$. Our principal goal of studying this model is to understand the role of $`b`$ during the cut and delete process. In order to extract the role of $`b`$, it is essential to check whether the relative strength between the cut and delete process changes as we change the $`b`$ value. Note that the breakup rate $`a(x)`$ for this model is $`a(x)=px^{2b1}`$ where $`p=[\mathrm{\Gamma }(b)]^2/\mathrm{\Gamma }(2b)`$. It is thus clear that as $`b`$ increases so does the relative strength of the cut process with respect to the delete process if we do not tune the $`m`$ value accordingly. That is, to judge the role of $`b`$ we ought to give equal weight to all the terms of Eq. (1) so that each of them can compete on an equal footing. This can only be done if we set $`q(x)=px^{2b}`$ (i.e. $`m=p`$) and hence the relative strength between the two terms of Eq. (1) stays the same as the value of $`b`$ increases. The rate equation for $`M_n(t)`$ then is

$$\frac{dM_n}{dt}=\mathrm{\Gamma }(b)\left(\frac{[\mathrm{\Gamma }(b)]}{\mathrm{\Gamma }(2b)}(n+1)\frac{2\mathrm{\Gamma }(n+b)}{\mathrm{\Gamma }(n+2b)}\right)M_{n+2b1}.$$
(14)
The problem now rests on finding the $`n`$ value for which the moment becomes a time independent quantity. From the previous fractal analysis we learnt that the real positive root of the polynomial equation (obtained by setting $`dM_n(t)/dt=0`$) is in fact the fractal dimension of the resulting set. Applying the scaling theory based on dimensional analysis, like we did in the previous section, gives $`\theta =(1+d_f)z`$ where the kinetic exponent is $`z=1/(2b1)`$. A detailed numerical survey reveals that the fractal dimension increases monotonically with increasing $`b`$ (see Fig. (2)). In order to find the fractal dimension in the limit $`b\mathrm{}`$, we can use the Stirlings approximation in the resulting polynomial equation. In doing so we obtain the following equation

$$\mathrm{ln}[n+1]+(n1)\mathrm{ln}[2]=0$$
(15)
and find $`d_f=0.45699956`$ is the only real positive root and hence is the fractal dimension $`d_f`$. In an attempt to search for the role of $`b`$ we study the following model

$$F(x,y)=(x+y)^\gamma \delta (xy).$$
(16)
It describes that at each time an interval is chosen and the cut is placed exactly in the middle to produce two equal sized smaller intervals. It is note worthy to mention that the homogeneity index gamma only tunes the rate $`a(x)=x^\gamma /2`$ and does not play any role in controlling the location of each cut on the interval vis-a-vis the fractal dimension. The exact value of $`\gamma `$ is not important provided $`\gamma >0`$ and therefore we can choose $`\gamma =2b1`$ without loss of generality. Once again the dimensional consistency requires that $`q(x)=x^{2b}/2`$ while the factor $`1/2`$ is to account for the equal weight in the cut and delete process. Introducing this model into Eq. (1) and inserting the definition of moment we obtain the following rate equation for $`M_n(t)`$

$$\frac{dM_n(t)}{dt}=\left(\frac{(n+1)}{2}\frac{1}{2^n}\right)M_{n+2b1}(t).$$
(17)
In order to find the fractal dimension for this process, we set the numerical factor on the right-hand side of this equation equal to zero and then take the natural logarithm on both sides to obtain the $`n`$ value for which $`M_n(t)`$ is time independent. In doing so, we arrive at the same functional equation for $`n`$ as we found for the Gaussian model with $`b\mathrm{}`$, namely Eq. (15).

This is in fact the main message of the present work. It shows that the kernel $`F(x,y)=(xy)^{b1}`$ behaves exactly in the same fashion as for $`F(x,y)=(x+y)^{2b1}\delta (xy)`$. We thus find that in the limit $`b\mathrm{}`$, the resulting distribution of points is a set with fractal dimension $`D_f=0.45699956`$, which is a strictly self-similar fractal as randomness is ceased by dividing the intervals into equal pieces. We are now in a position to give a physical picture of the role played by $`b`$. First of all, the model with $`b=1`$ describes that the frequency curve of placing cuts about the size of the fragmenting particles is Poissonian in nature. Consequently, the system enjoys the maximum randomness and the corresponding fractal dimension is $`d_f=0.414213`$. For $`b>1`$, the frequency curve of placing cuts about the size of the fragmenting particles is Gaussian in nature, meaning that as the value of $`b`$ increases, the particles are increasingly more likely to break in the middle than on either end. That is, as $`b`$ increases, the variance decreases in such a manner that in the limit $`b\mathrm{}`$, the variance of the frequency curve becomes infinitely narrow, meaning a $`\delta `$-function distribution for which the fragments are broken into two equal pieces. Therefore, there is a continuous spectrum of fractal dimensions between $`b0`$ when $`d_f=0.414213`$ and $`b\mathrm{}`$ when $`d_f=0.45699956`$ (see Fig. (2)).

## III Discussion

In this work we have presented a cut and delete process that can capture the notion of the Cantor set construction. It has potential applications both in the kinetics of fragmentation and in the kinetics of nucleation and growth phenomena. We have found that in the long time limit the system creates a geometric pattern which is statistically self-similar. To quantify the resulting geometric structure, we invoke the idea of fractal geometry and show that the fractal dimension is strictly dependent on the strength of the mass removal term $`m`$. Thus by tuning this $`m`$, we obtain a fractal dimension of any value between $`0<d_f1`$. We have further generalized our model allowing a parameter $`b`$ that can tune the degree of randomness by controlling the way the cuts are deposited. This is done in an attempt to find a precise answer of how fractal dimension changes as we tune the degree of order of the resulting geometric structure. We show that when $`b=1`$ the system enjoys the maximum degree of freedom in the sense that cut points are chosen with uniform probability density and so is the delete process, since the two processes are interlocked via dimensional consistency. For $`b>1`$, at each time event a cut point is more likely to be chosen in the middle than on either end of a given interval. That is, the frequency curve of the cut density about any given interval size is Gaussian in nature and as $`b`$ increases the middle point is increasingly more likely to be the cut point than any other points in the interval. The $`b`$ value is intimately related to the variance of the probability density of cut. As $`b`$ increases, the variance of the probability density curve decreases in such a manner that in the limit $`b\mathrm{}`$, the variance of the curve reduces to an infinitesimally narrow line. The equivalence of the $`\delta `$-function and the Gaussian model in the $`b\mathrm{}`$ limit does indeed prove this. Changing the $`b`$ value effectively means we are creating an increasingly ordered set. In this work we have derived a spectrum of fractal dimensions between $`d_f=0.414213`$ for $`b=1`$ and $`d_f=0.45699956`$ when $`b\mathrm{}`$. Fig. (2) and a detailed numerical survey confirms that the fractal dimension increases monotonically with increasing order and reaches its constant value when the randomness is completely ceased.

In summary, we have discussed the scaling theory of a particular class of cut and delete process, emphasizing on the dimensional analysis. In order to make the scaling theory more meaningful than what a simple scaling can provide we have invoked the concept of fractal analysis. We have identified a new set with a wide range of subsets produced by tuning the degree of randomness only. We have quantified the size of the resulting set obtained in this way by fractal dimension and have shown that the fractal dimension increases with increasing order and reaches its maximum value when the pattern described by the set is perfectly ordered, which is contrary to some recently found results kn.brilliantov . Finally, we argue on the basis of our findings that fractal dimension, degree of order and the extent of ramifications of the arising geometric patterns are intimately connected to each other.

Acknowledgment

M. Kamrul Hassan would like to thank Alexander von Humboldt foundation for financial support under the Georg Forster Fellowship.


text:
# Event rates and timescale distributions from realistic microlensing models of M31

## 1 Introduction

Gravitational microlensing as a means to detect compact objects in the Galactic halo was first considered by Paczyński (1986), but the basic idea is much older (Einstein 1936). This suggestion was realized as results from two surveys of microlensing events towards the Large and Small Magellanic Clouds (LMC and SMC) (Alcock et al. 2000; Ansari et al. 1996). These were consistent with a significant, but subdominant, contribution of microlensing masses to the Galactic dark matter halo. Nonetheless, these conclusions are still controversial, and the identity and location of the microlensing masses are still mysterious.

A decade ago, M31 was suggested as a promising venue where galactic microlensing might be explored in ways advantageous and distinctive from that in and around the Galaxy (Crotts 1992). Several papers (Jetzer 1994; Han & Gould 1996; Baltz & Silk 2000; Kerins et al. 2001) have confirmed that a substantial microlensing signal can be expected. Two collaborations, MEGA (preceded by the VATT/Columbia survey) and AGAPE have produced a number of microlensing event candidates involving stars in M31 (Crotts & Tomaney 1996, Ansari et al. 1999, Aurière et al. 2001, Uglesich 2001, Calchi Novati et al. 2002). Here we show the potential for these and future surveys to settle some of the outstanding questions regarding microlensing in spiral galaxies.

This is the second paper in a series. Paper I (Gyuk & Crotts 2000) provides optical depth maps for M31. While these are useful tools for certain purposes they are unfortunately not directly measurable. Event rates on the other hand are directly measurable and hence their magnitude and variation across the face of M31 are more meaningful in planning and evaluating surveys of microlensing. In this paper we extend Paper I to include event rate maps, both total and also differential rates with respect to the event timescale.

This paper is organized as follows. In §2 we briefly discuss the M31 models we used, including disk, bulge and halo components. Following this we present rate maps for various halo models, including the self lensing contribution, in §3. In §4 we provide differential rate distributions as a function of two different timescale measures, and we discuss how cuts in timescale can be used to separate self lensing from a MACHO halo contribution to the lensing rate. We conclude with a discussion of measuring lens masses and halo properties in §5.

## 2 Modeling

The execution of a microlensing survey of M31 is qualitatively different from those towards the Magellanic Clouds or our own Galactic bulge. For the latter, individual stars can often be resolved, and their unlensed fluxes measured, admitting a direct measurement of the magnification during a microlensing event. In contrast, M31 is more than ten times farther away, and individual stars are rarely resolved from the ground. Images are almost always highly crowded, and individual sources are almost always highly blended. This is the “pixel lensing” regime (Crotts 1992; Baillon et al. 1993). Image subtraction, necessitated by the high degree of crowding, has been very successful. While this technique allows the detection of variable objects of any kind, and is photon noise limited, the unlensed fluxes of lensed sources are unknown. Thus M31 microlensing surveys must make do with a full width at half maximum timescale and a flux increase at maximum as the useful fit parameters. The unlensed flux measured by classical microlensing surveys is unavailable under most circumstances.

Sources are taken to reside in a luminous two-component model of M31 consisting of a double exponential disk and a bulge. The disk model is inclined at an angle of 77 and has a scale radius of 5.98 kpc, a scale height of 400 pc, a central surface brightness of $`\mu _R=20`$, $`VR=0.67`$, and we ignore spiral arm effects (Walterbos & Kennicutt 1988; Gould 1994). We take an $`I`$ band surface brightness profile from the $`VR`$ color and the observation that $`VRRI`$ (Bessel 1979). The rotation velocity is taken to be 240 km s-1, with a linear rise from the center to 6 kpc. The bulge model is based on the “small bulge” of Kent (1989) with a central surface brightness of $`\mu _R=14`$, and $`RI=0.72`$ (Walterbos & Kennicutt 1988; Bessel 1979). This is an axisymmetric bulge with a roughly $`\mathrm{exp}(r^{0.4})`$ falloff in volume density with an effective radius of approximately 1 kpc and axis ratio, $`c/a0.8`$. Values of the bulge density are normalized to make $`M_{\mathrm{bulge}}=4\times 10^{10}M_{}`$. The velocity distribution of bulge stars is taken to be maxwellian ($`dN\mathrm{exp}(v^2/2\sigma ^2)d^3\stackrel{}{v}`$), with $`\sigma =150`$ km s-1. These quantities are fairly well known, and unlikely to change the results by a large amount if revised.

We explore a parameterized set of M31 halo models. Each model halo is an axisymmetric cored “isothermal sphere” determined by two parameters: the axis ratio or flattening $`q`$ (where $`q=1`$ indicates no flattening) and the core radius $`r_c`$. As we are concerned with lensing objects we also define the MACHO fraction ($`f_b`$) as the fraction of the halo mass that consists of lensing objects. Together, the mass density of halo lensing objects is then

$$\rho (x,y,z)=f_b\frac{V_c(\mathrm{})^2}{4\pi G}\frac{e/(q\mathrm{sin}^1e)}{x^2+y^2+(z/q)^2+r_c^2},$$
(1)
where $`e=\sqrt{1q^2}`$ and $`V_c(\mathrm{})=240`$ km s-1 is taken from observations of the M31 disk. The velocity distribution of the halo is taken to be maxwellian, with a circular velocity equal to $`V_c(\mathrm{})`$, making $`\sigma =170`$ km s-1.

We compute rate distributions with respect to two different timescales. The first is the more familiar Einstein time $`(t_E)`$, defined as the time to cross the full Einstein disk which has radius

$$R_E=\sqrt{\frac{4GM_{\mathrm{lens}}D_sx(1x)}{c^2}},$$
(2)
where $`x`$ is the fractional distance to the lens in relation to the distance to the source, $`D_s`$. The Einstein time is then given in terms of the lens velocity perpendicular to the line of sight,

$$t_E=\frac{2R_E}{v_{}}.$$
(3)
The second timescale we use is the full–width at half maximum of the lightcurve $`t_{1/2}`$, much more easily measured when the source fields are crowded, as is the case for M31. This timescale is simply related to the Einstein time and to the minimum impact parameter $`\beta `$ of the lens relative to the line of sight (taken in units of $`R_E`$). This timescale is given by (Gondolo 1999)

$$t_{1/2}=t_Ew(\beta ),$$
(4)
with the following definitions

$$w(\beta )=\sqrt{2f(f(\beta ^2))\beta ^2},f(x)=\frac{2+x}{\sqrt{x(4+x)}}1.$$
(5)
The function $`w(\beta )\beta `$ for all values of $`\beta `$, with the limiting behavior

$$w(\beta 1)\beta \sqrt{3},w(\beta 1)\beta \sqrt{\sqrt{2}1}.$$
(6)
Thus, $`t_{1/2}\beta t_E`$, and is more degenerate than $`t_E`$ due to the dependence on impact parameter. A further complication is that $`t_{1/2}`$ is essentially the measured timescale. Determining the Einstein time requires knowledge of the magnification, which can be very difficult in highly crowded fields as the source stars are highly blended. The Einstein time might be inferred with extra effort, either by high resolution imaging of the source star or statistically with the $`t_\sigma `$ technique (Baltz & Silk 2000). It is crucial to evaluate if these are necessary.

## 3 Rate maps

Using the model of §2, we compute the rate of detectable microlensing events. Two computer codes written entirely independently and described elsewhere (Baltz & Silk 2000; Gyuk & Crotts 2000) have been used, and produce nearly identical results. A numerical integration of the rate is performed, over the positions and velocities of the source (with fixed brightness) and lens (with fixed mass) along a given line of sight. The probability that an event with given parameters is detected is folded into the the integral. Mass functions for the lenses and luminosity functions for the sources are applied at the end.

With lines of sight spaced at one arcminute intervals, we construct contour maps of the event rate for “self lensing” by stars (assuming no binary lenses, at most a 15% contribution (Mao & Paczyński 1991; Baltz & Gondolo 2001)) as well as for halo lensing. The self lensing contribution has four logical contributions, taking sources and lenses from the bulge and disk. The bulge-bulge contribution dominates the self lensing near the center of the bulge, along the minor axis the bulge-disk and disk-bulge contributions dominate, and far from the bulge the disk-disk contribution dominates (though this is always dominated by the halo contribution for $`f_b>0.05`$). These four in sum give a self lensing rate that is nearly symmetric about both the major and minor axes of M31. The halo lensing contributions arise from both the M31 and Milky Way halos, lensing both disk and bulge stars. The Milky Way contribution has a nearly uniform optical depth. The M31 contribution is strongly asymmetric, with a significantly larger rate from the far edge of the minor axis.

We have assumed the following definition for a detectable microlensing event. The MEGA survey most frequently employs the MDM 2.4m telescope, thus we use its capabilities in the following. We assume that integrations totaling three hours are taken twice weekly during the M31 observing season: this is fairly conservative. We define an event as a deviation that has two consecutive samples four standard deviations above the baseline. We assume one arcsecond seeing, which is typical: the median MDM seeing is 0.95 arcsecond. To approximate the sensitivity of MDM, we assume that a star of $`R`$ magnitude $`m_R=25.2`$ or $`I`$ magnitude $`m_I=24.8`$ gives one photoelectron per second, and furthermore that the noise in the images is twice the photon counting noise (this is quite conservative). We take the distance modulus to M31 as $`D=24.5`$, a distance of 795 kpc. We assume a sky brightness of $`\mu _R=21`$ and $`\mu _I=20`$ mag arcsec-2. Here we reiterate the fact that the source stars are not resolved, but instead are typically highly blended. Only with difference imaging can the microlens variability be detected. The only measured “baseline flux” is the light falling within a resolution element due to several stars in the highly crowded fields.

The luminosity function of M31 sources must be specified. For the disk, we take $`R`$ and $`I`$ band luminosity functions from Mamon & Soneira (1982). For the bulge, we take the $`I`$ band from Terndrup, Frogel & Whitford (1990). Data for the $`R`$ band is scarce, so we average (in $`\mathrm{log}dN/dM`$) $`V`$ and $`I`$ band data for $`M_R>0`$ (using Terndrup et al. (1990) for $`I`$ and Holtzman et al. (1998) for $`V`$). For $`M_R<0`$ we attach a power law slope of 0.59 taken from the MACHO project data (Alves 2001).

The mass function of lenses, both stellar and the MACHO component, is somewhat problematic. For the stellar component, we use the exponential Chabrier (2001) mass function, down to $`0.01M_{}`$. This is steeply decreasing, so there is little mass in the lowest decade (the brown dwarfs). Varying the stellar mass function in acceptable ranges has little effect unless for example there is a large component of brown dwarfs (Baltz & Silk 2001). The mean stellar mass of about $`0.5M_{}`$ is what dominates the rate and timescale from self lensing. We have also investigated a Scalo (1986) mass function, which has a slightly higher mean mass, thus longer timescales and lower rates. The differences are not very large however. For lack of convincing evidence to the contrary, we take a fixed mass for MACHOs. Based on Alcock et al. (2000) we assume a mass of $`3/8`$ dex relative to solar $`(0.422M_{})`$, and we test values of $`0.1M_{}`$ and $`1.0M_{}`$ as well. We note here that MACHO mass can have a large effect on the expected timescales and rates. Larger masses would indicate a lower rate and longer events, while smaller masses give a higher rate but shorter events. This is discussed further in the next section.

In Figs. 1-6 we plot event rate contours for $`R`$ band observations. The contours for $`I`$ band are very similar. The only pronounced difference is in the bulge, where the $`I`$ band luminosity function is quite shallow, there is a significantly larger rate in the most central 4 square arcminutes. While it is desirable to detect events in two bands to test that the flux enhancement has a constant color (in time), as should be the case for gravitational lensing, we will not include this criterion. The separate event rates are very similar, thus we deem it unlikely that the joint event rate will be much different. In Fig. 1 we plot the rate contours for self lensing and for M31 halos with a 20% MACHO fraction c.f. Alcock et al. (2000) ($`f_b=0.2`$) and a core radius of 2 kpc. Both round ($`q=1`$) and flattened ($`q=0.3`$) halos are illustrated. These maps are given as contours of constant rate, in units of events yr-1 arcmin-2. In Fig. 2 we plot rate contours again for $`f_b=0.2`$, but compare the small ($`r_c=1`$ kpc) and large ($`r_c=5`$ kpc) core cases, for both round and flattened halos. In Fig. 3 we illustrate the event rate from Milky Way MACHOs, assuming both round and flattened cases again. The rate should be related to the surface brightness, also illustrated. In fact a naive calculation indicates that the rate from a constant optical depth of lenses (such as a Milky Way population, which wouldn’t vary much over the M31 fields) should be proportional to the square root of the surface brightness. This can be seen from the following simple argument. The number of monitored sources is proportional to the surface brightness (the proportionality constant depends on luminosity function), but the noise level increases as the square root of the surface brightness, necessitating an increase in magnification by the same factor to obtain an equivalent signal to noise. Since the peak magnification of the event is proportional to the inverse of the impact parameter, the total cross section goes down by this square root of surface brightness. Thus, the number of monitored sources times the cross section is proportional to the square root of the surface brightness. This argument breaks down at very high surface brightness where the magnification required implies events whose full width at half maximum timescales are too short to detect, however, this regime is not reached for M31 for the masses assumed. From Fig. 3 we see that this calculation is reasonable, as the shape of rate contours matches the shape of the surface brightness contours, and the rate of decline is roughly half on the logarithmic scale illustrated. The flattening of the Milky Way halo only changes the normalization of the rates, not the shape of the contours.

Combining all of these components, we illustrate the total expected event rate. In the bottom panel of Fig. 4, we plot contours of the total microlensing rate, from self lensing and 20% MACHO halos for both the Milky Way and M31. The self lensing dominates in the inner 5 kpc, but outside, the 20% halo provides most of the events. Changing the halo fraction $`f_b`$ has little effect on the shape of the contours, and only in the region where the self lensing and halo lensing are comparable (roughly speaking within the 0.3 events arcmin-2 yr-1 contour assuming a 20% halo). The overall normalization of the event rate away from the bulge can give a clear measurement of $`f_b`$. For example, comparing a 20% halo to the case of a 5% halo, the inner contours (within 1 event arcmin-2 yr-1) are not much affected since self lensing dominates and the shape of the outer contours is very similar since it is the halo that dominates there in both cases. Between the 1.0 and 0.3 contours the event rate drops more quickly in the 5% halo case since this is the region where the microlensing rate from the halo and stellar components are comparable. This is clear from inspecting the top panel of Fig. 1 (the self lensing rate) and comparing to the bottom panel (the M31 halo lensing rate) which is trivially rescaled according to MACHO fraction.

Fig. 4 shows the difference in the rate contours for round $`q=1`$ and flattened $`q=0.3`$ halos. The two cases appear to be easily distinguishable, especially along the major axis. In the bottom panel of Fig. 5, we again plot contours of the total microlensing rate, this time varying the halo core radius from $`1,2,5,10`$ kpc. The significant difference is on the far side, along the minor axis. A crude measurement of core radius may thus be possible.

Here we see the utility of using the raw event rates to determine parameters. We do not need to evaluate the optical depth to do so. Since in the end we are interested in the halo parameters like MACHO fraction and flattening, the optical depth is secondary.

## 4 Timescale distributions

The Einstein time is the fundamental timescale parameterizing the variability due to lensing, but it is notoriously difficult to measure from a pixel microlensing event. The directly measured parameter is the full-width at half-maximum timescale. We will illustrate the rate distribution of events in both of these timescales in Fig. 4, and we will discuss the extraction of halo parameters using both. The calculation of the differential rate with respect to $`t_{1/2}`$ follows Baltz & Silk (2000), based on Griest (1991) who discusses the differential rate with respect to $`t_E`$. We plot the differential rate with respect to the logarithm of the timescale, normalized by the differential rate at a fixed timescale $`\stackrel{~}{t}`$,

$$\frac{d\stackrel{~}{\mathrm{\Gamma }}}{d\mathrm{log}t}=\frac{d\mathrm{\Gamma }}{d\mathrm{log}t}/\frac{d\mathrm{\Gamma }}{d\mathrm{log}t}(\stackrel{~}{t}).$$
(7)
Taking $`\stackrel{~}{t}`$ longer than the the peak timescale then illustrates the variation in how rapidly the differential event rate falls at long timescales. In the upper panels of Fig. 4 we illustrate $`d\stackrel{~}{\mathrm{\Gamma }}/d\mathrm{log}t`$ for both $`t_E`$ and $`t_{1/2}`$. For lines of sight far from the bulge, the distribution in $`t_{1/2}`$ peaks at around 20 days, while the distribution of $`t_E`$ peaks at around 100 days. Near the bulge, the peaks of both timescale distributions are shorter.

Away from the bulge, where MACHO events dominate, the peak timescale is a measurement of the MACHO mass, as $`t_E\sqrt{M_{\mathrm{lens}}}`$. Taking a single mass, the $`t_E`$ distributions are $`0.5`$ dex wide in $`t_E`$, translating to $`1.0`$ dex wide in $`M_{\mathrm{lens}}`$. With 50 events, a $`\pm 20\%`$ measurement of the mass would be possible with known Einstein times. With only $`t_{1/2}`$, the situation is considerably worse, as the distribution in $`t_{1/2}`$ for a fixed mass is $`1.0`$ dex wide in $`t_{1/2}`$, thus $`2.0`$ dex wide in $`M_{\mathrm{lens}}`$. With the same 50 events, a mass measurement at the level of $`\pm 40\%`$ can be made. We thus see the utility of determining the Einstein times for the events.

We expect MACHO events to have longer timescales than stellar events simply due to the geometry, i.e. the source–lens distances can be much larger. To compare, the bulge is of order 1 kpc in size, while a cored halo does not drop appreciably in density for of order 10 kpc. Halo lenses have a larger velocity dispersion than bulge lenses, but not by a large factor. In all, we might expect halo lens timescales to be a factor of 2.5 or so larger. This fact can be used to determine cuts in timescale, either $`t_E`$ or $`t_{1/2}`$, that maximizes the possibility that MACHO events can be distinguished from stellar events. We want a cut that both excludes as many self lensing events as possible and includes as many MACHO events as possible. We thus wish to maximize the signal to noise, including both the Poisson counting of events, and a systematic error due to the uncertainty in the self lensing model. For a given area, observation time and minimum timescale, we expect $`N_{\mathrm{M31}}`$ events from M31 MACHOs, $`N_{\mathrm{MW}}`$ events from Milky Way MACHOs, and finally $`N_\mathrm{S}`$ events from stellar lenses, with an uncertainty of $`\mathrm{\Delta }N_\mathrm{S}`$ (which we take to be $`0.3N_\mathrm{S}`$). We minimize the quantity

$$\frac{1}{Q^2}=\frac{N_{\mathrm{M31}}+N_{\mathrm{MW}}+N_\mathrm{S}}{N_{\mathrm{M31}}^2}+\left(\frac{\mathrm{\Delta }N_\mathrm{S}}{N_{\mathrm{M31}}}\right)^2$$
(8)
as a function of minimum timescale to determine what, if any, timescale cut to use. These results are plotted in the top panels of Fig. 5. For most of the field, we find that no timescale cut should be made. Near the center of the bulge, within 6 arcmin on the minor axis and 12 arcmin on the major axis, the situation is different. If the Einstein times are known, it is advantageous to cut all events with $`t_E<75`$ days. This allows the largest signal to noise for separating the halo component. Without the cut, for a one year survey the separation can be done at the $`2.5\sigma `$ level over the interior region. With the cut, the significance rises to 4$`\sigma `$. If only $`t_{1/2}`$ is known, a similar timescale cut does not help much. We comment here that if we were to take a more conservative value for the model uncertainty (say $`\mathrm{\Delta }N_S=0.5N_S`$), the results are only affected near the center of the bulge, but the desired cutoff in Einstein time and the signal to noise for detecting a halo knowing the Einstein times are hardly affected even there. We also note that the utility of timescale cuts may change if the MACHO mass were to be significantly different from typical stellar masses. As it is now, the MACHO project mass value and the typical stellar mass are quite similar.

Over the full far–side MDM field, a microlensing halo for M31 can be detected at roughly $`7\sigma `$ significance with one full year of well–sampled data, or approximately three observing seasons. This assumes $`f_b=20\%`$, but the significance is roughly linear with $`f_b`$. For $`f_b=10\%`$, a microlensing halo can be detected at roughly $`4\sigma `$ significance, and with $`f_b=5\%`$, the significance drops to about $`2.5\sigma `$. Thus, $`f_b=5\%`$ is about the limit of detectability for a microlensing halo with three seasons of observing. This 5% value is also a plausible lower limit to halo lensing, caused by a much greater number of stars in M31’s halo/spheroid than the Galaxy’s (Reitzel et al. 1998, scaled to produce a microlensing prediction using Alcock et al. 1997, Table 10).

Finally, we discuss the effects of varying the MACHO mass on the total and differential rates. The stellar mass function will remain that of Chabrier (2001). In the bottom panel of Fig. 6 rate contours are plotted for fixed MACHO masses (for both the M31 and Milky Way contributions) of $`1.0M_{}`$, $`3/8`$ dex solar, and $`0.1M_{}`$. As expected, the lower masses give larger rates. We have also given the MACHOs a Chabrier (2001) function, identical to the stars. The resulting rates are similar to the $`0.1M_{}`$ case (as the Chabrier mass function peaks there), thus we do not illustrate them. In the top panels of Fig. 6 we show $`d\stackrel{~}{\mathrm{\Gamma }}/d\mathrm{log}t_{1/2}`$ for lens masses of $`1.0M_{}`$ and $`0.1M_{}`$. The shift in peak is clearly seen, as is the higher rate for smaller lenses. Taking a Chabrier (2001) mass function, the rate distribution is very similar to the fiducial case illustrated in Fig 4.

## 5 Discussion

We have provided maps of the microlensing event rate towards M31 in a model including M31 disk and bulge sources, and lenses from the M31 disk and bulge, as well as lenses in the dark halos of M31 and the Milky Way. We have varied the parameters of the dark halo, namely the core radius and the flattening, and studied the effects on the rate contours. The core radius affects the rates most in the galaxy center but these changes are partially obscured by the high rate of self lensing events there. The flattening does have a significant effect, especially along the major axis. The normalization of the rate contours away from the bulge gives a measure of the halo MACHO fraction.

The easily measured timescale parameter is the full width at half maximum, but the more physical timescale is the Einstein time. Knowing the Einstein times allows a much more accurate measurement of the lens mass, thus it is worth the extra trouble to try to measure the Einstein times of the detected events.

We have shown that a microlensing halo in M31 should be clearly distinguishable from self lensing if an appreciable event rate away from the M31 bulge is measured. We have explored the use of cuts in event timescale to separate the self lensing component from the halo lensing component of the event rate, and found that this helps near the bulge, but not further away. We have quantified the level of halo that is detectable, and found that a marginal detection of a 5% microlensing halo would be possible in three seasons of ground–based observations. Higher halo fractions can be detected more convincingly of course.

We wish to thank David Alves, Krzysztof Stanek and Larry Widrow for useful conversations. E.B. acknowledges support from the Columbia University Academic Quality Fund. G.G. wishes to acknowledge financial support from the Brinson Foundation. A.C. was supported by grants from NSF (AST 00-70882 and 98-02984) and STScI (GO-7376).


text:
# Convective Wavelength Shifts in the Spectra of Late-Type Stars

## 1. Introduction

Convection in late-type stars penetrates into the photosphere, producing inhomogeneities. Temperature variations of up to a few hundred Kelvin are apparent in optical images of the Sun, and a velocity field of several kilometers per second is observed in spatially resolved solar spectra. The warmer upflows appear as bright granules surrounded by narrower, cooler downflows. While present technology cannot resolve stellar granulation, there is evidence of its presence in all late-type stars. The temperature and velocity inhomogeneities produce absorption line profiles which are shifted and asymmetric, effects that are readily noticed in ultra-high dispersion stellar spectra.

An extensive literature exists on solar line asymmetries and shifts (see, e.g., Neckel & Labs 1990; Asplund et al. 2000; Pierce & Lopresto 2000). Modern stellar studies were pioneered by Dravins (1974, 1985) and Gray (1980, 1981). Dravins (1999) provides a short review of recent work. Line asymmetries and line shifts are consequences of the same phenomenon and, therefore, give largely equivalent information. When the number of lines measured is limited, but the available features are mostly unblended, line asymmetries carry more information. On the contrary, if the spectrum is heavily crowded, line shifts, which can be measured for many lines, are likely to provide more information. A limited wavelength coverage, largely dictated by the detector size and the scarcity of high-accuracy laboratory wavelengths, steered stellar studies toward the analysis of line asymmetries.

Line bisectors – a measure of the line asymmetry – vary smoothly across the HR diagram (Gray & Toner 1986; Dravins 1987; Gray & Nagel 1989; Dravins & Nordlund 1990). Other parameters, such as rotation, chemical composition, magnetic fields, or binarity, are likely to play a role, but the limited data available have prevented an in-depth study. A recent comprehensive work on the Fe i spectrum (Nave et al. 1994) has provided accurate laboratory wavelengths for 9501 lines. These data and improvements in astronomical instrumentation make it feasible and practical to turn to line shifts as a complement to line bisectors.

Gravitational and convective wavelength shifts systematically affect absolute determinations of radial velocities. Other effects can also introduce systematic radial velocity errors, but they are expected to be less important for most stars (Lindegren, Dravins, & Madsen 1999). Gravitational redshifts are proportional to the mass-to-radius ratio. For spectral lines formed in the photospheres of dwarf stars, the gravitational shift is the same for all lines, in the range between 0.4 and 2 km s-1. When accurate parallaxes and photometry are available, comparison with evolutionary models allows us to determine dwarf masses and radii within 8 % and 6 % , respectively (Allende Prieto & Lambert 1999). It is therefore possible to estimate the gravitational shift within 10 %, yielding worst-case uncertainties of $`0.2`$ km s-1. Convective line shifts may vary with spectral type. The velocity difference between weak and strong lines is about $`0.6`$ km s-1 for the Sun (Allende Prieto & García López 1998) and probably more than 1 km s-1 for the F5 subgiant $`\alpha `$ CMi (Procyon; Allende Prieto et al. 2002), but no systematic study across the HR diagram has been published. Obviously, accurate spectroscopic studies of absolute radial velocities cannot afford to neglect these shifts.

In this paper, we analyze optical spectra of several late-type stars obtained with the High Resolution Spectrograph (HRS) coupled to the Hobby-Eberly Telescope (HET). The high quality and large spectral coverage of the spectra allow us to measure line shifts for a large number of Fe i lines. Section 2 describes the observations and the data. In Section 3 we describe the analysis, and in §4 we discuss the results.

## 2. Observations

Our observations were obtained with the 9.2m HET at McDonald Observatory and the HRS (Tull 1998). The HRS saw first light on March 17, 2001, starting science observations shortly afterwards. This fiber-coupled echelle spectrograph uses an R-4 echelle, cross-dispersing gratings, and a mosaic of two 2k$`\times `$4k CCDs. In the selected mode, HRS provides a FWHM resolving power $`R120,000`$ and complete spectral coverage between 4094 and 7890 Å, with the exception of the gap (5977–5999 Å) between the CCDs. The observations were performed in queue mode in May and June 2001 (see Table 1). An observation consisted of two exposures that were extracted and combined. The final signal-to-noise (S/N) per pixel was $`400`$ at most wavelengths.

The spectra were processed with the tasks in the IRAF echelle package. We removed the bias signal, flatfielded and extracted the spectra, subtracting the low background in the wide interorder regions. In addition to the program stars observed with the HET, we have included in this study the spectra of the Sun (Kurucz et al. 1984) and $`\alpha `$ CMi (Allende Prieto et al. 2002).

## 3. Analysis

### 3.1. Stellar parameters

The effective temperature ($`T_{\mathrm{eff}}`$) has been estimated for all the program stars by one or more authors using the Infrared Flux Method (Blackwell & Lynas-Gray 1994; Alonso et al. 1996, 1999). The $`T_{\mathrm{eff}}`$ of $`\alpha `$ Boo (Arcturus) has been independently derived by di Benedetto (1998) and Griffin & Lynas-Gray (1999). We have adopted the values in Table 1. All our stars were analyzed by Allende Prieto & Lambert (1999), and we accepted the surface gravities derived there from the (B V) colors and the Hipparcos trigonometric parallaxes. The iron abundances in Table 1 are a straight average of previous spectroscopic analyses compiled by Cayrel de Strobel et al. (2001). A recent discussion of the preferred stellar parameters of $`\alpha `$ CMi can be found in Allende Prieto et al. (2002). Fe i lines at 4602.00 and 4602.94 Å with accurate $`gf`$-values and damping parameters were analyzed with synthetic spectra computed from flux-constant LTE model atmospheres to estimate the projected rotational velocities in Table 1.

### 3.2. Measuring the line centers

A variety of procedures have been used previously to determine the central wavelength of a spectral line. The most popular methods involve some type of least-squares polynomial fitting or polynomial interpolation (see, e.g., Neckel & Labs 1990; Allende Prieto & García López 1998). As the stellar spectral lines are asymmetric, the details of the measuring procedure affect the derived central wavelength. Polynomial least-squares fits are well-suited for error estimation, however, there is some arbitrariness in the selection of the ideal order and wavelength interval to use in the fits. For the slowly rotating stars analyzed here, we have found that a reasonable choice, the same for the whole sample, provides stable results regardless of the exact selection for the order and wavelength interval. Our particular choice is to fit 35 mÅ around the line center with a third order polynomial. Preceding the fitting by a cubic-spline interpolation with a step of 5 mÅ proved useful for the HRS spectra, but degraded slightly the quality of the results for the superior solar atlas.

To estimate the error in a derived central wavelength, we determine the horizontal scatter between the polynomial fit F($`\lambda `$) and the observed spectrum $`f(\lambda )`$. For each observed wavelength $`\lambda _i`$ considered in the fit, we find numerically the roots of the polynomial $`F(\lambda )f(\lambda _i)`$ using Laguerre’s method, select the appropriate root, and then derive the rms scatter between the observed and fit wavelengths. Assuming the errors are Gaussian, this scatter is divided by the square root of the number of points entering the polynomial fit to derive the uncertainty $`\sigma _\mathrm{s}`$. As an example, application to the solar flux spectrum provides a nearly normal distribution of uncertainties that peaks at 0.00025 Å and shows bumps on the long-value tail, signaling blends. Nave et al. (1994) classify their wavelengths into four categories, depending on their uncertainty ($`\sigma _\mathrm{l}`$), which ranges from 0.4 mÅ to more than 10 mÅ. Most lines have errors below 5 mÅ and about half of the lines have errors below 1 mÅ111This comparison indicates that an extension and improvement of the laboratory wavelengths for neutral iron and other species is desirable. A spectrograph used to obtain stellar spectra may be well suited for this purpose (see Allende Prieto 2001).

The uncertainties in the measured wavelength shifts, $`\sigma =\sqrt{\sigma _\mathrm{s}^2+\sigma _\mathrm{l}^2}`$, are then converted to velocity, and their distribution $`n(\sigma )`$ is determined. Then, we use $`n(\sigma )`$ to predict the distribution of errors ($`\delta `$) in the observed velocity shifts as

$$N(\delta )=\frac{1}{C}_0^{\mathrm{}}\frac{n(\sigma )}{\sigma }\mathrm{exp}\left(\frac{\delta ^2}{2\sigma ^2}\right)𝑑\sigma ,$$
(1)
where $`C`$ is chosen such that $`_{\mathrm{}}^{\mathrm{}}N(\delta )𝑑\delta =1`$. The derived $`N(\delta )`$ can be compared to the observations in order to constrain the intrinsic scatter. In all cases, the derived and measured histograms are similar.

### 3.3. Searching for correlations

The number of iron lines identified and measured for each star varies from 1340 to 1569222With the exception of $`\alpha `$ CMi, for which only 588 lines were measured. Solar analyses have shown clear correlations of the line shifts with line strength (Allende Prieto & García López 1998) and wavelength (Hamilton & Lester 1999). We explore here the variation of the shifts with line strength, as this is the dominant effect. The line equivalent width and the central line depth have been the most widely used line strength indicators. As pointed out by Pierce & Lopresto (2000), the use of the equivalent width avoids saturation for the strongest lines. On the other hand, equivalent width measurements are distorted by blends in the line wings which may hardly affect line depths. Blends are a particularly serious problem when dealing with cool stars, as some of those in our sample. Therefore, we have adopted the line depth.

The left-hand panels in Fig. 3.2 show the line velocity shifts for the stars in our sample as a function of the residual flux at the center of the line $`f_\mathrm{c}`$ ($`=1`$ line depth). In these panels, the shifts of the individual lines are marked with dots. Grouping the lines in bins of 0.1 in depth, and after applying a median filter, we derive the average values shown as asterisks with error bars. The solid lines are 3th-order polynomial fits to guide the eye, and the broken lines trace the original median values. We have not attempted to disentangle the stellar radial velocities from the gravitational shifts and the atmospheric motions. The absolute values on the vertical axis have been shifted to have a null median velocity shift for the strongest lines; a condition that the polynomial fits were as well forced to satisfy. The right-hand panels show the scatter of the line shifts about the average points. The scatter expected from Eq. (1), mainly the result of the uncertainties in the laboratory wavelengths, the finite width of the lines, and the presence of photometric noise, is also displayed (thick symmetric curve). The $`\sigma (\delta )`$ given in the right-hand panels corresponds to a Gaussian fit by least-squares to the observed scatter. The velocity span of the polynomial fits, $`S`$(v), is listed in Table 1.

A spectroscopic atlas of $`\alpha `$ Boo at a superior resolving power and S/N was published by Hinkle et al. (2000). The corresponding left-hand panel in Figure 3.2 compares the average shifts from our measurements from the atlas (filled circles) with those from the HRS spectrum (asterisks with error bars). The agreement is quite satisfactory. The $`\sigma `$ of the scatter in our spectrum and in the atlas agree within 10 m s-1.

## 4. Discussion

The pattern of the velocity shifts is similar for all the stars, and in qualitative agreement with solar results. The net convective blue shifts of the lines strengthen toward deeper photospheric layers, and therefore affect the weak lines more than the strong ones. Despite our limited sample, there is an indication of a smooth dependence of the average velocity shifts with spectral type. This effect is clearer in our analysis of line shifts than in previous studies using bisectors. As argued in the introduction, this is likely the result of line shifts measurements being less affected by blending than line bisectors in late-G and K type stars.

Our preliminary results indicate that the trend of line shifts as a function of line strength can be determined as a function of spectral type and gravity. At least for late-type dwarfs, assuming the strongest lines in the spectrum are free from convective shifts (as is the case for the solar photosphere), it is possible to correct for convective wavelength shifts to within $`0.2`$ km s-1.

The measured range of wavelength shifts for neutral iron lines is larger in $`\alpha `$ CMi than in cooler stars of classes IV and V. The range of wavelength shifts is also larger, reaching up to $``$ 1 km s-1, for the K giants than for the G-K dwarfs and subgiants in our sample, which is not surprising. Although their lower effective temperatures provide less flux to transport, their lower atmospheric pressure implies larger convective velocities. Fig. 1 suggests that, once we account for the observational errors, there is little room for intrinsic line-to-line scatter in the F-G-K stars sampled here. Convection is supposed to cease in main sequence stars with spectral types earlier than about F2, but other velocity fields must be present in those atmospheres, judging from the asymmetry of the spectral lines (e.g. Gray & Nagel 1989). Determining accurate absolute radial velocities demands an understanding of the wavelength shifts in the spectra of these stars too.

Solar observations have shown differences in the center-to-limb variation of the granulation along the central meridian and the equator (e.g. Beckers & Taylor 1980; Rodríguez Hidalgo, Collados, & Vázquez 1992). Observations of a large number of stars may then show peculiarities at a given spectral type depending on the orientation of the rotational axis. Stars with enhanced magnetic fields are also expected to be peculiar in terms of observed wavelength shifts, as the enhanced magnetic fields may hinder the convective motions. Furthermore, line shifts may vary with time for stars exhibiting an analog of the solar 11 year cycle.

Observations of line shifts for a significant number of stars with high quality are required for a deeper understanding of granulation in stellar atmospheres, its relationship with other stellar phenomena, and the role of surface convection in the structure and evolution of stars. Our first results show that measurements of convective line shifts are a must in order to derive accurate absolute radial velocities.

We thank the HET staff for their outstanding job making possible science observations with HRS since day one. The Hobby-Eberly Telescope is operated by McDonald Observatory on behalf of The University of Texas at Austin, the Pennsylvania State University, Stanford University, Ludwig-Maximilians-Universität München, and Georg-August-Universität Göttingen. NSO/Kitt Peak FTS data used here were produced by NSF/NOAO. This research was supported in part by the NSF (grant AST-0086321).


text:
# Acknowledgements

## Acknowledgements

I would like to thank E. Witten for helpful comments. Financial support by by the European Commission RTN programme HPRN-CT-2000-00131, in which I am associated with INFN, Sezione di Frascati, is gratefully acknowledged.