text:
# Theory of Spin Orientation of Semiconductor Carriers at a Ferromagnetic Interface

## I Introduction

Semiconductor spintronics spinreview  is a new field in which the manipulation of carrier spins is a central issue. The first device proposal in this field was by Datta and Das, datta  which consisted of (1) spin injection from a ferromagnet source contact into a semiconductor two-dimensional electron gas (2DEG), (2) manipulation of the spin with a gate bias through the Rashba spin-orbit effect, rashba  and (3) measurement of the spin with a ferromagnetic drain contact. Progress has recently been made both in spin injection fiederling ; ohno ; zhu ; hanbicki ; johnston  and in control via the Rashba effect nitta  but it appears that there are difficulties in practical implementation of this kind of device. The demand of sufficiently strong spin-orbit interaction requires a narrow gap semiconductor. More basically, the spin-orbit interaction driven by the electric field normal to the plane of the electron gas polarizes the electron spin in the plane normal to its wave vector. To reduce spin-cancellation effects among different wave-vectors, comment  a one-dimensional channel was suggested. datta  In this paper, we describe in theory a different approach to generate and control spin polarization in semiconductors which is based on the spin-dependent properties of thin semiconductor layers in close proximity to a ferromagnet.

Time-resolved Faraday rotation experiments over the last decade have given much insight into electron spin dynamics in semiconductors. faraday  In particular, it has been shown that in lightly doped semiconductors the electron spin coherence can persist for hundreds of nanoseconds, kikkawa  a packet of spin-coherent electrons can be dragged microns using an electrical bias drag  and can maintain their coherence through an interface between two semiconductors, malajovich  and electron spin coherence can produce large nuclear effectsnmr  in the semiconductor through the hyperfine interaction orientation . Recent pump-probe experiments have been performed on $`n`$-doped semiconductor epilayers in contact with a ferromagnet, demonstrating that the proximity of the ferromagnet can induce large nuclear fields fmimprinting  and spontaneous electron spin polarization fmcoherence  in the semiconductor. In a recent letter, ciuti1  we gave a theory of spin-dependent reflection at the semiconductor-ferromagnet Schottky barrier as an explanation of the origin of the spin polarization. In this long paper, we provide a detailed treatment of the spin reflection. We also include calculated results of applying electric bias between the semiconductor and the ferromagnet, stimulated by an ongoing experiment under bias.epstein  We also investigate the different spin reflection properties resulting from optical pumping and from electrical excitation of non-equilibrium carriers in the semiconductor.

The paper is organized as follows. Section II gives a description of the physics of semiconductor electron spin polarization upon reflection at a semiconductor-ferromagnet interface. In Section III we present the general scattering theory. In Section IV, we apply the theory to the case of Schottky junction, showing the role of the Schottky barrier, semiconductor doping and applied bias. Moreover, we compare the case of optical excitation to that of electrical excitation. Conclusions are drawn in Section V.

## II Spin Reflection Off a Ferromagnet

The generation of spin polarization by scattering of spin particles against a spin-polarized target has a long history.mott  We consider the situation where a non-equilibrium distribution of carriers is injected in the semiconductor (either electrically or optically) and determine the transient spin-dependent reflection dynamics. The semiconductor is assumed to be $`n`$-doped to ensure a reasonably thin Schottky barrier (for sizeable quantum-mechanical coupling between the semiconductor and the ferromagnet) and reduced spin relaxation. Thus, an initially spin-compensated group of excited electrons will be reflected by the ferromagnetic interface with a net spin polarization. Since the momentum relaxation of the non-equilibrium carriers is much faster than their spin relaxation, the Fermi sea in the semiconductor conduction band is left with a spin polarization. In this way, the spin reflection produces a ferromagnetic “imprinting” of the semiconductor electrons. More generally, if the non-equilibrium electrons have a pre-existing spin polarization, the ferromagnetic imprinting manifests itself as a tilting of the original polarization vector. Thus, the spin reflection not only generates a spin polarization, but can also rotate a pre-existing ensemble spin.

In the case of optical excitation one should also consider the spin dynamics of the excited holes in the valence band. However, we will neglect any spin effects from the valence band holes. The holes can gain polarization from the reflection at the interface exactly as the conduction electrons, but the corresponding spin polarization is known to decay very fast due to valence band-mixing.uenoyama  In particular the hole spin relaxation time is much shorter than the optical electron-hole recombination time. In the case of optical excitation, a moderate electron doping is essential to provide very long spin lifetimeskikkawa , because the Fermi sea acts as a spin reservoir. The holes will recombine with electrons from the Fermi sea and thus remove a strong source of electron relaxation via exchange with the holes.bap

## III Theory of spin-dependent reflection

To capture the essential physics of the problem we work with a simplified effective mass Hamiltonian in both the semiconductor and the ferromagnet.slon  The effective mass approximation, though suitable for semiconductor heterostructures, can not possibly account for all band properties of the ferromagnetic metal.monier  We shall investigate one improvement of the wave function matching between the dissimilar metal and semiconductor media.einevoll  The general results we present should remain valid with more realistic calculations.

The Hamiltonian of the metal/semiconductor junction is given by

$`H=`$ $``$ $`{\displaystyle \frac{\mathrm{}^2}{2}}{\displaystyle \frac{d}{dz}}\left[{\displaystyle \frac{1}{m^{}(z)}}{\displaystyle \frac{d}{dz}}\right]+U(z)`$ (1)
$`+`$ $`{\displaystyle \frac{\mathrm{\Delta }}{2}}𝝈\widehat{\text{M}}\mathrm{\Theta }(z)+{\displaystyle \frac{g^{}}{2}}\mu _\text{B}𝝈𝐁_{\text{tot}}\mathrm{\Theta }(z),`$
wher the $`z`$-axis is along the growth direction. The first term is the kinetic energy, where $`m^{}(z)`$ is the effective mass which is different for each region. The second term represents the spin-independent part of the potential energy. In the semiconductor region ($`z<0`$), the potential $`U(z)`$ produces the band-bending due to the space charge layer associated with the Schottky barrier. $`U(z)`$ can be tailored by proper heterostructure engineering (for example by inserting a delta-doped layer at the interface) or by applying in situ an electrostatic bias. As we will show in detail in Section IV, the profile of $`U(z)`$ can drastically modify the spin-dependent coupling between semiconductor and ferromagnet. The second line of Eq. (1) contains the spin-dependent part of the Hamiltonian, the exchange interaction operator in the ferromagnet ($`z>0`$) and the Zeeman energy in the semiconductor ($`z<0`$). $`\mathrm{\Delta }`$ is the exchange splitting energy between the majority and minority spin bands in the ferromagnet and $`\widehat{\text{M}}`$ is the unit vector along the direction of the ferromagnet magnetization. In the Zeeman term, $`g^{}`$ is the effective electron g-factor ($`g^{}=0.44`$ for GaAs), $`\mu _\text{B}`$ is the Bohr magneton, and the magnetic field $`𝐁_{\text{tot}}`$ is the sum of the external field and the induced nuclear field, $`𝐁_{\text{tot}}=𝐁+𝐁_\text{N}`$. The Zeeman splitting has a negligible effect on the spin polarization because it is typically several orders of magnitude smaller than the exchange splitting in the ferromagnet and the Schottky barrier height. However, a weak magnetic field is useful as a probe of the ferromagnetic imprinting, because it induces a Larmor precession of the reflection-induced spin polarization which can be detected, for example, by time-resolved Faraday rotation.fmimprinting  The origin of the nuclear field $`𝐁_\text{N}`$ is the hyperfine coupling between the electron and the nuclear spins. If the electron polarization has a component along the external magnetic field vector, the electron polarization is known to induce dynamically a nuclear spin polarization through the Overhauser effect.fmimprinting  The dynamically polarized nuclei produce an effective magnetic field that acts back on the electrons

$$𝐁_\text{N}\frac{g^{}}{|g^{}|}\frac{(𝐒𝐁)𝐁}{B^2+B_0^2}.$$
(2)
Depending on the orientation of the electron spin polarization vector $`𝐒`$ relative to $`𝐁`$, the nuclei can align either along or against the external applied field. $`B_0`$ is a phenomenological parameter to account for the fact that at low applied field the nuclei are unable to align with the applied field due to nuclear spin-spin interactions which tend to destroy the dynamic nuclear polarization.

We denote a majority spin in the ferromagnet as $`|+`$ and a minority spin in the ferromagnet as $`|`$, which are eigenstates of the exchange splitting operator, $`𝝈\widehat{\text{M}}|\pm =|\pm `$ (the magnetization is antiparallel to the net electron spin). In the following, we will work explicitly in this basis. In general, the reflection of semiconductor electrons at the ferromagnetic interface is represented by the reflection matrix $`\widehat{r}(𝐤)`$, which, in the ferromagnet spin basis, is

$$\widehat{r}(𝐤)=\left(\begin{array}{cc}r_{+,𝐤}& 0\\ 0& r_{,𝐤}\end{array}\right),$$
(3)
where $`r_{+,𝐤}`$ is the reflection coefficient for a semiconductor electron with its spin aligned with the majority spin band in the ferromagnet and $`r_{,𝐤}`$ likewise with the minority spin band. This can be expressed in the vectorial form

$$\widehat{r}(𝐤)=\frac{1}{2}\left[(r_{,𝐤}+r_{+,𝐤})𝟙+(𝕣_{,𝐤}𝕣_{+,𝐤})\widehat{𝐌}𝝈\right],$$
(4)
where $`𝟙`$ is the unit matrix.

Suppose that a short excitation pulse injects non-equilibrium electrons into the semiconductor. This perturbation can be applied either by a short optical pump pulse or electrically through lateral contacts or an STM tip (see Fig. 1). Initially, the non-equilibrium electron spin density matrix has the form

$$\widehat{\rho }^\mathrm{i}(𝐤,t=0)=\frac{1}{2}f^\mathrm{i}(k)(𝟙+𝐏^\mathrm{i}𝝈),$$
(5)
with $`f^\mathrm{i}(k)`$ the distribution of injected electrons (for optical excitation it depends on the laser spectrum). The initial polarization is determined by the perturbation. Optically, the polarization is determined by the usual selection rules, namely $`𝐏^\mathrm{i}=0`$ for linearly polarized light and $`𝐏^\mathrm{i}0`$ for elliptically polarized light. In the case of electrical excitation, an initial polarization can be created through injection from magnetic contacts through ordinary spin injection. Upon striking the interface, the non-equilibrium electrons will reflect and transmit, and, to fully account for the polarization in the semiconductor, we must calculate the effect of this reflection on the spin density matrix:

$`\widehat{\rho }^\mathrm{r}(𝐤,t)`$ $`=`$ $`\widehat{r}(𝐤)\widehat{\rho }^\mathrm{i}(𝐤,t)\widehat{r}^{}(𝐤)`$ (6)
$`=`$ $`f^\mathrm{i}(k,t){\displaystyle \frac{1}{2}}[R_0(𝐤)𝟙+𝐑(𝐤)𝝈],`$
where, with the $`𝐤`$-dependence understood,

$`R_0`$ $`=`$ $`{\displaystyle \frac{1}{2}}[(|r_{}|^2+|r_+|^2)+(|r_{}|^2|r_+|^2)\widehat{𝐌}𝐏^\mathrm{i}],`$ (7)
$`𝐑`$ $`=`$ $`{\displaystyle \frac{1}{2}}[(|r_{}|^2|r_+|^2)+(|r_{}|^2+|r_+|^2)\widehat{𝐌}𝐏^\mathrm{i}]\widehat{𝐌}`$
$`+`$ $`\mathrm{Re}(r_{}r_+^{})(\widehat{𝐌}\times 𝐏^\mathrm{i})\times \widehat{𝐌}\mathrm{Im}(r_{}r_+^{})\widehat{𝐌}\times 𝐏^\mathrm{i}.`$
In general the polarization after reflection $`𝐑(𝐤)`$ is different from the original polarization $`𝐏^\mathrm{i}`$ of the excited electrons in the semiconductor.

Since the electron distribution is not in equilibrium, the relaxation of the spin density matrix will be dominated by the relaxation of the hot carrier distribution,$`f^\mathrm{i}(k,t)=f^\mathrm{i}(k)exp(t/\tau _k)`$. This relaxation is spin-independent because it occurs on a much faster time scale than the spin-relaxation time. Thus, the reflection-induced spin-polarization will leave a spin excitation in the semiconductor electron sea. In order to quantitatively determine the effect, we need to calculate the current flow into the ferromagnet during the non-equilibrium transient,

$$\widehat{j}(t)=\widehat{j}^\mathrm{i}(t)+\widehat{j}^\mathrm{r}(t)=_{k_\mathrm{z}>0}\frac{d^3𝐤}{(2\pi )^3}[\widehat{\rho }^\mathrm{i}(𝐤,t)\widehat{\rho }^\mathrm{r}(𝐤,t)]v_z,$$
(9)
where $`v_z=\mathrm{}k_z/m_{\text{sc}}^{}>0`$ is the velocity component normal to the interface. This current flow will be spin-dependent, so that the net spin in the semiconductor from the reflection will be the negative of the spin transmitted into the ferromagnet,

$`𝐒^\mathrm{r}`$ $`=`$ $`\mathrm{Tr}\left\{{\displaystyle \frac{\mathrm{}}{2}}𝝈{\displaystyle 𝑑t[\widehat{j}^\mathrm{i}(t)+\widehat{j}^\mathrm{r}(t)]}\right\}`$ (10)
$`=`$ $`{\displaystyle \frac{\mathrm{}}{2}}{\displaystyle _{k_\mathrm{z}>0}}{\displaystyle \frac{d^3𝐤}{(2\pi )^3}}f^\mathrm{i}(k)\left[𝐑(𝐤)𝐏^\mathrm{i}\right]\tau _kv_z.`$
This is the spin density per unit area, so that the total spin density per unit area in the semiconductor after reflection is

$$𝐒=n^\mathrm{i}\frac{\mathrm{}}{2}𝐏^\mathrm{i}L+𝐒^\mathrm{r},$$
(11)
where $`n^\mathrm{i}=\frac{d^3𝐤}{(2\pi )^3}f^\mathrm{i}(k)`$ is the volume density of pumped electrons and $`L`$ is the semiconductor length perpendicular to the interface. As shown in Eq. (10), the contribution to the imprinted spin from each wavevector channel is proportional to the mean-free path $`\tau _kv_z`$. Only the fraction of non-equilibrium electrons within a mean-free path of the interface participate in the spin-dependent reflection. This implies that for increasing sample length $`L`$, the imprinted spin density per unit volume decreases as $`1/L`$. In the opposite limit, with the semiconductor length $`L`$ shorter than the mean-free path, the situation is different because multiple reflections occur and the electrons are quantum confined. This was the impetus behind a recent proposal for a spin valve with ferromagnetic gates ciuti2  and will be addressed in a future publication.

After the non-equilibrium transient, the imprinted spin of the semiconductor electron sea will decay with the long spin relaxation time. A weak magnetic field induces a Larmor precession which is useful for measuring the imprinted spin. The evolution of the spin in the semiconductor is governed by the Bloch equation

$$\frac{d𝐒}{dt}=\frac{g^{}\mu _\mathrm{B}}{\mathrm{}}𝐁_{\text{tot}}\times 𝐒(t)\frac{𝐒(t)}{T_2^{}}.$$
(12)
The component of $`𝐒`$ orthogonal to the applied field $`𝐁`$ can be extracted from the amplitude of the Larmor precession. In addition, by measuring the effective Larmor frequency it is possible to extract the nuclear field $`𝐁_\text{N}`$ (which is proportional to the component of $`𝐒`$ along the magnetic field). There are three basic vectors in the problem that determine the Faraday rotation: the external field $`𝐁`$, the ferromagnet magnetization $`𝐌`$, and the injected pump polarization $`𝐏^\mathrm{i}`$. Our theory can be tested by systematically changing the relative orientations of these vectors. We now discuss separately the cases of unpolarized and polarized excitation.

### III.1 Unpolarized Excitation

The excited electron population is unpolarized, $`𝐏^\mathrm{i}=0`$, when the non-equilibrium electrons are injected by a linearly polarized laser or a non-magnetic electrical contact. Although initially unpolarized, the non-equilibrium electrons will be polarized by the reflection process. This can be understood in simple terms. Since the majority and minority spin electrons have different wavevectors in the ferromagnet (i.e., the two spin bands are exchange split), the reflection coefficients will in general be different for the two spin channels. This will leave a net spontaneous polarization in the semiconductor.

The polarization process by spin-reflection for unpolarized excitation from a pump beam is shown in Fig. 2 using a simplified model of a parabolic conduction band and a single parabolic valence band. (a) The semiconductor is lightly $`n`$-doped to maximize the spin lifetime, so there is an unpolarized background Fermi sea of electrons. (b) A linearly polarized pump beam excites spin-compensated non-equilibrium electrons and holes. The electrons reflect off the interface. The different reflection coefficients of the two spin channels create a net spin polarization in the semiconductor. The holes undergo a similar process, but the spin-orbit coupling rapidly relaxes any polarization in the valence band. (c) Energy and momentum relaxation then drives the electrons and holes to the lowest available states, with a net spin polarization remaining in the conduction band. (d) After the electron-hole recombination time, all excess holes are gone, and a net polarization is left in the Fermi sea of the conduction band. (e) After the long spin-relaxation time, the electron Fermi sea relaxes back to its original unpolarized state.

The spin in the semiconductor after reflection will be

$$𝐒^\mathrm{r}=\frac{\mathrm{}}{4}\widehat{𝐌}_{k_z>0}\frac{d^3𝐤}{(2\pi )^3}f^\mathrm{i}(k)(|r_{,𝐤}|^2|r_{+,𝐤}|^2)\tau _kv_z,$$
(13)
and hence is determined by the difference in the spin-dependent reflectivities. The spin asymmetry is quantum mechanical in origin, so the sign of the polarization can be positive or negative. The coupling between the ferromagnet and semiconductor is more efficient for either the majority or minority electrons depending on the electronic structure of the ferromagnet and the form of the interface potential in the semiconductor. The spontaneous spin polarization is thus either parallel or antiparallel to the ferromagnet magnetization $`𝐌`$.

Our theory offers a possible explanation of the findings of the recent time-resolved Faraday experiments by Epstein et al..fmcoherence  After pumping non-equilibrium unpolarized electrons in the semiconductor, the experiment found that the system quickly acquired (in tens of picoseconds, consistent with the orbital relaxation time) a spontaneous spin polarization which then decayed with the long spin relaxation time (nanoseconds). The imprinted spin $`𝐒^\mathrm{r}`$ was aligned along $`𝐌`$ by varying the angle between the applied field $`𝐁`$ and the magnetization $`𝐌`$. Interestingly, different ferromagnetic materials gave different signs of the spin polarization relative to $`𝐌`$.

### III.2 Polarized Excitation

When the electrons are pumped with a net polarization $`𝐏^\mathrm{i}0`$, the reflected spin does not align with $`𝐌`$ or $`𝐏^\mathrm{i}`$. There is an extra component which is perpendicular both to $`𝐌`$ and $`𝐏^\mathrm{i}`$, i.e., the reflection off the ferromagnet produces a spin-torque. The geometry of the process is depicted in Fig. 3(a). The amplitude of the extra component is proportional to $`\mathrm{Im}(r_{}r_+^{})`$, implying that the spin torque is due to the different reflection phase-shift for minority and majority spin channels. A pre-existing spin polarization vector can be manipulated using this property of the reflection process, as is shown schematically in Fig. 3(b).

In principle, this kind of effect can be detected through time-resolved Faraday rotation experiments with circularly polarized light. Tilting of the optically injected polarization vector due to reflection would appear as a phase-shift in the Larmor precession. The spin torque term, which is aligned along $`𝐌\times 𝐏^\mathrm{i}`$, would change sign if the magnetization were switched, $`𝐌𝐌`$, resulting in a jump in the phase of the Larmor precession. Experimental studies of this phase shift have been observed. However, a systematic study as a function of the relative geometric orientation of optically-injected polarization, ferromagnet magnetization and magnetic field are necessary to fully test our theoretical predictions.

The spin torque term we describe is analogous to an effect predicted by Slonczewski slon2  and recently demonstrated in all-metallic systems.torque  A spin-polarized current injected into a ferromagnet can tilt the magnetization of the ferromagnet, provided that (1) the ferromagnet magnetization and the polarization of the current are not collinear, and (2) the spin-polarized current injected into the ferromagnet is quite large. However, in the system that we consider, the density of carriers in the semiconductor in so small compared to the density of carriers in the ferromagnet that the tilting of the ferromagnet magnetization is negligible. Instead we find complementary behavior in which the ferromagnet magnetization tilts the polarization of the current.

## IV Results for a Schottky barrier

In this section we investigate in detail the spin reflection at the Schottky barrier between a homogeneously $`n`$-doped semiconductor and a ferromagnetic metal. We define the zero of the potential energy at the Fermi level in the metal. Within the depletion layer approximation, the potential energy in the semiconductor space-charge region ($`z_\mathrm{b}<z<0`$) is given by

$$U(x)=VE_\mathrm{f}^{\mathrm{sc}}+(U_\mathrm{b}+E_\mathrm{f}^{\mathrm{sc}}V)\left(1+\frac{z}{z_\mathrm{b}}\right)^2,$$
(14)
where $`V`$ is the applied bias, $`E_\mathrm{f}^{\mathrm{sc}}`$ is the Fermi kinetic energy in the semiconductor, $`U_\mathrm{b}`$ is the Schottky barrier height, and $`z_\mathrm{b}=\sqrt{ϵ_0U_\mathrm{b}/(2\pi ne^2)}`$ is the depletion width. For high bias $`V>U_\mathrm{b}+E_\mathrm{f}^{\mathrm{sc}}`$, the semiconductor is considered at flat band. The depletion approximation could be replaced with a more realistic model, such as one constructed using the coupled Poisson and Thomas-Fermi equations, but we have explicitly calculated and verified that the deviation from the depletion approximation is negligible.

In Ref. ciuti1,  we gave an analytical approximation for the spin-dependent reflection using an effective rectangular barrier. In this paper, we show the exact numerical results for the realistic Schottky potential in Eq. (14). We also present results for the optical excitation of non-equilibrium electrons, in which pumping takes place in the barrier region, and contrast them with the results for excitation at the Fermi level (electrical injection). These are of interest since both transport and optical experiments have been performed on Schottky barriers in which spin reflection plays a role.

### IV.1 Unbiased Schottky barrier

Under no applied bias, $`V=0`$, the semiconductor and the ferromagnet are in equilibrium, so that the Fermi levels in the two systems are equal. We consider the spin-dependent reflection of electrons which are incident on the Schottky barrier from the semiconductor side with the Fermi kinetic energy $`E_\mathrm{f}^{\mathrm{sc}}=\mathrm{}^2(3\pi ^2n)^{2/3}/(2m_{\mathrm{sc}}^{})`$. The spin-dependence of the reflection coefficient arises from the spin-dependent Fermi velocities in the ferromagnet, which are, in the two-band model, $`v_+^{\text{fm}}=\sqrt{2E_\text{f}^{\text{fm}}/m_{\text{fm}}^{}}`$ and $`v_{}^{\text{fm}}=\sqrt{2(E_\text{f}^{\text{fm}}\mathrm{\Delta })/m_{\text{fm}}^{}}`$ for the majority and minority spin, respectively. $`E_\text{f}^{\text{fm}}`$ and $`m_{\text{fm}}^{}`$ are the Fermi energy and effective mass in the ferromagnet.

The numerical results for the reflection coefficients have been obtained through a finite-difference solution of the Schrödinger equation.liu  The input parameters for the semiconductor are those of bulk GaAs ($`m_{\mathrm{sc}}^{}=0.07m_0`$, $`ϵ=12.9`$). For the ferromagnet, we use the values $`m_{\mathrm{fm}}^{}=1`$, $`E_\mathrm{f}^{\mathrm{fm}}=2.5e\mathrm{V}`$, $`\mathrm{\Delta }=1.9e\mathrm{V}`$). We show the results for the spin reflection difference $`|r_{}|^2|r_+|^2`$ and the spin torque amplitude $`\mathrm{Im}(\mathrm{r}_{}^{}\mathrm{r}_+)`$ as functions of the semiconductor doping $`n`$ and Schottky barrier height $`U_\mathrm{b}`$ are shown in Fig. 4.

Fig. 4(b) shows the spin difference $`|r_{}|^2|r_+|^2`$, which is the spin polarization generated by reflection, directed along the ferromagnet magnetization, $`𝐌`$. Both the Schottky barrier height and the semiconductor doping have significant impact on this generated spin polarization. Variation of either property within experimentally accessible values can cause a change of direction of the spin polarization (indicated by the white-dashed line). The numerical results for the realistic Schottky potential show that our approximated solutions with an effective rectangular barrier in Ref. ciuti1,  are very close. The shape of the contours are identical and only minor discrepancies are present. Note that, even in the absence of a barrier, the velocity mismatch at the interface leads to spin-dependent reflection. At low doping $`n<10^{19}\text{cm}^3`$, the semiconductor velocity is better matched to the minority spin band velocity, but at high doping $`n>10^{19}\text{cm}^3`$, the semiconductor velocity is better matched to the majority spin band velocity (since the Fermi energy for minority spins is smaller than the Fermi energy for majority spins). This is the reason for a change of sign in the spin difference as a function of doping. This argument may be extended to interfaces with small Schottky barriers, but for large Schottky barriers tunneling plays a dominant role. Fig. 4(b) shows the sign change of the spin difference when the barrier height is larger than about $`0.1e\text{V}`$ that joins smoothly with the change of sign from the low-barrier region. For large Schottky barriers, the large wave-vector in the barrier region better matches the majority spin band velocity, inducing the change of sign.

Figure 4(c) shows the spin torque term for polarized incident electrons. This is the reflected polarization component orthogonal to both the ferromagnet magnetization $`𝐌`$ and the electron polarization $`𝐏^\mathrm{i}`$. We find that this term always has the same sign and is analogous to the Kramers-Kronig partner of the spin difference shown in Fig. 4(b).

In the effective mass model used above for both the metal and the semiconductor, we match the wave function and its derivative with the only accommodation of the effect the different crystals being the effective masses($`\psi ^{}/m^{}`$ is continuous). We have investigated the effect of the semiconductor band structure effect in the limit of small gap with the boundary condition einevoll  which modifies the envelope wave function slope on the semiconductor side by a factor of a half,

$$\frac{1}{2m_{\text{sc}}^{}}\frac{\psi _{\text{sc}}^{}}{\psi _{\text{sc}}}=\frac{1}{m_{\text{fm}}^{}}\frac{\psi _{\text{fm}}^{}}{\psi _{\text{fm}}}.$$
(15)
We have found that the spin property dependence on the barrier height and doping density qualitatively unchanged (see Fig. 5).

Our theory implies that the reflected spin polarization is very dependent on the band profile in the semiconductor. The Schottky barrier height $`U_\mathrm{b}`$ and barrier width as measured by the semiconductor doping $`n`$ have a large impact on both the sign and magnitude of the polarization. Both effects can be used to tailor the Schottky barrier to achieve the desired spin polarization through different ferromagnetic materials and different doping concentrations. Our results also show that spin polarization generation through the transmission or reflection of a barrier is a quantum mechanical phenomenon. As is evident from our plot, the reflection can not simply be modelled as depending on the ferromagnetic density of states, or as a spin-dependent resistance that is insensitive to the exact properties of the barrier.

### IV.2 The effect of an applied bias

Now we examine the effect of applying a bias across the Schottky barrier. We assume that applying a positive bias $`V>0`$ raises the semiconductor Fermi level above the ferromagnet Fermi level by an amount $`V`$, i.e., the entire voltage drop is across the barrier; the semiconductor bulk and the ferromagnet bulk are assumed to be flat. The bias may be viewed approximately as just changing the Schottky barrier height, so that for negative bias the Schottky barrier is larger by $`V`$ and for positive bias the Schottky barrier is smaller by $`V`$. For $`V>(U_\mathrm{b}+E_\mathrm{f}^{\mathrm{sc}})`$, the barrier disappears and the entire semiconductor region is a flat band. Moreover, since the tunnelling is a ballistic process, the bias implies that the transmission into the ferromagnet is not at the Fermi level in the ferromagnet. Thus the applied bias changes the majority and minority velocities and hence the matching of the semiconductor electrons (at the Fermi level in the semiconductor $`E_\mathrm{f}^{\mathrm{sc}}`$). In this calculation we choose the Schottky barrier height to be $`U_\mathrm{b}=0.7e\mathrm{V}`$, the value of Fe/GaAs.

In Fig. 6 we plot the spin reflection difference $`|r_{}|^2|r_+|^2`$ for unpolarized excitation. At $`V=0`$, this corresponds to a slice out of Fig. 4 at the Schottky barrier height $`U_\mathrm{b}=0.7e\mathrm{V}`$. As the figure shows, changing the bias makes a change of sign in the spin reflection difference in a similar way to the change of sign as a function of Schottky barrier height in Fig. 4(b). For negative bias $`V<0`$, the Schottky barrier increases in both height and width. The net effect is that the reflection coefficients for both spin channels are very nearly unity, so that the difference between them is negligible. The sign of the polarization is the same as the sign at zero bias, but the spin difference goes asymptotically to zero. The sign of the polarization changes abruptly for positive bias around $`V.5e\mathrm{V}`$, then rises rapidly to its maximum value. This is because the bias makes the Schottky barrier smaller until it disappears completely; the potential becomes a step and the velocity mismatch becomes the important factor. Hence, the high positive bias regime is similar to the low barrier regime in Fig. 4(b). As the bias increases further, the reflection difference decreases, due to the concomitant increase in electron energy in the ferromagnet, so that the exchange splitting $`\mathrm{\Delta }`$ becomes less effective.

These results suggest that a bias could be used to control the sign and magnitude of the spin polarization. For a fixed doping, the reflected spin polarization can be tuned depending on the applied bias. However, it is much easier to get high polarization for one sign of polarization than the other (in Fig. 6(b), the polarization is maximum above the change of sign, but is very nearly zero below the change of sign for low doping).

### IV.3 Special effects of optical excitation

Optical pumping creates special effects in spin reflection at a Schottky barrier not seen under electrical injection. In the latter case, non-equilibrium electrons are injected close to the Fermi level far from the barrier whereas optical excitation creates electrons everywhere in the semiconductor, including in the space charge region. This has a profound effect on the net spin polarization because not all electrons will see the same effective barrier height. Electrons very near the interface will see a much smaller Schottky barrier than electrons pumped far away from the interface. Since the electrons closest to the interface see a smaller and narrower barrier, they will couple to the ferromagnet much more efficiently, so that the polarization is dominated by these electrons. This effect is quantitatively very important because of the exponential dependence of the tunnelling coupling on the barrier height and width. The results of optical excitation are averaged over the different barrier heights.

To fully grasp the difference between optical and electrical excitation, we attempt to quantify the optical effect. First, we assume that the electrons are pumped homogeneously in the space charge region (no contribution from the bulk region is included). In the experiment,fmcoherence  the semiconductor is thin enough for this to be approximately correct. The laser bandwidth is taken to be approximately $`10\mathrm{m}e\mathrm{V}`$ (the typical bandwidth of a pulsed laser in time-resolved pump-probe experiments). We calculate the reflection for a Gaussian distribution of electron energies for each coarse-grained spatial neighborhood in the space charge region. We then average these effects over the different barrier heights due to the homogeneous pumping along the space charge region (see Fig. 7(a)).

Results for the spin refection difference $`|r_{}|^2|r_+|^2`$ from unpolarized excitation are shown in Fig. 7(b) as a function of doping and Schottky barrier height keeping the bias at $`V=0`$. Compare the plot with the case of excitation at the Fermi level, Fig. 4(b). There is not much change for the low barrier region because the electrons pumped high in the barrier do not see a drastically different barrier from the electrons pumped farther away. On the other hand, at high barrier heights, there are drastic changes. The change of sign of the polarization disappears for low doping concentrations. In the case of electrical excitation, at a high barrier all the electrons see a very high and thick barrier which allows for very little tunnelling. In the case of optical pumping, there are some electrons pumped high in the barrier which see a very small barrier and are efficiently coupled to the ferromagnet. The polarization at low barrier is of opposite sign to the polarization at high barrier, but, since the coupling for low barrier is so much more efficient, it dominates and the change of sign disappears.

Figure 7(c) shows the spin reflection difference at a constant Schottky barrier height, $`U_\mathrm{b}=0.7e\mathrm{V}`$, and varying doping and applied bias. In comparison with Fig. 6(b), the optical excitation drastically changes the qualitative aspects of the spin reflection difference. For low doping, the change of sign as a function of bias has been pushed down to negative bias for the same reason as the disappearance of the change of sign in Fig. 7(b) for low doping. The spin reflection difference is dominated by the electrons pumped high in the barrier, which have the opposite sign of polarization compared with the electrons near the Fermi level in the semiconductor.

There is a caveat to these special effects. A more realistic model would use an electron wavepacket to reflect off the interface in order to account for the significant scattering due to efficient optical phonon emission. The electrons pumped high in the barrier would be in a region with almost no background electron density, so that the spin lifetime may be shorter than the lifetime outside the depletion region. Our purpose here is to point out that electrical and optical experiments can give inconsistent results if it is assumed that the underlying mechanism is exactly the same in both cases. Electrically excited electrons stay close to the Fermi level, while optically excited electrons exist all over the barrier. The pumping in the barrier drastically alters the spin reflection difference which would be observed in in experiments on the interface, and hence care must be taken when comparing optically excitation experiments with electrical injection experiments. In particular, the optically excited electrons pumped high in the barrier can have different sign of polarization to electrons lower in the barrier, making it difficult to compare polarizations from experiments involving optical versus electrical injection.

## V Conclusions

We have examined the electronic properties associated with semiconductor electrons that are in proximity with a ferromagnet. Instead of focusing on the injection of ferromagnet electrons into the semiconductor through the interface, we have instead chosen to look at the properties of semiconductor electrons that can interact with a ferromagnetic epilayer. We have focused on the reflection of non-equilibrium electrons from the interface.

For doped semiconductors with ferromagnetic epilayers, we have calculated the spin-dependent reflectivities for electrons incident on the Schottky barrier from the semiconductor. We find that the shape of the barrier has a large impact on the asymmetry between reflection for the two spin channels, and significantly that the sign of the difference can change depending on the system parameters. We have included calculations that mimic the behavior of optical excitation, which can radically affect the spin reflection difference. The ability to control the sign and magnitude of the spin polarization by tuning the properties of the Schottky barrier (through the semiconductor doping or through an applied bias, for example) may be useful for device design.

In contrast to spin injection from a ferromagnet into a semiconductor, spin reflection has the advantage that the processes are kept in the semiconductor. For carriers confined near the interface, multiple reflections from the interface can enhance the single-reflection spin asymmetry, yielding larger polarizations. Such reasoning has led us to propose a spin-valve device with ferromagnetic gates.ciuti2  A detailed study of the coupling of the equilibrium and transport properties of a confined semiconductor electron system in contact with the ferromagnet will be the subject of another long paper.

###### Acknowledgements.

This work is supported by DARPA/ONR N0014-99-1-1096, NSF DMR 0099572, the Swiss National Foundation (for CC), and University of California Campus-Laboratories Cooperation project (for JPM). We thank D.D. Awschalom, R.J. Epstein and R. Kawakami for helpful discussions.


format: text

text:
# Top-Down Fragmentation of a Warm Dark Matter Filament

## 1 Introduction

The Cold Dark Matter crisis on small scales is far from being resolved. First, the highest resolution simulations available still favour ”cuspy” galaxy haloes (Power et al. 2003) with an inner slope for the density profile $`1.2`$, whereas high resolution observations of low surface brightness galaxies are best fit by halo ”cores” (de Block & Bosma 2002; Swaters et al. 2003). Second, the dearth of low luminosity galaxy satellites found $`N`$-body simulations remains, although squelching of dwarfs due to the reionisation of the universe at high redshift may provide a solution (e.g. Benson et al. 2002; Somerville 2002). However, it is unclear that such a process is the definite solution, because the low-mass dark satellite haloes left over still interact gravitationally with the luminous galaxies, exchanging energy and angular momentum with them. Such interactions could have a significant impact on these galaxies’ properties, heating their disks and reducing their sizes. Direct observational probing of dark matter halo substructure through anomalies in the flux ratio of multiple-image quasar lenses could in principle determine whether the scenario previously suggested is the correct one (Metcalf & Madau 2001; Dalal & Kochanek 2002), but there still are major difficulties to overcome in order to obtain reliable estimates of the amount of mass locked in dark substructures (Chen et al. 2003).

An alternate mechanism to solve the satellite problem, which does not suffer from these side effects, is to reduce the power of dark matter fluctuations on small scales. This can be accomplished by changing the nature of the dark matter, (e.g. Kaplinghat, Knox, & Turner 2000; Spergel & Steinhardt 2000; Colin et al. 2000; Bode, Ostriker & Turok 2001; Avila-Reese et al. 2001), by advocating an inflationary model with broken-scale invariance (Kamionkowski & Liddle 2000) or by simply assuming an “unusual” dip in the primordial CDM power spectrum on small scales (Little, Knebe & Islam 2003). The latest constraints from polarisation measurements of the Cosmic Microwave Background Radiation by WMAP (Bennett et al., 2003; Spergel et al., 2003; Kogut et al., 2003) seem to exclude such models, with the major caveat that the current theoretical understanding of the objects that produce the first ionising photons is at best crude. As an example, a reionisation redshift of 11 (on the lowest side of the WMAP published 95 % confidence level) is still reached in a WDM model for which the warmon mass is about 1 keV, provided the objects responsible for photoionising the universe at this redshift emit on average between 100 and 1000 times more ionising photons per baryon than normal stars, as is possible with a top–heavy IMF dominated by accreting black holes (Schneider et al. 2002). Moreover, the constraints become less restrictive as the mass of the warmon increases.

More sophisticated numerical simulations of high-redshift WDM–like universes, combined with a better understanding of the reionising photon emitters are required in order to set definitive constraints on the nature of the dark matter from polarization measurements. This important point set aside, previous simulations (see Paper I; Bode, Ostriker & Turok 2001; Avila-Reese et al. 2001) have indicated that deviation of WDM models from the CDM hierarchical structure formation scenario should be more pronounced in low density regions (e.g. filaments, voids). More specifically, in these simulations, most of the halos close to the resolution limit (i.e. less than 100 particles) seemed to be forming at smaller redshifts than in the standard CDM models, and in a top-down fashion through the break-up of filamentary structures. This phenomenon is of prime interest for galaxy formation and evolution models as it predicts that dwarfs properties strongly depend on the nature of the dark matter. It also gives rise to a series of interesting questions: how much do filaments fragment? When do they fragment the most? What is the mass distribution of the fragments? To tackle this issue we identified one of the filaments in the previous CDM and WDM low resolution runs presented in Paper I and resimulated it with a mass resolution 64 times higher than the original.

This paper proceeds as follows: in Section 2 we briefly give the model parameters and describe the $`N`$-body-simulations which we use in Section 3 to calculate the properties of gravitationally bound halos found in the filament and of the filament itself. Finally, we present our conclusions in Section 4.

## 2 $`N`$-BODY simulations

The $`\mathrm{\Lambda }`$WDM and the $`\mathrm{\Lambda }`$CDM cosmological parameters used in this paper are the same as those of the $`\mathrm{\Lambda }`$WDM2 and $`\mathrm{\Lambda }`$CDM runs presented in Paper I, namely, $`\mathrm{\Omega }_0=1/3`$, $`\lambda _0=2/3`$, $`\sigma _8=0.88`$, $`h=2/3`$, and $`m_{\mathrm{WDM}}=0.5`$keV. We note that they are in good agreement with the latest results from WMAP (Bennett et al. 2003).

As outlined in Paper I, initial conditions for our lower resolution simulations were generated by packing particles from higher resolution initial conditions into groups. For the simulations presented in this paper, we simply identified a filament at $`z=0`$ in our 25 $`h^1\mathrm{Mpc}`$ low resolution box, tracked back in time the Lagrangian volume its particles occupied at $`z=35`$, and unpacked the low-resolution particles into high-resolution particles in that region. In practice, our resimulation of the filament resulted in a 64-fold increase in mass resolution from $`m_p7\times 10^8`$$`h^1\mathrm{M}_{}`$ (1283 particles) in the low resolution runs to $`m_p10^7`$$`h^1\mathrm{M}_{}`$ (5123 particles) in the runs presented here. Each of our $`\mathrm{\Lambda }`$WDM and $`\mathrm{\Lambda }`$CDM filaments contains over 3 million high resolution particles and our force resolution is of the order 3 $`h^1\mathrm{kpc}`$ in the densest regions. The total length of a filament can be represented by the dimensions $`x_b=7`$, $`y_b=4`$, $`z_b=6`$ (in $`h^1\mathrm{Mpc}`$) of the rectangular box enclosing it. We took ”snapshots” of the simulations at the following redshifts: $`z`$=5, 4.5, 4, 3.5, 3, 2.5, 2, 1.75, 1.5, 1.25, 1, 0.5, and 0. These re-simulations were undertaken using the publicly available Adaptive Mesh Refinement code MLAPM (Knebe, Green & Binney, 2001).

Particle groups were identified using the Bound Density Maxima (BDM, Klypin & Holtzman 1997) method, in keeping with Paper I. The BDM code identifies local overdensity peaks by smoothing the density field on a particular scale of the order of the force resolution. We adopted a smoothing of 10$`h^1\mathrm{kpc}`$ for the runs under investigation. These peaks are prospective halo centres. For each of these halo centres we step out in (logarithmically spaced) radial bins until the density reaches $`\rho _{\mathrm{halo}}(r_{\mathrm{vir}})=\mathrm{\Delta }_{\mathrm{vir}}\rho _b`$, where $`\rho _b`$ is the background density and $`\mathrm{\Delta }_{\mathrm{vir}}340`$. This defines the outer radius $`r_{\mathrm{vir}}`$ of the halo. Once we know the outer radius we are able to calculate internal properties of the halos. Note that during this group identification procedure we also check if each individual particle truly belongs to the halo by comparing its velocity to the local value of the escape velocity. We exclude unbound particles from the haloes. This iterative method ensures that we only retain gravitationally bound objects.

## 3 Results

### 3.1 Dark Matter Density Field

We begin our investigation by looking at the large-scale structure density field within and around the filament. Fig. 1 shows the particle distribution for both models with each particle being grey-scaled proportionally to the logarithm of the local (over-)density. A striking feature of Fig. 1 is the marked granularity of the density field in the CDM model at all redshifts in contrast with WDM. This effect reflects the lack of filtering of power on small scales in CDM as opposed to WDM. However, while the WDM filament is obviously different from its CDM counterpart at redshift $`z`$=5, they are very similar at $`z`$=1. This similarity is all the more remarkable that the logarithmic scale used to plot the density fields shown in Fig. 1 artificially enhances the density contrast, thus highlighting differences between the two universes that would have gone unnoticed had the scale been linear.

A further difference of note is the appearance of individual (low-mass) halos during the course of the WDM filament evolution, whereas dense objects only seem to be merging in the CDM filament, therefore decreasing steadily in number.

We have not presented results at redshift $`z=0`$ because the picture at this epoch is qualitatively similar to that at reshift $`z=1`$: most of the interesting processes leading to different looking WDM and CDM filaments occur at high redshifts ($`z>1`$).

### 3.2 Halo Abundance in the Filament

To quantify the discussion in Section 3.1, we compute how the number of halos within the filament evolves with redshift. This number is plotted in Fig. 2, and contains all gravitationally bound particle groups above $`M>10^9`$$`h^1\mathrm{M}_{}`$ (more than 100 particles – thick lines) and $`M>10^{10}`$$`h^1\mathrm{M}_{}`$ (more than 1000 particles – thin lines), respectively. The former simply corresponds to the integral of the mass function over the entire halo mass range available given our resolution.

Fig. 2 demonstrates that independent of redshift, there are more halos in the CDM filament (as already expected from Fig. 1) than in the WDM filament. However, the ratio $`N_{\mathrm{halo}}^{\mathrm{CDM}}/N_{\mathrm{halo}}^{\mathrm{WDM}}`$, which appears to be nearly independent of the mass cut employed, drops significantly from about 50 at $`z=5`$ to approximately 7 at $`z=2`$. This ratio is a factor two smaller for the higher mass-cut of $`10^{10}`$$`h^1\mathrm{M}_{}`$. The depletion of halos in both models from $`z=2`$ to $`z=0`$ is readily explained when we consider that very few objects are created within the filament at low redshifts (see Fig. 5), and therefore, mergers must be the dominant process in determining the number of halos in that redshift range. Indeed, as can be seen in Fig. 1 haloes stream along the web-like structures to eventually become subsumed within larger objects (galaxy groups or clusters) located at the intersections of the filamentary network.

Fig. 3 displays the cumulative mass function of all halos identified within the filament (left panel) along with the mass function of ’new’ halos (right panel). The definition of a ’new’ halo is given in Section 3.4 where Fig. 3 is analysed in more detail, but for the moment one can simply view them as haloes which did not have any progenitor at the previous output. At this point of the analysis, we stress that the cumulative mass function of all halos in the CDM filament closely follows a power-law (with a slope close to -1) whereas in the WDM case one observes a deficit of low-mass objects, in agreement with the results presented in Paper I.

### 3.3 Mass of the Filament

We next addres the question of the fraction of the filament material which gets locked up in galactic halos and how much remains smoothly distributed in the filament. To answer this question, we calculate the ratio of particles belonging to halos with $`M>10^9`$$`h^1\mathrm{M}_{}`$ (identified using the BDM code, cf. Fig. 2) to the total number of particles in the filament. This requires a definition for the boundary of the filament which follows from a previous definition for voids: Einasto et al. (1998) define void particles as having $`\rho _p<\rho _b`$ where $`\rho _p`$ is the local density at the particle position and $`\rho _b`$ the cosmological background density. We adopt the same prescription and call those particles that follow $`\rho _p>\rho _b`$filament particles. The density at the particle position is interpolated from the density field calculated on a regular 2563 grid; as we restrict ourselves to only high-resolution particles we are at all times confined to the filament region (by definition).

The results of this analysis are plotted in Fig. 4. This figure highlights the overabundance of filament halos in CDM with respect to WDM, as the fraction of material in particle groups is at all times higher in the former than in the latter. Only at $`z<1`$ does the ratio $`M_{\mathrm{halos}}/M_{\mathrm{filament}}`$ become comparable in both models, which confirms the similarity of these structures previously noticed at $`z=1`$ (cf. Fig. 1). We note however that from Fig. 1 alone, one might have claimed that most of the material should already be in haloes at $`z=5`$ in the CDM filament, but Fig. 4 argues against it, as only about 10 % of the particles belong to haloes more massive than $`M>10^9`$$`h^1\mathrm{M}_{}`$, for this model, at this redshift. Of course, in the CDM case, this percentage is quite sensitive to the mass cut-off adopted, in the sense that if we had chosen to consider haloes composed of 50 particles or more, it would have increased to about 20 %. Nevertheless, this implies that the gravitational field within the filament is still dominated by the smoothly distributed dark matter component up to redshifts $`z`$ 1, both for the CDM and WDM filaments, and that accretion and merging are indeed important since the total number of halos decreases below $`z=2`$ in both structures whereas the mass of each filament which is locked up within those haloes increases.

### 3.4 New Halos in the Filament

This leads us to the fundamental question of this study, namely how and when does the WDM filament fragment to form haloes? In other words, how important is this ”top-down” fashion of forming haloes compared to the more traditional hierarchical ”bottom-up” scenarios of CDM-like models? To tackle this issue, we constructed the merger tree of each individual halo between consecutive redshift outputs and tagged a halo as ’new’ at a given redshift when no progenitor could be found at the previous output. The requirement for a progenitor to be identified as such was that it had to contain a minimum of 100 particles or 1/10 of the actual halo particles depending on which of these two numbers was the smallest. If such a progenitor was found the halo was marked ’old’. Fig. 5 shows the fraction of ’new’ halos formed per unit time as a function of redshift. As can be seen from this figure, the interesting redshift range for this analysis is from $`z=4.5`$ to $`z=2`$, as one observes a significant halo formation rate within the WDM filament during this time period (cf. Fig. 5).

In order to obtain results not affected by spurious formation or evaporation of haloes between timesteps, we have set a lower mass limit of $`10^9`$$`h^1\mathrm{M}_{}`$ (which translates into 100 particles). This means that if a halo is less massive than this cut-off it has been neglected in our analysis. However, this does not imply that its progenitor(s) had to contain 100 particles. The mass limit for the progenitor(s) is the minimum of 100 and 1/10 of the actual halo particles and hence can be as low as 10 particles (minimal mass for which a bound object can reasonably be detected). We emphasize that in the extreme case where a 10 particle progenitor of a halo is identified, it only affects the ’tag’ of the halo currently under consideration, meaning that it is classified as ’old’ instead of ’new’. More specifically, if such a progenitor turns out to be unphysical (a mere collection of 10 particles which happen to pass by one another quite slowly) it biases the results towards fewer new halos being created in the filament. In this sense, the curves plotted in Fig. 5 are lower limits on the true number of new haloes in the mass range $`10^9`$$`h^1\mathrm{M}_{}`$ and above, as we have made sure that the halos marked as ’new’ were not a part of any structure at the previous redshift. Fig. 5 alone strongly suggests that the bulk of the halo population is forming between $`z=4.5`$ and $`z=2.5`$ within the WDM filament whereas the same figure points at the majority of halos being already in place in the CDM structure by $`z=4`$. To be more explicit, this figure shows that the WDM halo formation rate at $`z=3`$ is still high enough that it could account for all halos present at this redshift, had it been constant in the past. But we also know from Fig. 2 that the maximum number of WDM halos is broadly peaked around $`z=2`$, and from Fig. 5 that the halo formation rate is not constant but steeply increasing with redshift for $`z>2`$. Furthermore, Fig. 5 and Fig. 7 tell us that the halo merger rate is negligible compared to the halo formation rate at $`z2`$. Therefore, we are forced to conclude that most halos were indeed ’born’ in the WDM filament between $`z=5`$ and $`z=2`$. The same analysis, applied to the CDM filament leads to the claim that most halos are in place before $`z=5`$.

This answers the question as to when the halos form out of the WDM filament, but Fig. 6 shows an example of how such a fragmentation occurs. The sequence of panels in that figure traces back in time all the particles found in a typical new halo of mass $`M7.5\times 10^9`$$`h^1\mathrm{M}_{}`$ at redshift $`z=2`$ (i.e. a halo which was identified to be a ’new’ gravitationally bound object at this redshift using the BDM algorithm). It is obvious from Fig. 6 that particles which end up within the virial radius of this halo at $`z=2`$ are part of the filamentary structure until the very moment when the collapse takes place. The quantitative information obtained by tracing the particles pertaining to the final bound halo back in time and provided in Table 1 backs up this qualitative analysis. In particular, we computed in Table 1 the averaged local density (in terms of the background density) at all 749 particle positions, as well as the mean extent of the region occupied by those particles. As a reference we also listed the virial overdensity value $`\mathrm{\Delta }_{\mathrm{vir}}`$ for each of the redshifts under consideration in the last column of Table 1.

To obtain a physical explanation for the fragmentation, we can calculate how the Jeans mass, $`M_J`$, varies with redshift in the WDM filament prior to its fragmentation, assuming its density field to be close to homogeneous on scales larger than the Jeans scale $`\lambda _J`$. These masses and scales are given by the following formulae (e.g. Binney & Tremaine 1987):

$$M_J=\frac{\pi ^{5/2}}{6G^{3/2}}\frac{\sigma _v^3}{\sqrt{(1+\delta )\rho _b}}$$
(1)
$$\lambda _J=\sqrt{\frac{\pi }{G\rho _b(1+\delta )}}\sigma _v$$
(2)
where $`\rho (1+\delta )\rho _b`$ is the density of the filament and $`\sigma _v`$ its averaged velocity dispersion. In the previous equations, $`\delta `$ stands for the (over-)density contrast within the filament with respect to the mean density of the background universe, $`\rho _b`$, at a given redshift. For our cosmological model, Eq. (1) and Eq. (2) become

$$M_J=5.1\times 10^7(1+z)^{3/2}(1+\delta )^{1/2}\left(\frac{\sigma _v}{\mathrm{km}/\mathrm{s}}\right)^3\mathrm{M}_{}$$
(3)
$$\lambda _J=0.13\times (1+z)^{3/2}(1+\delta )^{1/2}\frac{\sigma _v}{\mathrm{km}/\mathrm{s}}\mathrm{Mpc}$$
(4)
A rough estimate of the mass of the fragments can be calculated by assuming that the average density contrast of the filament, $`\delta 10`$ is a slowly varying function of redshift between $`z=5`$ and $`z=2`$, and that typical velocity dispersions are of the order of 20 km/s111This value is based on the averaged velocity dispersion of particles measured in the simulated filaments at local overdensities of around 10.. We then get $`8\times 10^9\mathrm{M}_{}M_J2\times 10^{10}\mathrm{M}_{}`$, to be compared to the mass functions for new haloes given in Fig. 3 (right panel) at redshifts 3.5, 2 and 1. We observe that the masses of the new WDM halo are actually close to these estimated values for the Jeans mass, thus supporting the assumption that they do indeed form through a fragmentation of the filament caused by gravitational instability. Now, the same range of values for $`\delta `$, $`z`$ and $`\sigma _v`$ yields $`53\mathrm{kpc}\lambda _J150\mathrm{kpc}`$ which validates our assumption for the homogeneity of the WDM filament density field on these scales since the free–streaming length $`R_f`$ of a 0.5 keV particle is $`R_f300\mathrm{kpc}`$ corresponding to a free streaming mass $`M_f4\times 10^{10}\mathrm{M}_{}`$. Note that this is not the case for the CDM filament as the free-streaming scale associated with the cold particle is much smaller than the Jeans length, and therefore we cannot consider the density field as homogeneous on these scales. Moreover, one consideration regarding the simulation technique needs to be borne in mind: our $`N`$-body code MLAPM achieves high spatial resolution primarily in high density regions. This means that the force resolution can smooth out the density field on scales above the Jeans mass in lower density regions, with the result that one would then tend to identify haloes as ’new’ that should have already collapsed on smaller scales at earlier times. Nevertheless this possible numerical fragmentation is not too much of a concern in the present study, because even if one was to consider all the new haloes identified in the CDM filament at all redshifts as numerical artifacts, one still would have to invoke physical fragmentation to explain how the more massive new WDM haloes formed (see Fig. 3).

### 3.5 Merger History in Filament

Parallel to fragmentation and creation of new haloes, we also measure the rate of major mergers within the filaments in both models. A major merger event is defined by the ratio of the two most massive progenitors: if this ratio is smaller than 3:1, then the merger is major, otherwise it is minor. As before, a lower mass cut-off of $`10^9`$$`h^1\mathrm{M}_{}`$ is applied for each halo. Fig. 7 indicates that at early times (high redshift) there is no major merger activity in the WDM filament whereas at later times the relative fraction of major mergers exceeds that found in the CDM scenario. This supports our previous claim that haloes at high redshift ($`z>3`$) do not form according to the traditional bottom-up hierarchy (i.e. via mergers) in the WDM filament. Moreover, since the mass distribution of the fragments formed through the top-down mechanism between $`z=5`$ and $`z=2`$ is broadly peaked around the Jeans mass (see Fig. 3), one naturally expects the major merger activity to reach a maximum around $`z=2`$ in WDM because this epoch combines a high number density of low mass halos with a small mass spread in their distribution. In other words, at redshift higher than 2, quite a large fraction of small halos are still forming from fragmentation, and at lower redshift most of the formation and merging have already taken place, broadening the mass distribution function. As a result of these effects, the major merger fraction is strongly suppressed both at low and high redshifts. On the other hand, in the CDM filament, the fraction of major mergers is almost constant in time and reaches more modest values. This is the result of most of the haloes being already in place at high redshift ($`z>5`$) which guarantees that the mass distribution of haloes will broaden smoothly over time as merging and accretion proceed.

## 4 Discussion and Conclusions

We presented the first high resolution $`N`$-body simulation of the top-down formation of dark matter halos through the fragmentation of a dark matter filament. This leads to a number of interesting results: the epoch of formation of small haloes is pushed back in time to lower redshifts ($`2<z<5`$) with respect to the standard CDM model, and the mass distribution of the fragments is broadly peaked around a Jeans mass of a few $`10^9`$ M, corresponding to a gravitational instability of smooth regions with an overdensity contrast around 10 at these redshifts. We need to undertake detailed hydrodynamic and star formation simulations to properly assess the impact of such a mode of formation on dwarf galaxies, but we anticipate the predictions in terms of clustering and the mean age of the stellar populations to be quite different from those of the standard CDM model. We defer this study to a companion paper. Reionisation of filamentary structures might not be a problem since we have shown that massive objects are already in place at z = 5. Of course, the higher the redshift the more difficult it will be to reionise the universe unless the decrease in the number of objects is compensated by an increase in the efficiency of the ionising sources. Our numerical simulations strongly suggest that spurious numerical fragmentation may occur in CDM models due to the softening of the gravity force. However, this is not to be confused with the real fragmentation observed in the WDM simulation presented here because in AMR codes such as MLAPM the force resolution is a function of the mass resolution. Therefore, two haloes with the same number of particles in the $`\mathrm{\Lambda }`$CDM and $`\mathrm{\Lambda }`$WDM simulations are resolved with the same force resolution. This implies that if we were dominated by artificial fragmentation in both $`\mathrm{\Lambda }`$CDM and $`\mathrm{\Lambda }`$WDM we should not observe different upper mass cut-offs in the mass distribution functions of new halos. Looking at Fig. 3 we clearly see that this is not the case: the mass distribution of new halos in $`\mathrm{\Lambda }`$CDM drops off sharply at $`2\times 10^9`$$`h^1\mathrm{M}_{}`$ whereas it extends to masses $`\genfrac{}{}{0pt}{}{_>}{^{}}10^{10}`$$`h^1\mathrm{M}_{}`$ in $`\mathrm{\Lambda }`$WDM. Furthermore, the resolution of our new WDM haloes is high enough that they should have been detected by the BDM algorithm at higher redshift had they been present at that time, since CDM haloes as massive are indeed classified as old ! Finally, we point out that Melott (1982) already showed that galactic haloes were forming via real fragmentation in Hot Dark Matter simulations where the resolution was several orders of magnitudes below what we achieve in this study. In light of these remarks, we conclude that while our results for the formation of massive (at least several hundred particles) WDM haloes through fragmentation are robust, a more detailed study is needed in order to really quantify numerical fragmentation effects on smaller scales, as they might be relevant for studies of dwarf galaxies in cosmological contexts where small haloes are not sufficiently resolved. More specifically, we believe these numerical effects will lead to a later formation of low mass haloes in moderately overdense regions.

## Acknowledgments

The simulations presented in this paper were carried out on the Beowulf cluster at the Centre for Astrophysics & Supercomputing, Swinburne University and the Oxford Supercomputer Centre. AK acknowledges the support of the Swinburne University Research Development Grants Scheme. JEGD acknowledges enlightening discussions with James Taylor about granularity, and with Greg Bryan about numerical resolution issues. The research of JEGD at Oxford is supported by a major grant from the Leverhulme trust. BKG acknowledges the support financial support of the Australian Research Council.


text:
# Avoiding Boundary Effects in Wang-Landau Sampling

## Acknowledgments

This research was partially supported by the Deutsche Forschungsgemeinschaft (DFG) under grant No. Bi314/16-2, the Bundesministerium für Bildung und Forschung (BMBF) under grant No. 03N6015, and by NSF grant DMR-0094422.


text:
# An improved 𝜀 expansion for three-dimensional turbulence: summation of nearest dimensional singularities

## Abstract

An improved $`\epsilon `$ expansion in the $`d`$-dimensional ($`d>2`$) stochastic theory of turbulence is constructed by taking into account pole singularities at $`d2`$ in coefficients of the $`\epsilon `$ expansion of universal quantities. Effectiveness of the method is illustrated by a two-loop calculation of the Kolmogorov constant in three dimensions.

The renormalization-group (RG) method in the theory of turbulence is based on the stochastic Navier-Stokes equation with a Gaussian random force Nelson ; Dominicis79 ; Adzhemyan83 . One of the central problems in this approach is calculation of the Kolmogorov constant $`C`$, i.e. the dimensionless amplitude in the scaling law Monin

$$S_2(r)=C(\overline{}r)^{2/3}$$
(1)
expressing the dependence of the second-order structure function $`S_2(r)`$ on the relative distance $`r`$ in the inertial range $`r_drL`$. Here, $`L`$ is the external length of turbulence, $`r_d`$ the dissipative length and $`\overline{}`$ the energy injection rate per unit mass (which, in the steady state, coincides with the dissipation rate).

Several attempts have been made in the past to solve this problem UFN  \- Jap , but they all suffer from ambiguities in connecting model parameters and observable quantities. As a consequence, there are significant discrepancies in the predicted numerical values for the Kolmogorov constant (the spread is about a factor of two). In this Letter we analyse reasons of this unsatisfactory situation and present results of a calculation based both on an expression of $`C`$ in terms of universal quantities and account of additional singularities arising near and below two dimensions. Rather unexpectedly, the analysis reveals that these singularities have a major effect on the numerical values of observable quantities well above two dimensions. We also show that a partial summation of these singularities is possible and significantly improves the numerical value obtained for $`C`$. To assess properties of the expansion produced within the RG approach, we have carried out the calculation in the two-loop approximation (the results of Refs. UFN  \- Jap  were obtained in the one-loop approximation).

In the RG approach to $`d`$-dimensional turbulence a powerlike correlation function of the random force is often used: $`ffD_0k^{4d2\epsilon }d_f(k)`$. In the RG framework various quantities may be calculated in the form of an $`\epsilon `$ expansion which subsequently must be extrapolated to the physical value $`\epsilon =2`$. For some important quantities the $`\epsilon `$ expansion breaks off, which for the function $`S_2(r)`$ yields the Kolmogorov exponent 2/3 \[as in Eq. (1)\] at $`\epsilon =2`$. To find the Kolmogorov constant the amplitude of this function has to be calculated, which, however, can be done only approximately, because its $`\epsilon `$ expansion does not break off. In calculation of the amplitude, apart from technical difficulties at two-loop order, a principal problem arises as well: the answer for $`S_2(r)`$ has to be expressed in terms of the energy injection rate $`\overline{}`$ \[as in Eq. (1)\] instead of the parameter $`D_0`$ of the powerlike forcing function. This problem has been treated in different ways in Refs. UFN  \- Jap  which has led to different one-loop values of $`C`$.

In Ref. JETP  (see also UFN ; turbo ) the connection between $`D_0`$ and $`\overline{}`$ was sought with the aid of the exact relation

$$\overline{}=\frac{(d1)}{2(2\pi )^d}𝑑𝐤d_f(k).$$
(2)
In the unphysical region $`\epsilon <2`$ this integral has to be cut off at wave numbers of the order of $`\mathrm{\Lambda }r_d^1`$. At fixed $`\overline{}`$ this procedure introduces, first, dependence on $`\epsilon `$ of the form $`D_0(2\epsilon )`$ in $`D_0`$ (which has to taken into account in the construction of the $`\epsilon `$ expansion), and second, an ambiguity connected with the possibility to replace the upper limit $`\mathrm{\Lambda }`$ by $`c\mathrm{\Lambda }`$ with an arbitrary coefficient $`c`$. The first feature is rather natural, because the powerlike forcing $`d_f(k)(2\epsilon )k^{4d2\epsilon }`$ reproduces in the limit $`\epsilon 2`$ the realistic forcing by infinitely large eddies: $`d_f(k)\delta (𝐤)`$. The second feature, however, introduces arbitrariness in the sought connection between $`D_0`$ and $`\overline{}`$ through the coefficient $`c^{2\epsilon 4}`$, which in turn renders the $`\epsilon `$ expansion of $`D_0`$-dependent quantities ambiguous (in Ref. JETP  the simplest choice $`c=1`$ was used ). This is a reflection of the fact that the physical content of the theory remains unaltered when $`D_0`$ is multiplied by an arbitrary function $`F(\epsilon )`$ with $`F(2)=1`$.

Another way to fix the connection between $`D_0`$ and $`\overline{}`$ has been used in Refs. 48  \- Jap . It amounts to the use of an exact relation (for the physical value $`\epsilon =2`$ of the falloff exponent) which allows to connect the spectral energy flux with an integral of a third-order correlation function, the latter being subsequently calculated in the form of an $`\epsilon `$ expansion. The use of this relation in the unphysical region $`\epsilon <2`$ is tantamount to a certain choice of the function $`F(\epsilon )`$ mentioned above.

Thus, the $`\epsilon `$ expansion of the Kolmogorov constant in the model with the powerlike forcing is not unambiguous. Therefore, a better or worse agreement with the experimental value of $`C`$ at one-loop level does not bear much meaning until a procedure for subsequent approximations has been pointed out and the stability of obtained results checked. On the other hand, since the real value of the expansion parameter $`\epsilon =2`$ is not small, it is difficult to expect good quantitative results without estimating – at least approximately – higher orders of the $`\epsilon `$ expansion.

In the model at hand, only quantities independent of $`D_0`$ have rigorously unambiguous dependence on $`\epsilon `$ (we will call them universal). Such quantities are, e.g., critical exponents and dimensionless ratios of structure functions $`S_n(r)`$, the skewness factor $`𝒮S_3/S_2^{3/2}`$ in particular. Calculation of universal quantities with the use of the RG method and the $`\epsilon `$ expansion yields unambiguous results and cannot lead to such ”paradoxes” as different one-loop values for the Kolmogorov constant.

In view of this we have pursued the goal of finding a suitable universal quantity the physical value of which would be simply connected with the Kolmogorov constant, and calculating this quantity with the aid of the RG. The skewness factor $`𝒮`$, connected with the Kolmogorov constant through an exact relation $`C=(4/5𝒮)^{2/3}`$ Monin ; Legacy , might serve as such a quantity. However, in the unphysical region $`\epsilon <2`$ the structure function $`S_2(r)`$ in the model with the powerlike forcing correlation $`d_f(k)k^{4d2\epsilon }`$ contains – at $`\epsilon <3/2`$ – an independent of $`r`$ UV-divergent additive term $`\mathrm{\Lambda }^{24\epsilon /3}`$ \[for $`S_3(r)`$ this problem is absent, see below\]. As a consequence, a straightforward generalization of the skewness factor $`𝒮S_3/S_2^{3/2}`$ to the region $`\epsilon <3/2`$ becomes pointless, because the powers of $`r`$ in this definition do not cancel due to the constant term in $`S_2(r)`$. Therefore, as the desired universal quantity we chose the ”nearest relative” of the skewness factor, the quantity

$$Q(\epsilon )\frac{rS_2(r)/r}{|S_3(r)|^{2/3}}=\frac{rS_2(r)/r}{(S_3(r))^{2/3}}$$
(3)
independent of $`r`$ in the whole region $`0<\epsilon <2`$ which allows to find the physical values of $`𝒮`$ and $`C`$ through the simple relations

$$𝒮=\left[\frac{2}{3Q(2)}\right]^{3/2},C=\frac{3}{2}\left[\frac{12}{d(d+2)}\right]^{2/3}Q(2).$$
(4)
The use of $`S_3(r)`$ in constructing universal quantities is advantageous because it can be found exactly from the spectral energy balance for all $`\epsilon <2`$:

$$S_3(r)=\frac{3(d1)\mathrm{\Gamma }(2\epsilon )(r/2)^{2\epsilon 3}D_0}{(4\pi )^{d/2}\mathrm{\Gamma }(d/2+\epsilon )}.$$

This expression allows, on one hand, to avoid calculation of graphs in construction of the $`\epsilon `$ expansion for $`S_3(r)`$, and, on the other, confirms that passing to the physical limit $`\epsilon 2`$, in which $`\mathrm{\Gamma }(2\epsilon )1/(2\epsilon )`$, requires the dependence $`D_0a(2\epsilon )`$ to arrive at a finite value of $`S_3(r)`$. The choice of the coefficient $`a`$ consistent with (2) yields the ”4/5 law” of Kolmogorov: $`S_3(r)=\frac{4}{5}\overline{}r`$ Monin ; Legacy .

In the usual $`\epsilon `$ expansion at $`d>2`$ the universal quantity $`Q(\epsilon )`$ (3) has the form slovac

$$Q(\epsilon )=\epsilon ^{1/3}\underset{k=0}{\overset{\mathrm{}}{}}Q_k(d)\epsilon ^k.$$
(5)
The RG method allows to find successively the coefficients of $`Q_k(d)`$ as a result of calculation of renormalization constants and scaling functions in perturbation theory (loop expansion). In Refs. UFN  \- Jap  only the one-loop approximation was used in the calculation of the Kolmogorov constant. The results of a two-loop calculation with the aid of the relations (4), (5) have been quoted in Ref. slovac . For the one-loop contribution to $`Q_0(d)`$ in Eq. (5) an analytic expression for all $`d`$ may be obtained:

$$Q_0(d)=(1/3)[4(d+2)]^{1/3}.$$
(6)
The two-loop contribution $`Q_1(d)`$ gives rise to integrals which may be evaluated numerically for any particular values $`d`$. For $`d=3`$ in the calculation of the Kolmogorov constant according to Eq. (4) the values $`C^{(1)}=1.47`$ (one-loop approximation) and $`C^{(2)}=3.02`$ (two-loop approximation) were obtained. Although the two-loop correction is not small, the recommended experimental value of the Kolmogorov constant $`C2.0`$ Monin ; Sreenivasan95  turned out to be in between the values given by the two approximations. Hardly any more could be expected in view of the fact that the value of the expansion parameter is not small. In what follows, we will show that the agreement with the experiment may be significantly improved by an approximate account of the high-order terms of the expansion (5).

Analysis of the dependence of the functions $`Q_k(d)`$ on the space dimension $`d`$ shows that they have singularities at $`d\mathrm{\hspace{0.17em}2}`$. In particular, $`Q_k(d)\mathrm{\Delta }^k`$ for $`2\mathrm{\Delta }d20`$. This means that in the course of $`d`$ tending to 2 the expansion(5) necessarily will become ”spoiled”, because the relative contribution of the high-order terms will grow without limit. In the present two-loop approximation this feature shows in that the ratio $`Q_1(d)/Q_0(d)`$ in the limit $`d\mathrm{}`$ (far away from all singularities) is about 1/20 and monotonically grows with decreasing $`d`$ assuming at $`d=3`$ the value $`1/2`$ of which the major part is brought about by graphs singular in the limit $`2\mathrm{\Delta }=d20`$. This gives rise to hope that summation of leading $`\mathrm{\Delta }`$ singularities in Eq. (5) allows to improve quantitative results of the RG theory.

In the theory of turbulence the space dimension $`d=2`$ is exceptional from both the physical point of view (additional conservation laws, inverse energy cascade) and the formal procedure of UV renormalization, because in the limit $`d2`$ new divergences appear in the graphs of the perturbation theory. These divergences show in the form of poles in $`\mathrm{\Delta }`$ in the coefficients $`Q_n(d)`$ for $`n1`$ in Eq. (5). A consinstent procedure to remove these divergences by an additional renormalization has been developed and gives rise to a two-parameter $`\epsilon ,\mathrm{\Delta }`$ expansion Nalimov . In this Letter it is not possible to dwell on details of calculations. However, we want to point out the following principal issue.

The use of the $`\epsilon ,\mathrm{\Delta }`$ expansion in the theory of turbulence was proposed in Ref. Ronis , whose author points out that an additional renormalization of the random force is required. This was carried out, following Ref. Nelson , by multiplicative renormalization of the random force. There is, however, a major difference between the models of Refs. Nelson  and Ronis . In Ref. Nelson  a local correlation function of the random force is considered ($`k^2`$ – model A, or $`const`$ – model B), wheras in Ref. Ronis  the correlation function is nonlocal $`k^{4d2\epsilon }`$. From general theory of renormalization it is known that counter terms may only be local (see, e.g. Coll ; Zinn ), which means that the renormalization adopted in Ref. Ronis  is inconsistent. A consistent procedure, which we have used, was put forward in Ref. Nalimov . Our two-loop calculation allowed to confirm directly the general statement and show that it is not possible to renormalize the present model by multiplicative renormalization of the random force. We emphasize that, although this fact begins to show in the two-loop approximation, it renders the corresponding one-loop result incorrect as well, since the use of the RG approach is based on the existence of certain relations between all terms of the perturbation expansion which break down in an incorrect renormalization. A detailed discussion of this issue we defer to a separate publication.

In Refs. Nalimov ; Ronis  this two-parameter renormalization procedure (in one-loop approximation) was considered an alternative to the usual $`\epsilon `$ expansion. We exploit it in a different manner – as a way to improve the expansion (5) by carrying out an approximate summation of the high-order contributions.

To single out the leading poles, we express the coefficients $`Q_k(d)`$ in the form

$$Q_k(d)=\mathrm{\Delta }^kq_k(\mathrm{\Delta }),\mathrm{\hspace{0.17em}\hspace{0.17em}2}\mathrm{\Delta }d2,$$
(7)
with a regular function

$$q_k(\mathrm{\Delta })=\underset{l=0}{\overset{\mathrm{}}{}}q_{kl}\mathrm{\Delta }^l.$$
(8)
Substitution of the expressions from Eqs.(7) and (8) in Eq. (5) leads for the quantity $`Q`$ to the representation

$$Q(\epsilon )=\epsilon ^{1/3}\underset{k=0}{\overset{\mathrm{}}{}}\underset{l=0}{\overset{\mathrm{}}{}}(\epsilon /\mathrm{\Delta })^kq_{kl}\mathrm{\Delta }^l.$$
(9)
The $`\epsilon ,\mathrm{\Delta }`$ expansion corresponds to the asymptotic regime $`\epsilon \mathrm{\Delta }0,\mathrm{\Delta }/\epsilon =const`$. Hence, the quantities $`(\epsilon /\mathrm{\Delta })^k`$ in Eq. (9) are not considered small and the powers $`\mathrm{\Delta }^l`$ play the rôle of a formal small parameter. The quantity $`Q`$ from Eq. (9) in the $`n`$th-order approximation $`(n1)`$ is the series

$$\epsilon ^{1/3}\underset{k=0}{\overset{\mathrm{}}{}}\underset{l=0}{\overset{n1}{}}(\epsilon /\mathrm{\Delta })^kq_{kl}\mathrm{\Delta }^lQ_{\epsilon ,\mathrm{\Delta }}^{(n)},$$
(10)
which corresponds to an approximate calculation of the coefficients (7) of the $`\epsilon `$ expansion (5) with the account of $`n`$ terms in the sum (8). For a RG calculation of the quantity $`Q_{\epsilon ,\mathrm{\Delta }}^{(n)}`$ in the $`\epsilon ,\mathrm{\Delta }`$-expansion scheme Nalimov  an $`n`$-loop approximation would be needed.

Let us assume for the moment that we have carried out an $`n`$-loop calculation in the usual $`\epsilon `$ expansion thus determining the following partial sum of the series (5)

$$\epsilon ^{1/3}\underset{k=0}{\overset{n1}{}}Q_k(d)\epsilon ^kQ_\epsilon ^{(n)},$$
(11)
and an $`n`$-loop calculation in the $`\epsilon ,\mathrm{\Delta }`$-expansion scheme as well, hence having determined the quantity $`Q_{\epsilon ,\mathrm{\Delta }}^{(n)}`$ of Eq. (10). Then we have the possibility to amend the sum (11) by an approximate contribution of all higher powers of $`\epsilon ^k`$ not taken into account in Eq. (11). The required information of this contribution is contained in the quantity $`Q_{\epsilon ,\mathrm{\Delta }}^{(n)}`$. To obtain the improved value of $`Q`$ we add the expressions (10)and (11), then subtract once the sum

$$\delta Q^{(n)}\epsilon ^{1/3}\underset{k=0}{\overset{n1}{}}\underset{l=0}{\overset{n1}{}}(\epsilon /\mathrm{\Delta })^kq_{kl}\mathrm{\Delta }^l$$

which enters twice in the sum of Eqs. (10) and (11). Thus, we arrive at the following $`n`$-loop approximation

$$Q_{eff}^{(n)}=Q_\epsilon ^{(n)}+Q_{\epsilon ,\mathrm{\Delta }}^{(n)}\delta Q^{(n)}$$
(12)
for $`Q`$. Our two-loop calculation yields the result

$`Q_{\epsilon ,\mathrm{\Delta }}^{(1)}`$ $`=2\left[{\displaystyle \frac{2(\epsilon +\mathrm{\Delta })^2\epsilon }{3(2\epsilon +3\mathrm{\Delta })^2}}\right]^{1/3},`$
$`{\displaystyle \frac{Q_{\epsilon ,\mathrm{\Delta }}^{(2)}}{Q_{\epsilon ,\mathrm{\Delta }}^{(1)}}}`$ $`=\left[1+\left(0.5181\epsilon +{\displaystyle \frac{1}{6}}\mathrm{\Delta }\right)\right]`$ (13)
for the quantities $`Q_{\epsilon ,\mathrm{\Delta }}^{(1)}`$ , $`Q_{\epsilon ,\mathrm{\Delta }}^{(2)}`$ with the subsequent expressions for $`\delta Q^{(1)}`$ , $`\delta Q^{(2)}`$:

$`\delta Q^{(1)}`$ $`={\displaystyle \frac{2}{3}}(2\epsilon )^{1/3},`$
$`{\displaystyle \frac{\delta Q^{(2)}}{\delta Q^{(1)}}}`$ $`=\left(1+{\displaystyle \frac{2\epsilon }{9\mathrm{\Delta }}}\right)\left[1+\left(0.5181\epsilon +{\displaystyle \frac{1}{6}}\mathrm{\Delta }\right)\right].`$ (14)
Calculating at $`d=3`$ the quantity $`Q_\epsilon ^{(n)}`$ from (11) with the aid of (6) and the value $`Q_2(3)0.4748`$ found in slovac , and substituting the result together with Eqs. (13) and (14) in Eq. (12) we find the quantity $`Q_{eff}`$ in first and second approximation: $`Q_{eff}^{(1)}=1.38`$, $`Q_{eff}^{(2)}=1.84`$. Substitution of these values in Eq. (4) at $`d=3`$ yields for the Kolmogorov constant and the skewness factor the values $`C_{eff}^{(1)}=1.79`$, $`C_{eff}^{(2)}=2.37`$, $`S_{eff}^{(1)}=0.33`$, $`S_{eff}^{(2)}=0.22`$.

In the Table 1 we have quoted for comparison the values of the Kolmogorov constant calculated according to Eq. (4) in the first and second order of the usual $`\epsilon `$ expansion ($`C_\epsilon `$), the double $`\epsilon ,\mathrm{\Delta }`$ expansion ($`C_{\epsilon ,\mathrm{\Delta }}`$), the contribution $`C_\delta `$ in Eq. (4) from the correction $`\delta Q^{(n)}`$ in Eq. (12) and the value $`C_{eff}`$ obtained from Eqs. (4) and (12). In all the cases quoted the recommended experimental value of the Kolmogorov constant $`C2.0`$ lies between the values of the first and second approximation. However, the difference between successive approximations is rather significant both in the $`\epsilon `$ expansion and in the $`\epsilon ,\mathrm{\Delta }`$ expansion, let alone the leading terms of the $`\epsilon `$ expansion of the latter. For the improved $`\epsilon `$ expansion, i.e. for the quantity $`C_{eff}=C_\epsilon +C_{\epsilon ,\mathrm{\Delta }}C_\delta `$ calculated according to Eqs. (12) and (4), however, this difference is about three times smaller leading to a far better agreement with the experimental data.

In conclusion, we have shown that a proper account of the ”nearest singularity” in the coefficients of the $`\epsilon `$ expansion (5) leads to a significant improvement of the results of the two-loop RG calculation at $`d=3`$. We have analysed the effect ot this procedure at other $`d`$ as well. It turned out to reduce significantly the relative contribution of the two-loop correction in the whole range considered $`\mathrm{}>d2.5`$. At the same time this contribution remained large at $`d=2`$, which we think to be an effect of singularities at the next exceptional dimension $`d=1`$.

Obviously, the proposed procedure of approximate summation of the $`\epsilon `$ expansion is applicable not only to the calculation of $`Q(\epsilon )`$, but all universal quantities such as dimensions of composite operators.

The authors thank N.V.Antonov for discussions. The work was supported by the Nordic Grant for Network Cooperation with the Baltic Countries and Northwest Russia No. FIN-6/2002. L.Ts.A., M.V.K. and A.N.V. were also supported by the program “Universities of Russia”. L.Ts.A. and M.V.K. acknowledge the Department of Physical Sciences of the University of Helsinki for kind hospitality.