text:
# Self-similar condensation of rotating magnetized self-gravitating isothermal filaments

## 1 Introduction

Understanding the physics of dense cores is a crucial issue for star formation. In this respect the shape of the cores is of great interest since it depends on the physics that drives their evolution.

Based on a statistical shape analysis, various studies have shown that the observational axis ratio distributions derived from a catalog of dense cores (e.g. Jijina et al. 1999) are not compatible with a distribution of oblate cores and compatible or marginally compatible with a distribution of prolate cores (David & Verschueren 1987, Myers 1991, Ryden 1996, Curry 2002). The same conclusion is reached by Hartmann (2002) who shows that the major axis of the cores in Taurus is preferentially aligned with the filaments in which they are embedded.

Interestingly Jones et al. (2001) and Goodwin et al. (2002) have shown recently that the apparent axis ratio distribution of the dense cores is compatible with the cores being triaxial but being more oblate than prolate. However the physical origin of such structures is presently unclear and the prolate assumption appears thus to be the simplest one which is compatible with the current observations.

It is the purpose of the present study to investigate further the collapse of such elongated structures.

### 1.1 Prolate dense cores

Static prolate cores have been found in various theoretical studies.

Curry (2000) has numerically obtained solutions of the Lane-Emden equations having a prolate shape. These solutions are periodic along the z-axis and reminiscent of the quasi periodicity of the cores sometimes observed in filaments (e.g. Dutrey et al. 1991).

Curry & Stahler (2001) have shown that prolate cores can be solutions of the static MHD equations if the structure of the poloidal magnetic field is such that it compresses the cloud.

Finally cores permeated by an helical magnetic field can also be prolate (Tomisaka 1991, Habe et al. 1991, Fiege & Pudritz 2000c). Fiege & Pudritz (2000c) also show that a large sample of cores permeated by an helical magnetic field have a shape compatible with the observational axis ratio distribution and that the polarization hole effect (Fiege & Pudritz 2000d), observed towards the densest part of the clouds, can be nicely interpreted as a purely geometrical effect due to cancellation of the contribution of the poloidal flux in the core and the toroidal flux in the envelope.

On the other hand, prolate cores can be the result of an instability in magnetized filaments as first investigated by Chandrasekar & Fermi (1953) and recently by Nakamura et al. (1993) who show that a filamentary molecular clouds permeated by an helical field of various intensity can be instable through a combination of the gravitational and the magnetized sausage instability leading to the formation of prolate cores. More recently, Fiege & Pudritz (2000a,b) extend this approach to molecular clouds confined by an external pressure.

### 1.2 Aim of the paper

In this paper we derive self-similar solutions that can describe the collapse of a magnetized, rotating, isothermal self-gravitating filament. We obtain four ordinary equations in $`X`$, the self-similar variable. These equations are very similar to the equations obtained by Larson and Penston for the collapse of a spherical non rotating and non magnetized cloud (Larson 1969, Penston 1969) and extensively studied by Hunter (1977), Shu (1977), Bouquet et al. (1985) and Whitworth & Summers (1985). Our equations admit as exact solutions, some of the explicit solutions which are obtained by Aburihan et al. (2001) and by the author (Hennebelle 2001). These solutions present some important restrictions since they are unbounded and the magnetic poloidal pressure vanishes. The present study demonstrates that some of these solutions are indeed bounded, since a bifurcation arise at the magnetosonic point. Semi-analytical solutions that have a non-vanishing magnetic poloidal pressure are also obtained.

Our goal is to understand better the gravitational contraction of a filament since this process could be relevant in the context of prestellar prolate core collapse, and to obtain the spatial profiles of the fields resulting from the dynamical collapse of a magnetized filament. These profiles could then be compared with observational data and guide a full numerical simulation, although the filamentary geometry used in the present study is a significant departure from a prolate configuration.

The second section of the paper presents the formalism and discusses the physical meaning of the equations. In the third section, we discuss explicit solutions of these equations and we carefully address the question of the magnetosonic critical points. In the fourth section, we solve numerically the ordinary differential equations and discuss the physical meaning of the solutions. The fifth section concludes the paper.

## 2 Self-similar solutions

### 2.1 Reduction to ordinary differential equations

We consider the perfect MHD equations of a self-gravitating gas in cylindrical coordinates. We thus ignore the ambipolar diffusion process.

We make the usual self-similar reduction and set, in the usual notations:

$`X`$ $`=`$ $`{\displaystyle \frac{\varpi }{C_st}},`$
$`\rho (r,t)`$ $`=`$ $`{\displaystyle \frac{1}{4\pi G}}{\displaystyle \frac{R(X)}{t^2}},`$
$`\mathrm{\Phi }(r,t)`$ $`=`$ $`C_s^2\varphi (X),`$
$`V_\varpi (r,t)`$ $`=`$ $`C_sv_\varpi (X),`$ (1)
$`V_\theta (r,t)`$ $`=`$ $`C_sv_\theta (X),`$
$`V_z(r,t)`$ $`=`$ $`\alpha {\displaystyle \frac{z}{t}},`$
$`B_\theta (r,t)`$ $`=`$ $`\sqrt{\mu _0C_s^2/(4\pi G)}{\displaystyle \frac{b_\theta (X)}{t}},`$
$`B_z(r,t)`$ $`=`$ $`\sqrt{\mu _0C_s^2/(4\pi G)}{\displaystyle \frac{b_z(X)}{t}},`$
where $`\alpha `$ is a real number. The radial magnetic component has to vanish in order to ensure zero magnetic divergence.

The problem considered in the present study is thus axisymmetric and independent of z but for the z-component of the velocity field. Also our analysis leaves aside the instability (e.g. gravitational, kink or sausage instability) that could develop and lead to a more complex evolution.

If $`t<0`$, these fields describe a contraction whereas they describe a rarefaction if $`t>0`$.

In the following, we will consider only the case $`t<0`$, i.e. the core contraction until the singularity formation, since during the subsequent evolution, it is expected that a disk will form, as it is the case in the numerical study of Nakamura et al. (1999), that cannot be described by the field stated by Eq. (1).

$`R,\varphi ,v_\varpi ,v_\theta ,b_\theta `$ and $`b_z`$ are dimensionless variables and are expected to be of order unity. The sound speed, $`C_s`$, is about 0.2 km s-1 for a dense core having a temperature of 10 K. Although in principle the self-similar solutions have no characteristic scales, realistic initial and boundary conditions limit the domain of validity. The characteristic spatial scale is given by $`r0.010.1`$ pc, the characteristic time scale by $`t10^{35}`$ years and the density is $`>10^{56}`$ cm-3. The typical magnetic intensity is about $`10\mu `$G for a density of $`10^3`$ cm-3. Then for higher densities, the relation :

$`{\displaystyle \frac{B}{B_0}}=\left({\displaystyle \frac{\rho }{\rho _0}}\right)^k`$ (2)
where $`1/3k1/2`$ is theoretically (Mouschovias 1976, Scott & Black 1980) and observationally inferred (Crutcher 1999).

$`B_z`$ is the poloidal magnetic field. It tends usually to support the cloud against the gravitational collapse. $`B_\theta `$ is the toroidal magnetic component. It has not been taken into account in most of the collapse studies of a magnetized dense core. Such fields can be produced through the stretching of the poloidal magnetic field. The most important effect of the toroidal component is that it usually compresses the gas radially. It can also be a support in the axial direction.

The MHD equations lead with the definitions stated by Eqs. (1) to the following equations.

Continuity:

$`2RXR^{}+{\displaystyle \frac{1}{X}}_X(Xv_\varpi R)+\alpha R=0.`$ (3)
Radial momentum conservation:

$`Xv_\varpi ^{}+v_\varpi v_\varpi ^{}{\displaystyle \frac{v_\theta ^2}{X}}`$ $`=`$
$`{\displaystyle \frac{R^{}}{R}}+\varphi ^{}`$ $``$ $`{\displaystyle \frac{b_\theta }{XR}}_X(Xb_\theta ){\displaystyle \frac{b_zb_z^{}}{R}}.`$ (4)
Azimuthal momentum conservation:

$`Xv_\theta ^{}+v_\varpi v_\theta ^{}+{\displaystyle \frac{v_\varpi v_\theta }{X}}=0.`$ (5)
Axial momentum conservation:

$`\alpha +\alpha ^2=0.`$ (6)
Poisson equation:

$`{\displaystyle \frac{1}{X}}_X(X\varphi ^{})=R.`$ (7)
Azimuthal induction equation:

$`b_\theta Xb_\theta ^{}+_X(v_\varpi b_\theta )+\alpha b_\theta =0.`$ (8)
Axial induction equation:

$`b_zXb_z^{}+{\displaystyle \frac{1}{X}}_X(Xv_\varpi b_z)=0.`$ (9)
Eq. (6) admits the two solutions: $`\alpha =0`$ and $`\alpha =1`$, which correspond respectively to a filament static along the z-axis and to a filament collapsing homologously along this axis. The first configuration has been investigated by Inutsuka & Miyama (1992) in the hydrodynamical case. In the magnetized case, a development around $`X0`$ (see Eqs. 16) implies that $`b_z`$, $`v_\theta `$ and $`b_\theta `$ cannot be all together different from zero which is an important restriction. Moreover, it is unlikely that a cloud will collapse in the radial direction only and not along the z-axis. Thus in the following, we will consider only the second case $`\alpha =1`$.

With $`\alpha =1`$, it is easily seen with Eqs. (3) and (7) that,

$`\varphi ^{}=R(v_\varpi X).`$ (10)
It is also easy to see with Eqs. (3) and (9), that the density, $`R`$, and the axial magnetic fields, $`b_z`$, are proportional: $`b_z=\mathrm{\Gamma }_zR`$, where $`\mathrm{\Gamma }_z`$ is a real number.

With these two last equations, the system of Eqs. (3)-(9), reduces after some algebra to:

$`v_\varpi ^{}`$ $`=`$ $`{\displaystyle \frac{v_\varpi X}{(v_\varpi X)^21b_\theta ^2/R\mathrm{\Gamma }_z^2R}}\times `$ (11)
$`\left({\displaystyle \frac{v_\theta ^2}{X}}+{\displaystyle \frac{1}{X}}+R(v_\varpi X)+{\displaystyle \frac{\mathrm{\Gamma }_z^2R}{X}}{\displaystyle \frac{b_\theta ^2}{XR}}\right),`$
$`R^{}={\displaystyle \frac{v_\varpi ^{}R}{v_\varpi X}}{\displaystyle \frac{R}{X}},`$ (12)
$`v_\theta ^{}={\displaystyle \frac{v_\varpi v_\theta }{X(v_\varpi X)}},`$ (13)
$`b_\theta ^{}={\displaystyle \frac{v_\varpi ^{}b_\theta }{v_\varpi X}}.`$ (14)
### 2.2 Asymptotic behavior

In this section we derive the asymptotic behavior of Eqs. (11)-(14) in the limit $`X\mathrm{}`$. The asymptotic behavior near the origin $`X0`$ will be presented and discussed in some details in Sect. 3.

From Eqs. (11)-(14), it is seen that in the limit $`X\mathrm{}`$, we have:

$`R(X)`$ $``$ $`A_1/X,`$
$`v_\varpi (X)`$ $``$ $`(A_1+(b_\theta ^{\mathrm{}})^2/A_1)\mathrm{ln}(|X|),`$ (15)
$`v_\theta (X)`$ $``$ $`v_\theta ^{\mathrm{}}\mathrm{exp}\left((A_1+(b_\theta ^{\mathrm{}})^2/A_1){\displaystyle ^{|X|}}\mathrm{ln}(u)/u^2𝑑u\right),`$
$`b_\theta (X)`$ $``$ $`b_\theta ^{\mathrm{}}\mathrm{exp}((A_1+(b_\theta ^{\mathrm{}})^2/A_1)/X).`$
where $`A_1`$ is a negative real number.

Consequently at large distance, the azimuthal velocity and the toroidal magnetic component tend to constant values, whereas density (and poloidal magnetic field) decreases as $`1/X`$. The radial velocity tends slowly to infinity. This is a consequence of the infinitely long filamentary geometry. Let us recall that the velocity of the Larson-Penston solutions tends to a finite value at infinity. In case of strong toroidal fields, $`v_\theta `$ and $`b_\theta `$ converge quickly to their asymptotic values.

In the limit $`X\mathrm{}`$, we have:
$`b_z/b_\theta \mathrm{\Gamma }_zA_1b_\theta ^{\mathrm{}}/X0`$ and the dominant forces are the gravitational force and the toroidal magnetic one.

### 2.3 Physical discussions

#### 2.3.1 Physical interpretations

The system of Eqs. (11)-(14) has a clear physical meaning. It is very similar to the equations derived by Larson and Penston (Larson 1969, Penston 1969). As for the Larson-Penston equations, thermal pressure (second term of the right-hand side of Eq. 11) and gravity (third term) are included. The solutions also include rotation (first term), magnetic poloidal forces (fourth term) and magnetic toroidal force (fifth term).

The system of Eqs. (11)-(14) describes the contraction (or expansion) of a filament with an homologous velocity along the z-axis. The gravitational force is zero along the z-axis which is consistent with an infinitely long filament.

In a more realistic situation, the cloud geometry and thus the the potential should evolve with time as for the solutions obtained by Lin et al. (1965) for a cold and unmagnetized cloud.

The velocity along the z-axis does not depend on $`\varpi `$, and consequently the solutions will have no clear physical meaning in the limit $`\varpi \mathrm{}`$. Asymptotic behaviors in this limit, will then not be considered in the following. Because of the homologous velocity field along the z-axis, the validity of the solutions is also restricted to a limited domain not too far from the $`z=0`$ plan, e.g. $`|z|`$ few $`|X_c|`$ ($`X_c`$ is the position of the critical point, see next sections).

The magnetic poloidal field is proportional to density which simply means that ratio between the mass per unit length ($`m=2\pi \rho r𝑑r`$) and the poloidal magnetic flux ($`\psi =B_zr𝑑r`$) is constant into the cloud. It is worth noting that $`m/\psi \mathrm{}`$ when $`t0`$. This is due to the fact that part of the gas flows along the field line.

#### 2.3.2 Magnetosonic singularities

Eqs. (11)-(14) become singular if: $`(v_\varpi X)^21b_\theta ^2/R\mathrm{\Gamma }_z^2R0`$. This is the critical point already obtained in the Larson-Penston equations. The first term is the square of the fluid velocity relatively to the similarity profile, the second one is the sound speed, the two last terms represent the Alfvén velocity, $`C_a`$. The physical meaning is that when the fluid velocity relatively to the similarity profile reaches the fastest velocity at which waves can travel, a singularity occurs. It is only the solutions that pass this singular point that have a physical meaning. The difference from the equations obtained by Larson and Penston is that since the magnetic field is included, the fastest waves are the fast MHD waves and travel with a velocity equal to: $`\sqrt{C_s^2+C_a^2}`$.

In their study, Bouquet et al. (1985) show that for any value of $`\gamma `$, the critical point reached by the homologous solutions is a node rather than a saddle. In the isothermal case, Whitworth & Summers (1995) (see also Bouquet 1984) were able to carry out a general study of the critical point (not restricted to the homologous solution) and show that if $`|X_c|<1`$, the critical point is a saddle and the solution cannot pass it, whereas in the other case, the critical point is a node and the solutions are physical.

## 3 Explicit homologous solutions and magnetosonic critical point

### 3.1 Explicit solutions: the homologous sub-Alfvénic core

In this section, we derive explicit solutions of Eqs. (11), (12), (13), (14). These solutions have already been obtained by the author (Hennebelle 2001) in an extended formulation of self-similarity (Eq. (63) with $`\mathrm{\Gamma }_z=0`$). For completeness and clarity we rederive them here.

We look for homologous solutions and assume the following spatial dependence:

$`R(X)`$ $`=`$ $`R_0,`$
$`v_\varpi (X)`$ $`=`$ $`v_{\varpi ,0}^{}X,`$
$`v_\theta (X)`$ $`=`$ $`v_{\theta ,0}^{}X,`$ (16)
$`b_z(X)`$ $`=`$ $`b_{z,0}=\mathrm{\Gamma }_zR_0,`$
$`b_\theta (X)`$ $`=`$ $`b_{\theta ,0}^{}X=\mathrm{\Gamma }_\varphi R_0X.`$
where, as in Fiege & Pudritz (2000a), $`\mathrm{\Gamma }_z`$ and $`\mathrm{\Gamma }_\varphi `$ are respectively the poloidal and toroidal magnetic flux to mass ratio.

It is straightforward to show that Eqs. (12)-(14) require that:

$`v_{\varpi ,0}^{}`$ $`=`$ $`1/2,`$ (17)
whereas Eq. (11) leads to:

$`1/4+(v_{\theta ,0}^{})^2`$ $`=`$ $`R_0/2+2(b_{\theta ,0}^{})^2/R_0.`$ (18)
In Eq. (11) the terms independent of $`X`$ cancel out and lead to no new constraints, whereas the terms proportional to $`X^2`$ lead to the condition stated by Eq. (18). Consequently, the fields stated by Eqs. (16) with the constraint stated by Eq. (17) constitute the asymptotic form of the solutions of Eqs. (11)-(14) at $`X0`$. The homologous solutions require two constraints (Eqs. 17 and 18) and there are five parameters, which leaves three independent degrees of freedom (e.g. the poloidal and toroidal magnetic flux to mass ratio and the density).

Eq. (18) has a clear physical meaning. It is simply the radial momentum conservation in the pressureless limit. The first term of the left-hand side is the square of $`v_{\varpi ,0}^{}`$, and represents the advection term, the second one is the centrifugal force, the first term of the right-hand side is the gravitational force whereas the second one is the toroidal magnetic pinching. At this point there is no explicit poloidal support since the axial magnetic component is uniform. However, there is an implicit effect due to the critical point (see Sect. 3.2 and 3.3).

If $`b_{\theta ,0}^{}=0`$ and $`v_{\theta ,0}^{}=0`$, then $`R_0=1/2`$ and $`v_{\theta ,0}^{}=1/2`$. This is very similar to the Larson-Penston solution that describes a spherical cloud having uniform density collapsing homologously, for which $`R_0=2/3`$ and $`v_{\varpi ,0}^{}=2/3`$.

If $`R_0<1/2`$, small values of $`b_{\theta ,0}^{}`$ are not allowed, whereas if $`R_0>1/2`$, small values of $`v_{\theta ,0}^{}`$ are forbidden. In the first case, gravity is too weak and cannot produce the homologous collapse by itself, the condensation is induced by the toroidal pinching. In the second case, gravity is too strong to produce the homologous collapse and the cloud must be supported by the centrifugal force.

These solutions are valid until the critical point is reached and describe a sub-Alfvénic (relatively to the similarity profile) dense core contracting homologously.

### 3.2 Critical point

A detailed analysis of the critical point can be achieved for the solutions stated by Eqs. (18). This has already been investigated for the Larson-Penston solutions by Whitworth & Summers (1985) in the isothermal case and by Bouquet et al. (1985) for any values of $`\gamma `$.

For the solutions stated by Eqs. (18), we have:

$`(v_\varpi X)^21b_\theta ^2/R\mathrm{\Gamma }_z^2R`$ $`=`$
$`({\displaystyle \frac{1}{4}}(b_{\theta ,0}^{})^2/R_0)X^2`$ $`1\mathrm{\Gamma }_z^2R_0.`$ (19)
Thus if $`\sqrt{R_0}/2<b_{\theta ,0}^{}<\sqrt{R_0}/2`$, the critical point occurs at:

$`X_c=\sqrt{{\displaystyle \frac{1+\mathrm{\Gamma }_z^2R_0}{1/4(b_{\theta ,0}^{})^2/R_0}}},`$ (20)
whereas if $`(b_{\theta ,0}^{})^2>R_0/4`$, there is no singular point since the velocity relatively to the similarity profile is always smaller than the Alfvén speed.

The topological nature of the critical point can be studied by introducing a new parameter, $`s`$, such that:

$`{\displaystyle \frac{dX}{ds}}`$ $`=`$ $`D,`$
$`{\displaystyle \frac{dR}{ds}}`$ $`=`$ $`RN{\displaystyle \frac{R}{X}}D,`$
$`{\displaystyle \frac{dv_\varpi }{ds}}`$ $`=`$ $`(v_\varpi X)N,`$ (21)
$`{\displaystyle \frac{dv_\theta }{ds}}`$ $`=`$ $`{\displaystyle \frac{v_\theta v_\varpi }{X(v_\varpi X)}}D,`$
$`{\displaystyle \frac{db_\theta }{ds}}`$ $`=`$ $`b_\theta N,`$
where,

$`N`$ $`=`$ $`\left({\displaystyle \frac{v_\theta ^2}{X}}+{\displaystyle \frac{1}{X}}+R(v_\varpi X)+{\displaystyle \frac{\mathrm{\Gamma }_z^2R}{X}}{\displaystyle \frac{b_\theta ^2}{XR}}\right),`$ (22)
$`D`$ $`=`$ $`(v_\varpi X)^21b_\theta ^2/R\mathrm{\Gamma }_z^2R.`$
The linearisation of the solutions stated by Eqs. (18) near the critical point, stated by Eq. (20) leads to a relation $`d𝐘/ds=[M]𝐘`$, where $`𝐘=[X,R,v_\varpi ,v_\theta ,b_\theta ]`$ and where the matrix $`[M]`$ is given by the value of $`(d𝐘_i/ds)/𝐘_j`$ at the critical point. The matrix $`[M]`$ is easily calculated with Eqs. (21) and the derivative of $`N`$ and $`D`$. We have:

$`D/𝐘`$ $`=`$ $`[X_c,{\displaystyle \frac{b_\theta ^2}{R^2}}\mathrm{\Gamma }_z^2,X_c,0,{\displaystyle \frac{2b_\theta ^c}{R_0}}],`$ (23)
$`N/𝐘`$ $`=`$ $`[{\displaystyle \frac{3R_c}{2}},{\displaystyle \frac{X_c}{2}}+{\displaystyle \frac{\mathrm{\Gamma }_z^2}{X_c}}+{\displaystyle \frac{(b_\theta ^c)^2}{X_cR_c^2}},R_c,{\displaystyle \frac{2v_\theta ^c}{X_c}},{\displaystyle \frac{2b_\theta ^c}{X_cR_c}}]`$
where $`R_c`$, $`v_\varpi ^c`$, $`v_\theta ^c`$ and $`b_\theta ^c`$ are the values of the fields at the critical point.

The study of the eigenvalues of $`[M]`$ allows to determine the topological nature of the critical point. The characteristic polynômial $`P`$, of $`[M]`$ is:

$`P(\lambda )=\lambda ^3\left(\lambda ^2X_c\lambda +{\displaystyle \frac{X_c^2}{4}}\left(1{\displaystyle \frac{8(b_{\theta ,0}^{})^4}{R_0^2}}\right)\right)=0.`$ (24)
Its roots are zero (third order) and

$`\lambda _1`$ $`=`$ $`{\displaystyle \frac{1}{2}}X_c\left(1+2\sqrt{2}{\displaystyle \frac{(b_{\theta ,0}^{})^2}{R_0}}\right),`$
$`\lambda _2`$ $`=`$ $`{\displaystyle \frac{1}{2}}X_c\left(12\sqrt{2}{\displaystyle \frac{(b_{\theta ,0}^{})^2}{R_0}}\right).`$ (25)
0 is eigenvalue of $`[M]`$ at the third order since the vectors $`((d𝐘_i/ds)/𝐘)(X_c)`$ are linear combinations of $`(D/𝐘)(X_c)`$ and $`(N/𝐘)(X_c)`$ (see Eqs. 21).

In the neighbourhood of the critical point, one has:

$$𝐘(s)=\mathrm{\Sigma }_{i=1,5}\alpha _i𝐕_i\mathrm{exp}(\lambda _is)$$
(26)
where $`𝐕_i`$ are the eigenvectors and $`\alpha _i`$ are real numbers and where $`\lambda _3,\lambda _4`$ and $`\lambda _5`$ are the three vanishing eigenvalues.

Consequently, since $`\lambda _1`$ and $`\lambda _2`$ are both negative (remembering that $`(b_{\theta ,0}^{})^2/R_0<1/4`$), the critical points are nodes rather than saddles. The solutions are thus able to cross the critical point. It does not necessarily mean that the solutions are always physical, since after this first critical point, it could be possible that the solutions reach a second critical point through which it could be unable to pass. Only a numerical study seems capable to give an answer to this question.

From Eqs. (21), we have the relations:

$`dv_\theta `$ $`=`$ $`{\displaystyle \frac{v_\theta ^cv_\varpi ^c}{(v_\varpi ^cX_c)X_c}}dX,`$
$`db_\theta `$ $`=`$ $`{\displaystyle \frac{b_\theta ^c}{v_\varpi ^cX_c}}dv_\varpi ,`$ (27)
$`dR`$ $`=`$ $`{\displaystyle \frac{R_c}{v_\varpi ^cX_c}}dv_\varpi {\displaystyle \frac{R_c}{X_c}}dX,`$
where $`v_\theta ^c=v_\theta (X_c)`$, $`v_\varpi ^c=v_\varpi (X_c)`$ and $`R_c=R(X_c)`$.

Thus, when the solutions cross the critical point, $`v_\theta `$, $`b_\theta `$ and $`R`$ vary according to these relations.

When the solution crosses the critical point, $`D`$ becomes positive. This leads to the condition:

$`dV_\varpi >{\displaystyle \frac{X_c^2+R_0^2\mathrm{\Gamma }_z^2(b_{\theta ,0}^{})^2X_c^2/R_0}{X_c^2+2R_0^2\mathrm{\Gamma }_z^2(b_{\theta ,0}^{})^2X_c^2/R_0}}dX.`$ (28)
According to Eqs. (27) and (28), and as pointed out by Whitworth & Summers (1985), there is a one dimensional space of possible solutions for $`X<X_c`$ which are compatible with a given solution for $`X>X_c`$.

### 3.3 Expected behaviors

The homologous solutions stated by Eqs. (16) and (18) allow us to anticipate some trends that we will recover numerically in the next section.

From Eq. (20), it is seen that, $`|X_c|`$ increases with $`\mathrm{\Gamma }_zR_0=b_z`$ and $`b_{\theta ,0}^{}`$. From Eqs. (16), it is seen that the value of $`v_\varpi (X_c)`$, $`v_\theta (X_c)`$ and $`b_\theta (X_c)`$ increases with $`|X_c|`$, whereas $`R(X_c)`$ and $`b_z(X_c)`$ do not depend on $`|X_c|`$.

Consequently, the following trends are expected:

If the density increases, $`|X_c|`$, the position of the magnetosonic point, decreases (remembering that $`b_z=\mathrm{\Gamma }_zR`$). This is simply due to the fact that the Alfvén speed decreases and thus the velocity of the fastest waves decreases as well. With Eqs. (18), it is seen that the rotation ($`v_{\theta ,0}^{}`$) increases. This means that when density increases, the gravitational force increases as well and must be counterbalanced by the centrifugal force.

If the poloidal magnetic intensity increases, $`|X_c|`$ increases and the value of the radial velocity and toroidal magnetic field at this point increase as well.

Consequently, the trend is that the larger the magnetic poloidal intensity is, the larger $`|v_\varpi |`$, $`|v_\theta |`$ and $`|b_\theta |`$ will be in the super-Alfvénic part ($`X<X_C`$). This is simply due to the fact that the magnetic poloidal pressure leads to an effective sound speed equal to: $`1+\mathrm{\Gamma }_z^2R`$ (see Eq. 11) and that $`V_\varpi ,V_\theta `$ and $`B_\theta `$ are proportional to the sound speed.

If the toroidal field increases, the position of the critical point increases and the value of the radial velocity as well, which simply means that the cloud is more compressed by the toroidal pinching. As we already said, if $`(b_{\theta ,0}^{})^{1/2}>R_0/2`$ there is no critical point since the velocity of the fastest waves is always greater than the velocity of the fluid relatively to the similarity profile. Rotation must increase (Eq. 18) in order to counterbalance the toroidal pinching.

An important parameter is the value of the radial velocity at large $`X`$. According to Eq. (15), the radial velocity slowly diverges when $`X\mathrm{}`$ but, as we already said, the physical meaning of the solutions is unclear in this limit. On the other hand, Eqs. (17) and (20) give for $`v_\varpi (X_c)`$ the value $`X_c/21`$ (in sound speed unit). Consequently, it is expected that the typical radial velocity at few $`X_c`$ should be as high as few times the value of the sound speed.

The smallest value of $`|v_\varpi (X_c)|`$ is 1 and is reached for the unmagnetized solution. In that case $`X_c=2`$. Let us recall that for the Larson-Penston solution, $`X_c=3`$ and $`v_r(X_c)=2`$. The present unmagnetized solutions are thus less dynamical, in the inner homologous part and in the radial direction, than the Larson-Penston solution.

## 4 Numerical results

### 4.1 Method

In this section we investigate numerically the system of Eqs. (11)-(14). These equations are four ordinary differential equations that can be easily integrated by a fourth order Runge-Kutta method.

At $`X=0`$, the system is singular and one has to develop the solutions to first order and to start the integration at $`X0`$. The asymptotic development towards $`X0`$ is given by Eqs. (16) and $`v_{\varpi ,0}^{}=1/2`$ (Eq. 17). Consequently, the number of free parameters is four ($`R_0`$, $`v_{\theta ,0}^{}`$, $`\mathrm{\Gamma }_z`$ and $`b_{\theta ,0}^{}`$) and there is one more free parameter when the solutions cross the critical point, leading to a total of five free parameters. It is thus difficult to explore all the parameter space and we will restrict the investigation to a limited parameter range. We start at $`X=10^5`$ and use a spatial step $`dX=10^3`$.

The integration through the critical point requires some care. As we already said, the solutions for $`X<X_c`$ compatible with a solution for $`X>X_c`$ are not unique and are fixed by Eqs. (27) and by the conditions stated by Eq. (28). We proceed as follow. When the critical point is reached ($`D`$ reverse sign), we select a value of $`dv_\varpi `$ satisfying Eq. (28) and give to the other fields the values stated by Eqs. (27). In the following we restrict our study to $`dv_\varpi =0`$.

We will first consider homologous cores (thus following Eq. 18) and then relax this assumption and consider the more general case of non-homologous cores.

### 4.2 Homologous cores and non homologous super-Alfvénic envelope

The homologous solutions stated by Eqs. (16) and (18) play an important rôle since except for the geometry they are equivalent to the Larson-Penston solution that has been found to be in reasonable agreement with numerical simulations in previous studies (Larson 1969, Hunter 1977, Blottiau et al. 1988, Foster & Chevalier 1993). Moreover, these solutions are explicit and allow us to check carefully the validity of the numerical integration.

In this section we present the numerical integration of Eqs. (11)-(14) with the initial density $`R_0=1/2`$. This value is important since according to Eq. (18), gravity is exactly equal to the advection term and all values of $`v_{\theta ,0}^{}`$ and $`b_{\theta ,0}^{}`$ are allowed.

The other parameters are: $`v_{\varpi ,0}^{}=1/2`$ and in order to explore a large range of parameters $`b_z=0,0.5`$ and $`1`$ corresponding to $`\mathrm{\Gamma }_z=0,\mathrm{\hspace{0.17em}1}`$ and $`2`$ respectively and $`b_{\theta ,0}^{}`$ equal to $`\alpha \sqrt{R_0/4}`$, whith $`\alpha =0.1,0.5`$ and 0.9. $`v_{\theta ,0}^{}`$ is calculated according to Eqs. (18). We integrate until $`X=5X_c`$ and display the solutions as a function of $`X/X_c`$ (see Eq. 20). Thus, the bifurcation always occurs at $`X/X_c=1`$.

Before presenting the numerical results, we briefly discuss the field intensities.

#### 4.2.1 Field intensity

For orientation, we estimate the typical values of the poloidal field in the subsonic homologous core and the value of the toroidal component at the critical point as a function of $`\alpha `$. With Eq. (1)-(20), one has:

$`B_z=\sqrt{\mu _0C_s^2\rho }\mathrm{\Gamma }_z\sqrt{R_0},`$ (29)
$`B_\theta =\sqrt{\mu _0C_s^2\rho }\sqrt{{\displaystyle \frac{\alpha ^2}{1\alpha ^2}}}\sqrt{1+\mathrm{\Gamma }_z^2R_0}.`$ (30)
If a canonical density of $`10^6`$ cm-3 in the collapsing filament is adopted and for $`R_0=1/2`$, $`\mathrm{\Gamma }_z=1`$, then typically in the weak, intermediate and strong toroidal field cases, toroidal field strengths of $``$15, $``$90 and $``$320 $`\mu G`$ are associated with poloidal field of about $`90\mu G`$.

#### 4.2.2 Weak toroidal magnetic field

Fig. 1 displays the results for $`b_{\theta ,0}^{}=0.1\sqrt{R_0/4}`$.

For $`\mathrm{\Gamma }_z=0`$, the solution is weakly magnetized (since the toroidal field is weak) and except for the geometry it is equivalent to the Larson-Penston solution. As pointed out by Whitworth & Summers (1985) this solution describes a strong compression wave propagating inwards into the cloud, the bifurcation point being the head of this compression wave.

It is seen that for $`X/X_c<1`$ the solution is homologous as predicted by Eqs. (18). At $`X/X_c>1`$ the density decreases, whereas $`|v_\varpi |`$, $`|v_\theta |`$ and $`|b_\theta |`$ increase. The density decreases less rapidly for large values of $`b_z`$ ($`\mathrm{\Gamma }_z`$) and as discussed in Sect. 3.3, the velocity and the toroidal magnetic field increase with $`b_z`$ (see Eq. 16 and 20).

The solutions are weakly peaked with $`R(0)/R(5X_c)2`$. The velocity and toroidal magnetic fields continue to decrease rapidly after the critical point. We have $`v_\varpi (X_c)/v_\varpi (5X_c)0.25`$.

The radial velocity is roughly $`8`$ times larger than the azimuthal one. When $`b_z=0.5`$ and $`1`$ the poloidal magnetic intensity dominates the toroidal magnetic intensity in the core (by a factor $`5`$ at $`X_c`$ for $`\mathrm{\Gamma }_z=1`$, and $`7`$ for $`\mathrm{\Gamma }_z=2`$). At $`X5X_c`$, the two components are almost equal for $`\mathrm{\Gamma }_z=1`$ whereas for $`\mathrm{\Gamma }_z=2`$, $`b_z/b_\theta 1.5`$.

The radial velocity at $`5X_c`$ is around $`45C_s`$ for low and intermediate poloidal magnetic intensity ($`\mathrm{\Gamma }_z=0`$ and 1) and it is much higher ($`7`$) for a strong poloidal magnetic field ($`\mathrm{\Gamma }_z=2`$).

#### 4.2.3 Intermediate toroidal magnetic field

The case $`b_{\theta ,0}^{}=0.5\sqrt{R_0/4}`$ is displayed in Fig. 2. The toroidal magnetic force in the core ($`2b_{\theta ,0}^{}X/R_0`$) is about half of the gravitational force ($`R_0X/2`$). Thus the toroidal field is important but not dominant.

It is seen that the bifurcation at the magnetosonic point is stiffer. The density decreases more rapidly ($`R(0)/R(5X_c)2.5`$) than in the previous case. This is due to the fact that the cloud is compressed by the toroidal magnetic field.

The velocity and the magnetic field evolution is much flatter in the envelope than in the core, we have $`v_\varpi (X_c)/v_\varpi (5X_c)0.3`$. Thus $`v_\varpi (5X_c)`$ is even slightly lower than for the case with weak toroidal magnetic field in spite of the higher value of $`v_\varpi (X_c)`$.

It is also seen that the density variation with the poloidal magnetic intensity is small compared to the variation of the velocity and the toroidal magnetic component ; we discuss this further in the next section.

The amplitude of the radial velocity is about $`2/3`$ the amplitude of the azimuthal one. The poloidal magnetic component dominates the toroidal one in the core (ratio 0.5-0.75 at $`X=X_c`$). The two component become comparable at $`X1.5X_c`$ which is earlier than for the weak toroidal magnetic field. At $`X/X_c=5`$ $`b_z/b_\theta 3`$.

#### 4.2.4 Strong toroidal magnetic field

The case $`b_{\theta ,0}^{}=0.9\sqrt{R_0/4}`$ is displayed in Fig. 3. The toroidal force is larger than the gravitational force (factor $`1.62`$) and the solution describes the inward propagation of a strong magnetic compression wave into the cloud.

The most striking features of the solutions with strong toroidal magnetic field is that the density is almost independent on the poloidal magnetic intensity, whereas the other fields are roughly proportional to $`\mathrm{\Gamma }_z`$. This is due to the fact, that since the toroidal field is strong ($`b_\theta ^2(X_c)/R>3`$), the magnetic and centrifugal forces are large compared to the thermal pressure. Thus in Eq. (11), the thermal pressure can be neglected (term $`1`$ in $`D`$ and $`1/X`$ in $`N`$). In this limit, it is then seen that Eq. (11) is invariant under the rescaling:

$`R`$ $``$ $`R,`$
$`v_\varpi `$ $``$ $`Kv_\varpi ,`$
$`v_\theta `$ $``$ $`Kv_\theta ,`$ (31)
$`b_z`$ $``$ $`Kb_z,`$
$`b_\theta `$ $``$ $`Kb_\theta .`$
Another interesting feature is that the radial velocity and the toroidal magnetic field are almost constant just after the critical point. They then decrease significantly again around $`X/X_c1.5`$. The radial velocity then decreases according to Eqs. (15) whereas $`v_\theta `$ and $`b_\theta `$ decrease more slowly and rapidly converge toward a constant value.

The solutions are even more peaked than in the previous cases, with $`R(0)/R(5X_c)5`$ whereas the radial velocity is flatter with $`v_\varpi (X_c)/v_\varpi (5X_c)0.5`$.

The azimuthal velocity is about $`1.31.5`$ larger than the radial velocity. The toroidal magnetic field rapidly dominates the poloidal one since at $`X/X_c0.5`$, the two field intensities are comparable. At $`X/X_c`$, $`b_\theta /b_z2`$ and at $`X/X_c=5`$, we have $`b_\theta /b_z10`$.

For $`\mathrm{\Gamma }_z=0`$ or $`\mathrm{\Gamma }_z=1`$, $`v_\varpi (5X_c)45C_s`$.

#### 4.2.5 Varying the density: forbidden values

In this section, we vary the density of the homologous core in order to explore the dependence of the fields on this parameter. We look carefully at two particular cases and briefly comment on other values. In order to make comparison with previous cases easier we choose $`\mathrm{\Gamma }_z=1`$ (dotted line of previous section) and $`b_\theta =0.1\sqrt{R_0/4},\mathrm{\hspace{0.17em}0.5}\sqrt{R_0/4}`$ and $`0.9\sqrt{R_0/4}`$, i.e weak, intermediate and strong magnetic toroidal intensity.

The two cases, $`R_0=1/4`$ and $`R_0=1`$ are more closely considered. In the first case, gravity alone cannot explain the collapse (see Eq. 18). The collapse is thus assisted by the toroidal magnetic field. In the second case, gravity is too strong and must be partially counterbalanced by the centrifugal force.

The results are respectively displayed in Fig. 4 and Fig. 5. The full line is $`b_\theta =0.1\sqrt{R_0/4}`$ (this case is not possible for $`R_0=1/4`$ according to Eqs. 18), the dotted line $`b_\theta =0.5\sqrt{R_0/4}`$ and the dashed line $`b_\theta =0.9\sqrt{R_0/4}`$.

From a comparison between Fig. 4 and 5, it is seen that the density is more peaked whereas the radial velocity is flatter for small densities ($`R_0=1/4`$) than for large one ($`R_0=1`$).

For $`R_0=1/4`$, the azimuthal velocity vanishes according to Eq. (18) when $`b_\theta =0.5\sqrt{R_0/4}`$. When, $`b_\theta =0.9\sqrt{R_0/4}`$, it is seen that the radial velocity and the toroidal field increase slightly just after the critical point, then reduce at $`X/X_c1.5`$.

Smaller values of density ($`R_0<1/4`$), require values of $`b_{\theta ,0}^{}`$ closer to $`\sqrt{R_0/4}`$ in order to have positive values of $`|v_\theta |`$ (see Eqs. 18). Similar behaviors to the previously presented cases are obtained.

For larger values of density ($`R_0>1`$), we find that with $`b_{\theta ,0}^{}=0`$, the solutions oscillate around the homologous solution, i.e. after the critical point the density decreases and increases alternately. When $`b_{\theta ,0}^{}0`$ (e.g. $`b_{\theta ,0}^{}>0.1`$) then the solutions cross another critical point and we stop the numerical integration. This behavior is somehow similar to the results obtained by Whitworth & Summers (1985) who find that the allowed values of the density at the origin in the Larson-Penston equations are quantized. It could thus be possible that larger values of $`R_0`$ are allowed in Eqs. (11)-(14). Only a systematic numerical search could answer this question.

### 4.3 Non-homologous cores

In this section, we investigate the case of non homologous cores, i.e. we start with initial conditions that do not satisfy the second conditions stated by Eqs. (18). Let us introduce the parameter, $`\delta v_\theta ^{}=v_\theta ^{}v_{\theta ,0}^{}`$, where $`v_{\theta ,0}^{}`$ is the value stated by Eqs. (18). The smaller $`|\delta v_\theta ^{}|`$ is, the closer we are to the homologous solutions.

We consider the two cases $`R_0=0.25`$ ; and $`R_0=0.5`$ with $`\mathrm{\Gamma }_z=1`$ and $`b_{\theta ,0}^{}=0.9\sqrt{R_0/4}`$ for three values of $`\delta v_\theta ^{}`$ ($`0.05,0.15,0.25`$ for $`R_0=0.25`$ and $`0.02,0.06,0.12`$ for $`R_0=0.5`$). The results are displayed in Figs. 6 and 7. Full lines represent the smallest value of $`|\delta v_\theta ^{}|`$, dotted line, the intermediate value, and dashed line the largest.

It is seen that the cores are no longer homologous and become less flat when $`|\delta v_\theta ^{}|`$ increases. Also, the transition between the sub-Alfvénic part and the super-Alfvénic one, is smoother. The poloidal force does not vanish any more in the cores and acts explicitly everywhere.

For $`R_0=0.25`$, the density and radial velocity values at $`5X_c`$ do not vary significantly with $`\delta v_\theta ^{}`$, whereas the values of the azimuthal velocity and the toroidal magnetic field decrease significantly when $`|\delta v_\theta ^{}|`$ increases.

For larger values of $`|\delta v_\varpi ^{}|`$ the solutions (for $`R_0=0.25`$ and $`R_0=0.5`$) cross another critical point. We consider them as not physical and we stop the numerical integration. It is then not possible for these two cases to obtain solutions with a vanishing azimuthal velocity as suggested by Eq. (18).

For larger values of $`R_0`$, the largest allowed value of $`\delta v_\varpi ^{}`$ decreases rapidly (for $`R_0=1`$, we find that the largest allowed $`|\delta v_\theta ^{}|`$ is between $`0.03`$ and $`0.04`$). This is consistent with the fact that the homologous solutions with large $`R_0`$ are not physical (see Sect. 4.2.5).

It is worth noting that the value of the radial velocity at $`5X_c`$ is around 2.5 to 3 which is lower than the values obtained previously. The smallest value of $`v_\varpi (5X_c)`$ we found is around 2 and is obtained for $`R_0<1/8`$, $`b_\theta \sqrt{R_0/4}`$, $`\mathrm{\Gamma }_z<1`$ and $`|v_{\theta ,0}^{}|<0.2`$.

This value is still 3 to 4 times higher than the value observed by Tafalla et al. (1999) in the starless cores L1544. However, it is known from previous study (e.g Foster & Chevalier 1993) that the Larson-Penston solution (in which the velocity tends to $`3.3C_s`$ at large $`r`$) is in good agreement with the inner collapsing part of the cloud, the external part being significantly different. This is due to the fact that the self-similar solutions make the assumption that the system is infinite and is self-similar everywhere. In particular, the self-similar assumption does not take into account the external medium. Self-similar solutions thus appear in the inner part of the cloud when the boundary conditions have been forgotten by the system.

## 5 Conclusion

In this paper, we have derived self-similar solutions able to describe the inner part of a collapsing, rotating, magnetized, self-gravitating, isothermal filaments. The filaments is collapsing homologously along the z-axis, the slope of the axial velocity being two times the slope of the radial one at the origin. The four ordinary equations obtained are similar to the equations derived by Penston (1969) and Larson (1969). They present a critical magnetosonic point that induces a bifurcation.

For some of the homologous solutions already obtained previously, (Aburihan et al. 2001, Hennebelle 2001) we carried out a study of the critical point and found that there is a critical value of $`b_{\theta ,0}^{}`$ (slope of the toroidal field at the origin) above which no bifurcation occurs. By studying the eigenvalues of the linearized equations in the neighbourhood of the critical point, we demonstrate that these solutions are able to cross the singularity since it is a node rather than a saddle point.

We have then explored the system numerically and obtained a series of density, velocity and magnetic fields profiles that could help to understand the observational data and could be used as benchmarks for a full numerical simulation.

The solutions obtained in this paper have some restrictions. First, they are valid not too far from the equatorial plane since they have an homologous axial velocity field that quickly diverges and not too far (few $`|X_c|`$) from the z-axis since the axial component of the velocity field does not depend on $`\varpi `$.

Nevertheless, these solutions are the first semi-analytical solutions describing the condensation of a magnetized rotating filament and having complex (non homologous) spatial profiles. They can be used for future analytical or numerical studies of the gravo-magnetic condensation.

## 6 acknowledgment

I thank Anthony Whitworth and Derek Ward-Thompson for stimulating discussions. I thank D. Moss, the referee, for his help in improving the manuscript. I gratefully acknowledge the support of an European Commission Research Training Network under the Fifth Framework Programme (No. HPRN-CT2000-00155).


format: text

text:
# Electron-Phonon Interaction in 𝑁𝑏𝐵₂ : A Comparison with 𝑀𝑔𝐵₂

## Abstract

We present a comparison of electron-phonon interaction in $`NbB_2`$ and $`MgB_2`$, calculated using full-potential, density-functional-based methods in $`P6/mmm`$ crystal structure. Our results, described in terms of (i) electronic structure, (ii) phonon density of states $`F(\omega )`$, (iii) Eliashberg function $`\alpha ^2F(\omega )`$, and (iv) the solutions of the isotropic Eliashberg gap equation, clearly show significant differences in the electron-phonon interaction in $`NbB_2`$ and $`MgB_2`$. We find that the average electron-phonon coupling constant $`\lambda `$ is equal to $`0.59`$ for $`MgB_2`$ and $`0.43`$ for $`NbB_2`$, leading to superconducting transition temperature $`T_c`$ of around $`22K`$ for $`MgB_2`$ and $`3K`$ for $`NbB_2.`$

The lack of success in finding superconductivity in other diborides with superconducting transition temperature, $`T_c`$, close to that of $`MgB_2`$  underscores the complex nature of interaction responsible for superconductivity in $`MgB_2.`$ In $`MgB_2`$ the complexity is further compounded by the presence of multi-faceted Fermi surface  and a highly anisotropic electron-phonon coupling, $`\lambda (𝐤,𝐤^{}),`$ over the Fermi surface . The dependence of superconducting properties on such details has ensured that we do not know, as yet, the exact nature of interaction leading to superconductivity in $`MgB_2.`$

Within Eliashberg-Migdal theory  of superconductivity, a reliable description of the superconducting state requires an accurate knowledge of $`\lambda (𝐤,𝐤^{})`$ and the renormalized electron-electron interaction, $`\mu ^{}`$, which are used as input to the fully anisotropic gap equation. The present computational capability allows us to evaluate $`\lambda (𝐤,𝐤^{})`$ accurately using density-functional-based methods but, unfortunately, $`\mu ^{}`$ cannot be evaluated. However, it is reasonable to assume that $`\mu ^{}`$ varies between $`0.1`$ to $`0.2`$ . Thus, the electron-phonon coupling $`\lambda (𝐤,𝐤^{})`$, which is a normal state function, must contain signatures of superconducting state.

In an attempt to identify some of the unique features of electron-phonon interaction in $`MgB_2`$ *vis-a-vis* other diborides we have studied (i) the electronic structure, (ii) the phonon density of states (DOS), (iii) the Eliashbrg function, and (iv) the solutions of the isotropic Eliashberg gap equation for $`NbB_2`$ and $`MgB_2`$ in $`P6/mmm`$ crystal structure.

The choice of $`NbB_2`$ has been motivated by the recent reports  of superconductivity, with $`T_c`$ going up to $`9.2K`$, under pressure in hole-doped $`Nb_xB_2`$. Earlier experiments have shown superconductivity in stoichiometric $`NbB_2`$  as well as Boron-enriched $`NbB_2`$  samples. The reported $`T_c`$ for stoichiometric $`NbB_2`$ varies from $`0.62K`$  to $`5.2K`$ , while for Boron-enriched $`NbB_{2.5}`$ the $`T_c`$ is found to be $`6.4K`$ . We also note that, recently, Kaczorowski *et al.*  did not find any superconductivity in $`NbB_2`$ down to $`2K.`$

We have calculated the electronic structure of $`NbB_2`$ and $`MgB_2`$ in $`P6/mmm`$ crystal structure with optimized lattice constants $`a`$ and $`c`$, as given in Table I. The lattice constants $`a`$ and $`c`$ were optimized using the ABINIT program  based on pseudopotentials and plane waves. For studying the electron-phonon interaction we used the full-potential linear response program of Savrasov , and calculated the dynamical matrices and the Hopfield parameter. These were then used to calculate the phonon DOS, $`F(\omega )`$, the electron-phonon coupling $`\lambda (𝐤,𝐤^{})`$, and the Eliashberg function, $`\alpha ^2F(\omega )`$, for $`NbB_2`$ and $`MgB_2`$. Subsequently, we have numerically solved the isotropic Eliashberg gap equation  for a range of $`\mu ^{}`$ to obtain the corresponding $`T_c.`$

Based on our calculations, described below, we find significant differences in the phonon DOS and the Eliashberg functions of $`NbB_2`$ and $`MgB_2`$. In particular, we find that the average electron-phonon coupling constant is equal to $`0.59`$ for $`MgB_2`$ and $`0.43`$ for $`NbB_2`$, leading to superconducting transition temperatures of around $`22K`$ for $`MgB_2`$ and $`3K`$ for $`NbB_2.`$

Before describing our results in detail, we provide some of the computational details of our calculation. The structural relaxation was carried out by the molecular dynamics program ABINIT  with Broyden-Fletcher-Goldfarb-Shanno minimization technique using Troullier-Martins pseudopotential  for $`MgB_2`$ and Hartwigsen-Goedecker-Hutter pseudopotential  for $`NbB_2`$, 512 Monkhorst-Pack  $`k`$-points and Teter parameterization for exchange-correlation. The kinetic energy cutoff for the plane waves was $`110Ry`$ for $`MgB_2`$ and $`140Ry`$ for $`NbB_2.`$ The charge self-consistent full-potential LMTO  calculations were carried out with the generalized gradient approximation for exchange-correlation of Perdew *et al*  and 484 $`k`$-points in the irreducible wedge of the Brillouin zone. For $`MgB_2`$, the basis set used consisted of $`3\kappa `$ panels and $`s,`$ $`p,`$ $`d`$ and $`f`$ orbitals at the $`Mg`$ site and $`s,`$ $`p`$ and $`d`$ orbitals at the $`B`$ site. In the case of $`NbB_2`$, we included $`2\kappa `$ panels and $`s`$, $`p`$ and $`d`$ orbitals at the $`Nb`$ site. In all cases the potential and the wave function were expanded up to $`l_{max}=6`$. The muffin-tin radii for $`Mg`$, $`B`$, and $`Nb`$ were taken to be $`2.4,`$ $`1.66,`$ and $`2.3`$ atomic units, respectively.

The calculation of dynamical matrices and the Hopfield parameters for $`MgB_2`$ were carried out using a $`6\times 6\times 6`$ grid while for $`NbB_2`$ we used a $`4\times 4\times 4`$ grid resulting in $`28`$ and $`12`$ irreducible $`𝐪`$-points, respectively. For Brillouin zone integrations in $`MgB_2`$ we used a $`6\times 6\times 6`$ grid while for $`NbB_2`$ we used $`8\times 8\times 8`$ grid of $`𝐤`$-points. The Fermi surface was sampled more accurately with a $`24\times 24\times 24`$ grid of $`𝐤`$-points using the double grid technique as outlined in Ref. .

Here, we like to point out the reasons for carrying out the linear response calculation for $`MgB_2`$ in spite of earlier calculations by Kong *et al*.  and Choi *et al.* . The linear response calculation by Kong *et al*. is similar to the present approach, while Choi *et al*. used pseudopotentials and frozen phonon method to evaluate the electron-phonon coupling $`\lambda (𝐤,𝐤^{})`$. The present approach differs from the work of Kong *et al* in the selection of q-points and the Brillouin zone integrations. As a result Kong *et al*. find the average electron-phonon coupling constant $`\lambda =0.87\pm 0.05`$, which is much higher than the value of $`\lambda =0.61`$ as reported by Choi *et al*. , as well as the experimental values of $`0.58`$  and $`0.62`$  as deduced from specific heat measurements. Our calculated value of $`\lambda `$ is equal to $`0.59,`$ in close agreement with the work of Choi *et al*..

A comparison of site- and $`l`$-resolved electronic density of states of $`NbB_2`$  and $`MgB_2`$  at the Fermi energy is given in Table II. A further decomposition of the densities of states ($`st/Ry`$) in terms of cubic harmonics reveals the dominance of $`B`$ $`p`$ electrons at the Fermi energy in $`MgB_2`$ $`(p_{x(y)}=0.9,p_z=1.55)`$ than in $`NbB_2`$ $`(p_{x(y)}=0.43,p_z=0.48)`$. In addition, in $`NbB_2`$, the $`Nb`$ $`d`$-electrons $`(d_{xy(x^2y^2)}=1.83,d_{yz(zx)}=1.66,d_{3z^21}=2.56)`$ are present in substantial amount, indicating a more active role for $`Nb`$ in determining the possible superconducting properties of these materials than played by $`Mg`$ in $`MgB_2.`$

In Fig. 1 we show the phonon DOS $`F(\omega )`$ of $`NbB_2`$ and $`MgB_2`$ calculated using the full-potential linear response program as described earlier. For $`MgB_2`$ we can clearly identify four significant peaks in the phonon DOS at $`33,`$ $`53,`$ $`79`$ and $`96`$ $`meV`$, respectively. The peak at $`33`$ $`meV`$ is related to the van Hove singularity in the acoustical mode , and it involves the motion of $`Mg`$ atom and $`B`$ atoms separately. However, the region around the peak at $`53`$ $`meV`$ results from the motion of both $`Mg`$ and $`B`$ atoms. The phonon DOS around $`79`$ and $`96`$ $`meV`$ peaks are due to the coupled motion of $`BB`$ atoms in the $`xy`$ plane. In particular, the peak at $`79`$ $`meV`$ corresponds to the in-plane $`BB`$ bond stretching mode, and in Ref.  it is located at $`77meV`$. Similarly, for $`NbB_2`$, the peak in $`F(\omega )`$ at $`32`$ $`meV`$ is dominated by the motion of $`Nb`$ atom, while the region around $`6570`$ $`meV`$ results from the coupled motion of $`Nb`$ and the two $`B`$ atoms. Not surprisingly, $`81`$ $`meV`$ peak in the phonon DOS of $`NbB_2`$ corresponds to the in-plane $`BB`$ motion. In contrast with the phonon DOS in $`MgB_2,`$ the phonon DOS in $`NbB_2`$ around $`106`$ $`meV`$ results from the displacements of both $`Nb`$ and the two $`B`$ atoms.

To see the strengths with which the different modes of the ionic motion couple to the electrons, and thus are capable of influencing the superconducting properties the most, we show in Fig. 2 the Eliashberg function $`\alpha ^2F(\omega )`$ of $`NbB_2`$ and $`MgB_2`$ calculated as described earlier. The most striking feature of Fig. 2 is the overall strength of the electron-phonon coupling in $`MgB_2`$ as compared to $`NbB_2`$. We find that the average electron-phonon coupling constant $`\lambda `$ is equal to $`0.59`$ for $`MgB_2`$ and $`0.43`$ for $`NbB_2`$, which clearly shows that $`MgB_2`$ is more likely to show superconductivity with a higher $`T_c`$ than $`NbB_2.`$

Further analysis of the Eliashberg function, as shown in Fig. 2, reveals the importance of the in-plane $`BB`$ bond-stretching optical phonon mode in $`MgB_2`$, which gives rise to the dominant peak at $`79meV`$. The other peaks in the phonon DOS of $`MgB_2`$, such as the peaks at $`33`$ and $`53`$ $`meV`$, couple weakly with the electrons at the Fermi energy. Thus the motion of $`Mg`$ atom plays a relatively insignificant role in determining the superconducting properties of $`MgB_2.`$ In contrast, in the case of $`NbB_2`$ the phonon modes with peaks at $`32`$, $`81`$ and $`106`$ $`meV`$ couple to the electrons with almost equal strength, albeit much smaller than in $`MgB_2,`$ as can be seen from Fig. 2. We could have expected this because of the significant presence of the $`Nb`$ $`d`$ electrons at the Fermi energy. In Table III we have listed the Hopfield parameter $`\eta ,`$ the electron-phonon coupling constant $`\lambda `$, and the various averages of the phonon frequencies for $`NbB_2`$ and $`MgB_2`$. The values listed in Table III for $`MgB_2`$ are in good agreement with the corresponding results of Choi *et al.* .

To examine the superconducting transition temperature, *if* any, of $`NbB_2`$ and $`MgB_2`$ we have used the calculated Eliashberg function $`\alpha ^2F(\omega )`$ to solve numerically the isotropic gap equation , and the results are shown in Fig. 3 for a range of values of $`\mu ^{}.`$ From Fig. 3 we find that for $`\mu ^{}=0.1`$ the $`T_c`$ for $`MgB_2`$ is equal to $`23K`$, while for $`NbB_2`$ it is equal to $`4K`$. Thus, our calculation shows that $`NbB_2`$ is superconducting with a possible $`T_c`$ of around $`3K.`$ It is worthwhile to point out that to obtain a $`T_c`$ close to that of $`39K`$ for $`MgB_2`$, as found experimentally , one has to solve the *anisotropic* gap equation . The need to solve the anisotropic gap equation for $`MgB_2`$ rather than the isotropic gap equation arises due to the highly anisotropic electron-phonon coupling $`\lambda (𝐤,𝐤^{})`$  over the Fermi surface. In the case of $`NbB_2`$ the electron-phonon coupling is *neither as strong nor as anisotropic*, and thus the results obtained with the isotropic gap equation are reliable.

As indicated earlier, the experiments show superconductivity in hole-doped $`Nb_xB_2`$ under pressure and Boron-enriched $`NbB_2.`$ Hole-doping and Boron enriching both lead to a relative increase in the Boron population at the Fermi energy and, probably, enhances the peak around $`81meV`$ leading to an increase in $`T_c.`$ Of course a more quantitative investigation is needed to pinpoint the exact nature of changes in $`NbB_2`$ which lead to superconductivity.

In conclusion, we have studied the electron-phonon interaction in $`NbB_2`$ and $`MgB_2`$, using full-potential, density-functional-based methods in $`P6/mmm`$ crystal structure. We have described our results in terms of (i) electronic structure, (ii) phonon density of states, (iii) Eliashberg function, and (iv) the solutions of the isotropic Eliashberg gap equation, which clearly show significant differences in the electron-phonon interaction in $`NbB_2`$ and $`MgB_2`$. We find that the average electron-phonon coupling constant is equal to $`0.59`$ for $`MgB_2`$ and $`0.43`$ for $`NbB_2`$, leading to superconducting transition temperature of around $`22K`$ for $`MgB_2`$ and $`3K`$ for $`NbB_2.`$


text:
# Atom-Atom Scattering Under Cylindrical Harmonic Confinement: Numerical and Analytic Studies of the Confinement Induced Resonance

## Abstract

M. Olshanii \[Phys. Rev. Lett. 81, 938 (1998)\] recently solved the atom-atom scattering problem with a pseudopotential interaction in the presence of transverse harmonic confinement, i.e. within an ‘atom waveguide’, deriving an effective one-dimensional coupling constant that diverged at a “confinement induced resonance” (CIR). Here, we report numerical results for finite range potentials that corroborate this resonance. In addition, we now present a physical interpretation of this effect as a novel type of Feshbach resonance in which the transverse modes of the waveguide assume the roles of ‘open’ and ‘closed’ scattering channels.

Recently, there have been a number of experiments in which ultracold atoms and/or Bose-Einstein condensates have been loaded into magnetic or optical ‘atom waveguides’ Bongs ; Hansch ; Denschlag ; Hinds ; Prentiss ; Greiner ; Arlt ; DWeiss . One goal of such experiments is to reach the ‘single-mode’ or quasi-1D regime, where only the ground state of transverse motion is significantly populated at thermal equilibrium. This regime is of great practical interest due to the potential for ultra-sensitive rotation and gravitational gradient detection with guided single-mode atom interferometers. In addition to such applications, reaching the quasi-1D regime is of significant theoretical interest as the 1D delta-interacting boson gas represents one of the few known fully integrable quantum field theories. In a finite system with infinitely strong (hard-core) delta-function interactions, boson many-body states in 1D have been shown to correspond via a one-to-one mapping with the highly-correlated states of the corresponding non-interacting Fermi gas Girardeau . The properties of this Tonks-Girardeau gas have been a topic of significant current theoretical interest Tonks ; Girardeau ; GWT ; DLO ; Petrov ; Stoof ; DGW  in anticipation of future atom-waveguide experiments. Additionally, the homogeneous 1D Bose with arbitrary-strength delta-function interactions, known as the Lieb Liniger model, is also a fully integrable system LL .

To make the connection between experiments in tightly confining waveguides and theoretical models in 1D, it is necessary to know the relationship between the effective 1D coupling constant, $`g_{1D}`$, and the 3D scattering length, $`a`$. This problem was first addressed rigorously in MOPRL , where it was predicted that a ‘confinement induced resonance’ (CIR) modifies the effective interaction, resulting in an effective 1D coupling strength which can be tuned from $`\mathrm{}`$ to $`+\mathrm{}`$ by varying the transverse width of the waveguide, $`a_{}`$ over a small range in the vicinity of the resonance at $`a_{}=Ca`$, where $`C=1.4603\mathrm{}`$. Hence this resonance clearly has significant implications for atom-waveguide experiments. Until now, however, there has been no convincing physical explanation for the effect, thus raising questions concerning its appearance in systems with finite-range interactions.

The primary goal of this Letter, therefore, is to present numerical calculations of scattering in the presence of a cylindrical harmonic potential using finite-range atom-atom potentials, confirming the existence of the CIR. In addition, we provide a much-needed physical interpretation of the effect as a Feshbach-type resonance involving bound states of the energetically closed transverse modes (‘channels’). Analogous numerical scattering studies in spherically symmetric traps have been reported Tiesingapp ; Bolda , while scattering in the presence of harmonic confinement in one dimension is discussed in PHS ; PS .

We begin our analysis by considering a collision between two atoms initially in the ground state of transverse motion. In the presence of harmonic transverse confinement the center-of-mass motion and relative motions are separable, with the wavefunction of the relative coordinate satisfying an effective single-particle model with a stationary scatterer at the origin. Furthermore, if the longitudinal kinetic energy in the center-of-mass frame is less than the transverse level spacing then the atoms remain asymptotically frozen in the ground state. Low-energy scattering in this regime can then be modeled in the pseudopotential approximation Demkov  by the 1D Hamiltonian

$$H_{1D}=\frac{\mathrm{}^2}{2\mu }\frac{^2}{z^2}+g_{1D}\delta (z),$$
(1)
where $`z`$ is the longitudinal atomic separation and $`\mu `$ is the reduced mass. In the pseudopotential approximation, $`g_{1D}`$ is typically obtained as in DGW  by assuming that the wavefunction of the relative coordinate $`𝐫=𝐫_1𝐫_2=z\widehat{z}+\rho \widehat{\rho }`$ factorizes as $`\mathrm{\Psi }(𝐫)=\varphi _0(\rho )\psi (z)`$, where $`\varphi _0(\rho )=\mathrm{exp}(\rho ^2/2a_{}^2)/(a_{}\sqrt{\pi })`$ is the transverse ground state, and $`a_{}=\sqrt{\mathrm{}/\mu \omega _{}}`$ is the transverse harmonic oscillator length for the relative atomic motion, $`\omega _{}`$ being the transverse trap frequency. The effective 1D potential is then defined via $`2\pi \rho 𝑑\rho |\varphi _0(\rho )|^2\frac{2\pi \mathrm{}^2a}{\mu }\delta ^3(𝐫)=g_{1D}\delta (z)`$, where $`a`$ is the 3D scattering length, which leads to

$$g_{1D}_0^{\mathrm{}}2\pi \rho 𝑑\rho =\frac{2\mathrm{}^2a}{\mu a_{}^2}.$$
(2)
From an exact solution of the three-dimensional scattering problem with a zero-range s-wave interaction, however, it was recently predicted MOPRL  that a “confinement induced resonance” (CIR) modifies the effective interaction. By matching the low-energy scattering amplitude of the exact solution to that of Eq. (1), it was found that:

$$g_{1D}=\frac{2\mathrm{}^2a}{\mu a_{}^2}\frac{1}{(1Ca/a_{})},$$
(3)
where $`C=\zeta (1/2)=1.4603..`$. We note that the use of Eqs. (1) and (3) requires $`ka_{}1`$, where $`\mathrm{}k`$ is the collision momentum. The appearance of the resonance term in the denominator implies that to obtain an infinite delta-function interaction (thereby accessing the Tonks-Girardeau regime), it is sufficient to satisfy the CIR condition $`a_{}Ca`$. From the naive formulation (2) one would conclude that accessing this regime requires the more extreme condition $`ana_{}^2`$, where $`n`$ is the linear density.

Before addressing the physical interpretation of the CIR, we first describe our numerical results pertaining to low-energy scattering and the effective 1D coupling constant with finite-range interactions. The relation between $`g_{1D}`$ and the 1D scattering amplitude may be found by assuming that the scattering eigenstates $`\mathrm{\Psi }(z)`$ of Hamiltonian (1) take the form given in MOPRL

$$\mathrm{\Psi }(z)=e^{ikz}+f_ee^{ik|z|};E=\mathrm{}^2k^2/2\mu ,$$
(4)
where $`f_e`$ is the coefficient for the even part of $`\mathrm{\Psi }`$, while the odd part can be shown to vanish by continuity arguments. From Eq. (4), we find $`\mathrm{\Psi }^{\prime \prime }(z)=k^2\mathrm{\Psi }(z)+2ikf_e(k)\delta (z)`$, which must equal $`(2\mu /\mathrm{}^2)g_{1D}\mathrm{\Psi }(z)`$ in the limit $`k0`$. This gives

$$g_{1D}=\underset{k0}{lim}\frac{\mathrm{}^2kf_{eR}(k)}{\mu f_{eI}(k)}$$
(5)
where $`f_{eR}`$ and $`f_{eI}`$ are the real and imaginary parts of $`f_e`$. To extract analogous values for $`g_{1D}`$ from numerical scattering calculations, we obtain eigenfunctions of the Hamiltonian $`\widehat{H}=\widehat{H}_z+\widehat{H}_{}+\widehat{V}`$, where

$`\widehat{H}_z={\displaystyle \frac{\mathrm{}^2}{2\mu }}{\displaystyle \frac{^2}{z^2}};\widehat{H}_{}={\displaystyle \frac{\mathrm{}^2}{2\mu }}\left[{\displaystyle \frac{^2}{\rho ^2}}+{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{}{\rho }}\right]+{\displaystyle \frac{\mu }{2}}\omega _{}^2\rho ^2.`$ (6)
In the present work we restrict ourselves to the case of zero azimuthal angular momentum $`m`$=0, as there is negligible s-wave scattering for $`m0`$. For the atom-atom potential, we will study two cases: $`V(r)=C_{12}/r^{12}C_6/r^6(r^2=z^2+\rho ^2)`$ and the spherical square-well, $`V(r)=\overline{V}S(br)`$, $`S(r)`$ being the unit step function.

With application to the case of Cs atoms in a 1D optical well DWeiss  in mind, we consider atoms with the mass of 133Cs, and with the $`C_6`$ coefficients as determined recently Leo  to be $`C_6`$(Cs) = 6890 a.u. For Cs atoms in the $`(F,M)`$ = (3,3) state, the scattering length is determined only as an upper bound, $`a(Cs)<`$ -140 nm, as compared with $`a_{}`$ = 29.5 nm reported in citeDWeiss. In order to study a more general situation, we allow $`C_{12}`$ to vary. To simplify the numerics we consider the regime of just 1 to 3 $`J=0`$ bound states, rather than the 47 bound states in the actual Cs2 $`{}_{}{}^{3}\mathrm{\Sigma }_{u}^{+}`$ state, and we neglect other terms in the dispersion potential. The 6-12 potential may be characterized by $`R_e`$, the minimum of the potential well, where $`R_e=(2C_{12}/C_6)^{1/6}`$. For $`C_6=C_6`$(Cs) and for a series of values of $`C_{12}`$ we obtain the free-space $`s`$-wave scattering length, $`a`$, by Numerov integration of the Schrödinger equation. When a resonance state passes through threshold $`a`$ exhibits a simple pole, thus the 6-12 potential provides the full range of 3D scattering lengths. For the second case of a spherical well potential the scattering length for a well of depth $`\overline{V}`$ and range $`b`$ is given by $`a=b\mathrm{tan}(b\eta )/\eta `$, where $`\eta ^2=2\mu \overline{V}/\mathrm{}^2`$.

To solve the scattering problem for such central potentials plus a transverse harmonic potential, we employed a numerical mesh in $`\rho `$ and $`z`$, and found eigenfunctions in a cylindrical box of finite length. The box was sufficiently long in $`z`$ that the asymptotic form of the wavefunction as $`|z|\mathrm{}`$ could be determined. The eigenfunctions were of odd or even parity, with the odd parity functions exhibiting negligible scattering effects. For even functions of energy $`E=\mathrm{}\omega _{}+\mathrm{}^2k^2/2\mu `$, the asymptotic form is

$$\mathrm{\Psi }(\rho ,z)\stackrel{|z|\mathrm{}}{}𝒩\left[(1+f_{eR})\mathrm{cos}(kz)f_{eI}\mathrm{sin}(k|z|)\right]\varphi _0(\rho ).$$
(7)
Values for $`f_{eR}`$, $`f_{eI}`$ and $`𝒩`$ were extracted from the coefficients of $`\mathrm{sin}(kz)`$ and $`\mathrm{cos}(kz)`$ in Eq. (7), determined over a range of $`z`$ values for which $`V(r)`$ is negligible, and from conservation of probability current. From $`f_{eR}`$ and $`f_{eI}`$, $`g_{1D}`$ is then obtained from Eq. (5) via extrapolation of finite $`k`$ data to $`k=0`$. The numerical mesh was provided by the discrete variable representation (DVR). In order to increase the density of points near $`z`$=0, the $`z`$ coordinate was scaled by $`z=U(y)`$, where $`U(y)=a\mathrm{cosh}(z/b)`$, analogous to scaling used in Tiesinga . A uniform mesh in $`y`$ was then used Colbert  along with a Laguerre DVR in $`\rho `$ BayeHeenan . Eigenfunctions in first iteration were found by exact diagonalization over a relatively small mesh ($`<`$10,000 mesh points). As additional mesh points were added, sparse matrix diagonalization techniques were used.

Numerical results for $`g_{1D}`$ with the 6-12 potential and with the spherical well are shown in Fig. 1, in comparison with the analytic result (3). The numerical calculations clearly exhibit a singularity in $`g_{1D}`$ at $`a/a_{}=1/C`$. For $`a`$ less than this value the calculated $`g_{1D}`$ values agree well with the analytic expression. There do appear to be systematic deviations between the exact numerical results and the pseudopotential theory that are greater for the 6-12 potential than for the spherical well for $`a/a_{}>1/C`$. This effect is connected to the Feshbach resonance interpretation and occurs because the bound state of the pseudopotential does not agree precisely with the upper bound-state of the finite-range potentials having the same scattering length $`a`$.

These numerical results validate the pseudopotential analytic result MOPRL , demonstrating that the CIR is indeed a physical phenomena rather than an artifact of the pseudopotential approximation. It is therefore important to establish the physical origin of the CIR. In what follows we are going to show that the CIR is in fact a zero-energy Feshbach resonance, occurring when the energy of a bound state of the asymptotically closed channels (i.e. the excited transverse modes) coincides with the continuum threshold of the open channel (lowest transverse mode).

To illustrate this interpretation we consider the total Hamiltonian $`\widehat{H}=\widehat{H}_z+\widehat{H}_{}+\widehat{V}_{pseudo}`$, where the pseudopotential operator is defined by $`𝐫|\widehat{V}_{pseudo}|\psi =\frac{2\pi \mathrm{}^2a}{\mu }\delta ^3(𝐫)\frac{}{r}r𝐫|\psi `$. Because our interpretation of the CIR relies heavily on the behavior of two-atom bound states, it is important to recall that the Huang-Fermi pseudopotential, $`\widehat{V}_{pseudo}`$, supports a single bound state in free space with energy $`\mathrm{}^2/(\mu a^2)`$, provided $`a>0`$. With the existence of this bound state in mind, we proceed by formally splitting the Hamiltonian onto ‘ground’, ‘excited’, and ‘ground-excited coupling’ parts according to

$`\widehat{H}`$ $`=`$ $`\widehat{H}_g+\widehat{H}_e+\widehat{H}_{ge}`$ (8)
$`=`$ $`\widehat{P}_g\widehat{H}\widehat{P}_g+\widehat{P}_e\widehat{H}\widehat{P}_e+(\widehat{P}_e\widehat{H}\widehat{P}_g+\text{h.c.}),`$
where $`\widehat{P}_g=|00|`$, $`\widehat{P}_e=_{n=1}^{\mathrm{}}|nn|`$, are the corresponding projection operators, $`|n`$ being the eigenstate of the transverse two-dimensional harmonic oscillator with radial quantum number $`n`$ and zero axial angular momentum (we recall that in defining $`\widehat{H}_{}`$ we have already projected out the states with non-zero axial angular momentum). The corresponding eigenvalues of the transverse Hamiltonian are given by $`\widehat{H}_{}|n,=\mathrm{}\omega _{}(2n+1)|n`$.

The ‘ground’ Hamiltonian has a 1-d coordinate representation of the form of Eq. (1), corresponding to the motion of a one-dimensional particle in presence of a $`\delta `$-barrier with the ‘bare’ coupling constant given by Eq. (2). The spectrum of $`\widehat{H}_g`$ is continuous for energies above the threshold energy $`E_{C,g}=\mathrm{}\omega _{}`$ ($`=0`$) and contains a single bound state at energy $`E_{B,g}<E_{C,g}`$ for the case of positive 3D scattering length. Likewise, the spectrum of the ‘excited’ Hamiltonian is clearly continuous for energies $`E_{C,e}=3\mathrm{}\omega _{}`$ ($`=1`$), but as we will see below, $`\widehat{H}_e`$ supports one bound state of an energy $`E_{B,e}<E_{C,e}`$ for all values of the 3D scattering length, $`a`$.

It is clear that the remaining coupling between the ‘ground’ and ‘excited’ subspaces will renormalize the bare coupling (2). Furthermore, according to the Feshbach scheme, one would predict a resonance in the renormalized $`g_{1D}`$ for a set of parameters, such that the energy of the bound state of $`\widehat{H}_e`$ coincides with the continuum threshold of $`\widehat{H}_g`$:

$`E_{B,e}=E_{C,g}\text{CIR}`$ (9)
As we will see below this scheme indeed predicts a position of the CIR exactly.

The energy $`E_{B,e}`$ of the bound state of $`\widehat{H}_e`$ can be found using the following two step procedure. First we identify the bound state energy of the full Hamiltonian $`\widehat{H}`$ as a pole of the scattering amplitude on the physical Riemann sheet. Second, we make use of the peculiar property of the two-dimensional harmonic oscillator that ‘excited’ Hamiltonian $`\widehat{H}_e`$and the full Hamiltonian $`\widehat{H}`$ can be transformed to each other via a simple unitary transformation. This allows for simple relation between their energy spectra, in particular the energy of the bound states.

The even-wave one-dimensional scattering amplitude $`f_e`$ at an energy $`E`$ ($`E_{c,g}E<E_{c,e}`$), as defined in (4), has been derived in MOPRL . We have now obtained a closed-form analytic expression for the previously derived result, so that the scattering amplitude can be expressed as

$$f_e(k)=\frac{2i}{ka_{}[\frac{a_{}}{a}+\zeta (1/2,(ka_{}/2)^2)]}$$
(10)
where

$`\zeta (s,\alpha )=\underset{N\mathrm{}}{lim}{\displaystyle \underset{n=0}{\overset{N}{}}}{\displaystyle \frac{1}{(n+\alpha )^s}}{\displaystyle \frac{(N+\alpha )^{s+1}}{s+1}};`$
$`Re(s)>0,Im(s)=0,`$
$`z^s=|z|^se^{is(\text{Arg}(z)2\pi )},\mathrm{\hspace{0.17em}\hspace{0.17em}0}<\text{Arg}(z)2\pi `$ (11)
is the Hurwitz Zeta function Apostol , and the wave-vector $`k`$ is given by $`E=E_{C,g}+\mathrm{}^2k^2/2\mu `$. The bound state energies of the full Hamiltonian $`\widehat{H}`$ will be given by the poles, $`\overline{k}`$, on the positive imaginary axis of the analytic continuation of $`f_e(k)`$: $`E_{B,full}=\mathrm{}^2Im^2(\overline{k})/2\mu `$. One can see that in order to avoid crossing the branch-cuts of the Zeta function, the continuation should be performed inside the $`0\text{Arg}(k)\pi /2`$ quadrant of the complex plane. We find a single pole corresponding to the following implicit equation for the bound state energy: $`\zeta (1/2,E_{B,full}/(2\mathrm{}\omega _{})+1/2)=a_{}/a`$.

Notice now, that the full Hamiltonian $`\widehat{H}`$ and the excited Hamiltonian $`\widehat{H}_e`$ are connected via a simple transformation: $`\widehat{H}_e=\widehat{A}^{}\widehat{H}\widehat{A}`$, where $`\widehat{A}^{}=_{n=0}^{\mathrm{}}|n+1,n|`$. Note that both $`\widehat{H}_e`$ and $`\widehat{H}`$ include interactions, thus the above property is highly nontrivial and stems from the fact that the $`m=0`$ eigenfunctions of the two-dimensional harmonic oscillator all have the same value, $`1/\sqrt{\pi }a_{}`$, at the origin. Thus the three dimensional $`\delta `$-interaction has the same matrix elements between all the harmonic oscillator states so that interaction matrix is unaffected by the shift operator.

From the above we conclude that the bound state energy $`E_{B,e}`$ of the ‘excited’ Hamiltonian is related to the bound state energy $`E_B`$ of the full Hamiltonian via $`E_{B,e}=E_B+2\mathrm{}\omega _{}`$, and thus obeys an implicit equation $`\zeta (1/2,E_{B,e}/(2\mathrm{}\omega _{})+3/2)=a_{}/a`$. The CIR condition (9) can be now explicitly formulated as $`\zeta (1/2,E_{C,g}/(2\mathrm{}\omega _{})+3/2)=a_{}/a`$. Using $`E_{C,g}=\mathrm{}\omega _{}`$, we finally arrive at the exact CIR condition $`\frac{a_{}}{a}=\zeta (1/2,1)=\zeta (1/2,0)=C`$. A similar effect is associated with resonance behavior in harmonically confined 2D scattering for $`a<0`$ PHS . We note that this resonance would most-likely be observed via changes in the macroscopic properties of the ground-state of a many-atom system, i.e. density distribution, as described in DLO .

This Feshbach scheme is illustrated in Figure 2, where we plot the bound-state energies $`E_{B,e}`$ (dark solid line) and $`E_B`$ (thin solid line) as a function of the ratio $`a_{}/a`$. The continuum thresholds $`E_{C,e}`$ and $`E_{C,g}`$ are also indicated, illustrating that the CIR occurs when the bound state of the manifold of closed channels, $`E_{B,e}`$, crosses the continuum threshold of the open channel, $`E_{C,g}`$. In addition, we have plotted the bound state energies of the full Hamiltonian as determined numerically for the $`612`$ and spherical well potentials, showing good agreement with the pseudopotential result. As the bound-state energy deepens, we start to see quantitative disagreement between the bound state energies of the finite range potentials and the pseudopotential. This disagreement is consistent with the discrepancy in the position of the CIR shown in Fig. 1, showing that it is the bound state energy, and not the scattering length which determines the location of the CIR. Lastly, we note that while in free space a weakly bound state exists only for $`a>0`$, we see that in the waveguide such a state exists for all $`a`$. These bound states may be of significant interest, allowing the formation of dimers via a modulation of the waveguide potential at the frequency $`(E_{C,g}E_B)/\mathrm{}`$. This may lead to an atom-waveguide based scheme for forming ultracold dimer molecules, as well as the possibility to use molecular spectroscopy as a sensitive probe of the atomic field inside the waveguide.

This work was supported by the National Science Foundation through a grant for the Institute for Theoretical Atomic and Molecular Physics at Harvard University and Smithsonian Astrophysical Observatory. The work at Stony Brook was supported by the Office of Naval Research, and the work at USC by the NSF grant PHY-0070333. We thank D. S. Weiss for discussions that provided part of the motivation for this study and R. Shakeshaft for discussions on Feshbach resonances.


text:
# Coherent processing of a light pulse stored in a medium of four-level atoms

## Abstract

It is demonstrated that the properties of light stored in a four-level atomic system can be modified by an additional control interaction present during the storage stage. By choosing the pulse area of this interaction one can in particular continuously switch between two channels into which light is released.
PACS numbers: 42.50.Gy, 03.67.-a

In recent years a number of nonlinear optical phenomena in weak fields have been intensively investigated. They are connected with a coherent excitation of an atomic medium, the optical properties of which can then be drastically modified. An important example is the electromagnetically induced transparency , possibly additionally controlled in time, which is manifested as light slow-down or even its storage and a controlled release . The simplest realization of those effects occurs in atomic systems with three active states in a $`\mathrm{\Lambda }`$ configuration. Extending such a configuration by adding a coherently coupled fourth level opens new possibilities of an external control of such processes  . In our previous paper  we have shown that in a double $`\mathrm{\Lambda }`$ system it is possible to change the light frequency of the stored light or even to release two pulses of different frequencies by applying two control fields, properly chosen and delayed in time.

Light storing in the form of atomic coherences joins the advantages of the efficiency of light as an information carrier and of an atomic medium as an information store. Thus it might be a question of a practical importance how to modify in a controlled way the properties of the released light by processing the atomic medium during the storage stage. In the present work we investigate two possibilities of controlling the released pulse or pulses by modifying the atomic coherence due to the stopped light. In the case (a) the lower, initially empty state of a typical $`\mathrm{\Lambda }`$ system is additionally coupled to a fourth state by another laser. In the case (b) the initially occupied state of the $`\mathrm{\Lambda }`$ system is coupled by some kind of an effective interaction (e.g., magnetic or two-photon electric coupling) with another state of the same parity. We show that the Rabi oscillations due to the new interactions modify in a coherent way the properties of the released light.

We consider a quasi one-dimensional medium of four-level atoms with three lower metastable states $`|b>`$ and $`|c>`$ and $`|d>`$ and an upper state $`|a>`$ (Fig.1). The position of an atom is described by the variable $`z`$, which is considered continuous. The states $`b,a`$ and $`c`$ constitute a typical $`\mathrm{\Lambda }`$ system with the weak signal field 1 and a strong control field 2. In the case (a) the state $`c`$ is additionally coupled with another state $`d`$ by a laser field 4. In the case (b) the state $`a`$ is coupled with a fourth state $`d`$ by a weak signal field 3 while the states $`b`$ and $`d`$ are coupled by some effective coupling 4. The interaction Hamiltonian in the case (a) is $`V=\widehat{d}_{j=1,2,4}ϵ_j\mathrm{cos}\varphi _j`$, with $`\varphi _j=\omega _jtk_jz`$, $`ϵ_j=ϵ_j(z,t)`$ being slowly varying envelopes and $`\widehat{d}`$ \- the dipole moment operator; we have assumed that all the fields have the same linear polarization. In the case (b) the Hamiltonian reads $`V=\widehat{d}_{j=1,2,3}ϵ_j\mathrm{cos}\varphi _j+[(iU)|b><d|(iU)|d><b|]\mathrm{cos}\varphi _4`$, where for numerical reasons we have made a physically insignificant assumption that the matrix element $`iU`$ of the effective interaction 4 is imaginary. The matrix elements of the dipole moment $`d_1=(\widehat{d})_{ab}`$, $`d_2=(\widehat{d})_{ac}`$, $`d_3=(\widehat{d})_{ad}`$, $`d_4=(\widehat{d})_{cd}`$ are taken real. Resonant conditions concerning all the couplings are assumed.

The evolution equation $`i\mathrm{}\dot{\rho }=[H,\rho ]`$ for the density matrix $`\rho =\rho (z,t)`$ for an atom at position $`z`$, after making the rotating-wave approximation, transforming-off the rapidly oscillating factors: $`\rho _{ab}=\sigma _{ab}\mathrm{exp}(i\varphi _1)`$, $`\rho _{ac}=\sigma _{ac}\mathrm{exp}(i\varphi _2)`$, $`\rho _{bc}=\sigma _{bc}\mathrm{exp}[i(\varphi _1\varphi _2)]`$, $`\rho _{db}=\sigma _{db}\mathrm{exp}(i\varphi _4)`$, $`\rho _{dc}=\sigma _{dc}\mathrm{exp}[i(\varphi _2\varphi _3)]`$, $`\rho _{ad}=\sigma _{ad}\mathrm{exp}[i(\varphi _3)]`$, $`\rho _{ii}=\sigma _{ii}`$, and after adding relaxation terms describing the spontaneous emission within the system, takes the form in the more complicated case (b)

$`i\dot{\sigma }_{aa}`$ $`=`$ $`{\displaystyle \frac{1}{2\mathrm{}}}ϵ_1d_1(\sigma _{ba}\sigma _{ab})+{\displaystyle \frac{1}{2}}\mathrm{\Omega }_2(\sigma _{ca}\sigma _{ac}){\displaystyle \frac{1}{2\mathrm{}}}ϵ_3d_3(\sigma _{da}\sigma _{ad})i(\mathrm{\Gamma }_b^a+\mathrm{\Gamma }_c^a+\mathrm{\Gamma }_d^a)\sigma _{aa},`$ (1)
$`i\dot{\sigma }_{bb}`$ $`=`$ $`{\displaystyle \frac{1}{2\mathrm{}}}ϵ_1d_1(\sigma _{ab}\sigma _{ba})+{\displaystyle \frac{1}{2\mathrm{}}}iU(\sigma _{bd}+\sigma _{db})+i\mathrm{\Gamma }_b^a\sigma _{aa},`$ (2)
$`i\dot{\sigma }_{cc}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{\Omega }_2(\sigma _{ac}\sigma _{ca})+i\mathrm{\Gamma }_c^a\sigma _{aa},`$ (3)
$`i\dot{\sigma }_{dd}`$ $`=`$ $`{\displaystyle \frac{1}{2\mathrm{}}}ϵ_3d_3(\sigma _{ad}\sigma _{da}){\displaystyle \frac{1}{2\mathrm{}}}iU(\sigma _{bd}+\sigma _{db})+i\mathrm{\Gamma }_d^a\sigma _{aa},`$ (4)
$`i\dot{\sigma }_{ab}`$ $`=`$ $`{\displaystyle \frac{1}{2\mathrm{}}}ϵ_1d_1(\sigma _{bb}\sigma _{aa})+{\displaystyle \frac{1}{2}}\mathrm{\Omega }_2\sigma _{cb}{\displaystyle \frac{1}{2\mathrm{}}}ϵ_3d_3\sigma _{db}+{\displaystyle \frac{1}{2\mathrm{}}}iU\sigma _{ad}{\displaystyle \frac{i}{2}}(\mathrm{\Gamma }_b^a+\mathrm{\Gamma }_c^a+\mathrm{\Gamma }_d^a)\sigma _{ab},`$ (5)
$`i\dot{\sigma }_{ac}`$ $`=`$ $`{\displaystyle \frac{1}{2\mathrm{}}}ϵ_1d_1\sigma _{bc}+{\displaystyle \frac{1}{2}}\mathrm{\Omega }_2(\sigma _{cc}\sigma _{aa}){\displaystyle \frac{1}{2\mathrm{}}}ϵ_3d_3\sigma _{dc}{\displaystyle \frac{i}{2}}(\mathrm{\Gamma }_b^a+\mathrm{\Gamma }_c^a+\mathrm{\Gamma }_d^a)\sigma _{ac},`$ (6)
$`i\dot{\sigma }_{ad}`$ $`=`$ $`{\displaystyle \frac{1}{2\mathrm{}}}ϵ_1d_1\sigma _{bd}+{\displaystyle \frac{1}{2}}\mathrm{\Omega }_2\sigma _{cd}{\displaystyle \frac{1}{2\mathrm{}}}ϵ_3d_3(\sigma _{dd}\sigma _{aa}){\displaystyle \frac{1}{2\mathrm{}}}iU\sigma _{ab}{\displaystyle \frac{i}{2}}(\mathrm{\Gamma }_b^a+\mathrm{\Gamma }_c^a+\mathrm{\Gamma }_d^a)\sigma _{ad},`$ (7)
$`i\dot{\sigma }_{bc}`$ $`=`$ $`{\displaystyle \frac{1}{2\mathrm{}}}ϵ_1d_1\sigma _{ac}{\displaystyle \frac{1}{2}}\mathrm{\Omega }_2\sigma _{ba}+{\displaystyle \frac{1}{2\mathrm{}}}iU\sigma _{dc},`$ (8)
$`i\dot{\sigma }_{bd}`$ $`=`$ $`{\displaystyle \frac{1}{2\mathrm{}}}ϵ_1d_1\sigma _{ad}+{\displaystyle \frac{1}{2\mathrm{}}}ϵ_3d_3\sigma _{ba}+{\displaystyle \frac{1}{2\mathrm{}}}iU(\sigma _{dd}\sigma _{bb}),`$ (9)
$`i\dot{\sigma }_{cd}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{\Omega }_2\sigma _{ad}+{\displaystyle \frac{1}{2\mathrm{}}}ϵ_3d_3\sigma _{ca}{\displaystyle \frac{1}{2}}iU\sigma _{cb},`$ (10)
where $`\mathrm{\Gamma }_b^a`$ is the decay rate of the state $`|a>`$ to $`|b>`$, etc., and $`\mathrm{\Omega }_2=ϵ_2d_2/\mathrm{}`$ is the Rabi frequency corresponding to the driving field 2.

The corresponding equations for the case (a) are obtained from the above set (1) by interchanging the indices $`bc`$, $`12`$, by setting $`ϵ_3=0`$, $`d_3=0`$, $`\mathrm{\Gamma }_d^a=0`$, by replacing the effective coupling $`iU`$ by $`=ϵ_4d_4`$ (when multiplied by $`\sigma _{dc}`$, $`\sigma _{da}`$, $`\sigma _{ac}`$, $`\sigma _{db}`$, $`\sigma _{bc}`$, $`\sigma _{cc}`$ and $`\sigma _{dd}`$) and by changing the sign of $`\omega _4`$.

The propagation equations for the signal field 1 (in the case (a)) and for both signal fields (1) and (3) (in the case (b)), in the slowly varying envelope approximation and in the conditions of the resonance read

$`{\displaystyle \frac{ϵ_1}{z}}+{\displaystyle \frac{1}{c}}{\displaystyle \frac{ϵ_1}{t}}=iNd_1{\displaystyle \frac{\omega _1}{ϵ_0c}}\sigma _{ba},`$ (11)
$`{\displaystyle \frac{ϵ_3}{z}}+{\displaystyle \frac{1}{c}}{\displaystyle \frac{ϵ_3}{t}}=iNd_3{\displaystyle \frac{\omega _3}{ϵ_0c}}\sigma _{da},`$ (12)
where $`N`$ is the atom density and $`ϵ_0`$ is the vacuum electric permittivity. Similarly as in earlier papers, we have neglected propagation effects for the driving fields, i.e. we take $`ϵ_{2,4}=ϵ_{2,4}(t)`$.

Eqs (1) and (2) have been solved numerically in the moving window frame of reference: $`t^{}=tz/c,z^{}=z`$. Switching the driving field 2 on and/or off was modeled by a hyperbolic tangent, while the additional pulse 4 was taken rectangular. The initial probe pulse was taken as the sine square shape
$`ϵ_1(0,t)=ϵ_{10}\mathrm{sin}^2[\pi (t\tau _1)/(\tau _2\tau _1)]\mathrm{\Theta }(t\tau _1)\mathrm{\Theta }(\tau _2t)`$, while the initial condition for the atomic part was $`\sigma _{bb}(z,0)=1`$, with other matrix elements equal to zero.

We have performed model computations for data being of realistic orders of magnitude, however without making attempt to imitate any real atom. The atomic energies were $`E_a`$= -0.10 a.u., $`E_b=`$-0.20 a.u. $`E_c=`$-0.18 a.u. with $`E_d=`$ -0.22 a.u (in the case (a)) and $`E_d=E_b+10^7`$ a.u. in the case (b) (the latter value is of order of a magnetic energy splitting). The relaxation rates for the spontaneous emission from the level $`E_a`$ to $`E_b`$, $`E_c`$ and in the case (b) also to $`E_d`$ were taken equal to $`2.4\times 10^9`$ a.u., from which the dipole moments have been calculated. The dipole moment for the electric transition $`cd`$ was taken $`2.74\times 10^1`$ a.u., which corresponded to a negligible width of the level $`E_c`$. The length of the atomic sample was $`2.5\times 10^7`$ a.u. (1.3 mm) in the case (a) and $`3\times 10^7`$ a.u. (1.6 mm) in the case (b) and its density $`3\times 10^{13}`$ a.u. $`(2\times 10^{12}`$ cm-3). The initial signal pulse length was $`10^{11}`$ a.u. (2.4 $`\mu `$s) and $`ϵ_{10}=10^{10}`$ a.u. (which corresponded to the power density of 3.5$`\times 10^4`$ Wcm-2); the maximum value of the amplitude of the control field 2 was $`1.2\times 10^9`$ a.u. (50 mWcm-2). The values of the $`bd`$ effective coupling $`U`$ in the case (b) were of order of $`10^{10}`$ a.u. while in the case (a) we took a coupling with $`ϵ_4=2\times 10^9`$ a.u.

In Fig.2 we show the released part of the pulse 1 as a function of the local time $`t^{}`$ for different values of the area of the pulse 4 (case(a)). The pulse can be lowered, completely damped or its sign reversed depending on the final phase of the Rabi oscillations between the levels c and d. Of course the final results do not depend on particular time instants of switching the interaction 4 on and off, provided that the pulse arrived after the signal pulse 1 has been stored and before the release stage has started. The presence of the Rabi oscillations becomes clearly visible in the situation in which the control pulses 2 and 4 partially overlap. In this case, with the latter pulse being now by an order of magnitude stronger than before, the restored pulse is constructed of parts freed in those intervals of the Rabi period in which the coherence $`\sigma _{bc}`$ differs significantly from zero. The Rabi oscillations between the levels c and d are thus imposed on the leaving signal pulse (see Fig.3).

The Rabi oscillations due to the additional control field may be used not only to destroy in a reversible way the atomic coherence $`\sigma _{ab}`$ necessary to release the pulse 1 (case (a)) but also to create a new coherence $`\sigma _{ad}`$ which can be converted into a new pulse 3 (case (b)). In Fig.4 we show the shapes of the two signal pulses 1 and 3 for different values of the area of the control pulse 4, switched on and off in the storage stage. If the area is a multiple of $`\pi `$ only the pulse 1 is released, with its sign being changed in the case of an odd multiple. If the area is an odd multiple of $`\frac{\pi }{2}`$ only the pulse 3 appears, alternatively with a changed sign. For pulse areas being not a particular multiple of $`\frac{\pi }{2}`$ both pulses 1 and 3 are released, their heights being under control.

As in the previous papers the problem can be analyzed in terms of dark state polaritons. Such an analysis allows one to describe the whole process of light storing in a single $`\mathrm{\Lambda }`$ system in terms of a shape preserving solution of the Maxwell-Bloch equations (Eqs (1,2)), the components of which, i.e. the signal field and the atomic coherence, adiabatically turn one into another. In the case of a four-level system the evolution could not in general be fully adiabatic, which means that bright-state polaritons must appear at some stage of the process and are later damped . Thus the dark-state polaritons at the initial and final stages are not identical.

The approach of Ref.  generalized in our case (b) leads to the following results. One can attempt to solve Eqs (1) and (2) perturbatively (as concerns signal fields), in an adiabatic and relaxationless approximation. The stage of light stopping occurs as in the case of a three-level system: the polariton solution

$$\mathrm{\Psi }=\frac{\mathrm{\Omega }_2ϵ_1+\frac{2\omega _1Nd_1}{ϵ_0}\sigma _{bc}}{\sqrt{\mathrm{\Omega }_2^2+\frac{2\omega _1Nd_1^2}{\mathrm{}ϵ_0}}}[\frac{d_2}{|d_2|}]$$
(13)
describes an adiabatic conversion of the pulse 1 into the coherence $`\sigma _{bc}`$ and the sign correction guarantees that for large $`|\mathrm{\Omega }_2|`$ we get $`\psi =ϵ_1`$ (we assume that $`ϵ_2>0`$). After pulse stopping, say at time instant $`t_1`$, the control pulse 4 is switched on and is present up to the time instant $`t_2`$. As a consequence the density matrix evolves and at $`t=t_2`$ we obtain in the case (b) $`\sigma _{bb}=\mathrm{cos}^2\theta `$, $`\sigma _{dd}=\mathrm{sin}^2\theta `$, $`\sigma _{bd}=\mathrm{sin}\theta \mathrm{cos}\theta `$, $`\sigma _{bc}=\sigma _{bc}(t_1)\mathrm{cos}\theta `$, $`\sigma _{dc}=\sigma _{bc}(t_1)\mathrm{sin}\theta `$, where $`\theta =\frac{U(t_2t_1)}{2\mathrm{}}`$ is the pulse area.

At time instant $`t_3`$ ($`t_3>t_2`$) the control field 2 is switched on in order to release the trapped pulse. In the assumed approximations $`\sigma _{bb}`$, $`\sigma _{dd}`$ and $`\sigma _{bd}`$ do not change any more. The pulses 1 and 3 satisfy the equations

$$(\frac{}{t}+c\frac{}{z})ϵ_j=\frac{1}{\mathrm{\Omega }_2}\frac{}{t}\frac{1}{\mathrm{\Omega }_2}\underset{k}{}M_{jk}ϵ_k,$$
(14)
where $`j,k=1,3`$ and $`M_{11}=\frac{2N\omega _1d_1^2}{ϵ_0\mathrm{}}\mathrm{cos}^2\theta `$, $`M_{13}=\frac{2N\omega _1d_1d_3}{ϵ_0\mathrm{}}\mathrm{sin}\theta \mathrm{cos}\theta `$, $`M_{31}=\frac{2N\omega _3d_3d_1}{ϵ_0\mathrm{}}\mathrm{sin}\theta \mathrm{cos}\theta `$, $`M_{33}=\frac{2N\omega _3d_3^2}{ϵ_0\mathrm{}}\mathrm{sin}^2\theta `$.

Eqs (4) can be decoupled by a linear transformation. One of the solutions can be shown, similarly as in previous papers, to be a shape-preserving solution traveling with a time-dependent velocity

$$v(t)=c\frac{1}{1+\frac{2N(d_1^2\omega _1\mathrm{cos}^2\theta +d_3^2\omega _3\mathrm{sin}^2\theta )}{\mathrm{}ϵ_0\mathrm{\Omega }_2^2}}.$$
(15)
(The other solution is zero due to the initial conditions at $`t=t_2`$.) The polariton, being a combination of two fields and two coherences, has the form

$`\mathrm{\Psi }=\sqrt{d_1^2\omega _1\mathrm{cos}^2\theta +d_3^2\omega _3\mathrm{sin}^2\theta }{\displaystyle \frac{\sqrt{\omega _1})}{\sqrt{1+\frac{2N}{ϵ_0\mathrm{}\mathrm{\Omega }_2^2}(d_1^2\omega _1\mathrm{cos}^2\theta +d_3^2\omega _3\mathrm{sin}^2\theta )}}}{\displaystyle \frac{d_1}{|d_1|}}`$ (16)
$`\times [{\displaystyle \frac{d_1ϵ_1\mathrm{cos}\theta d_3ϵ_3\mathrm{sin}\theta }{d_1^2\omega _1\mathrm{cos}^2\theta +d_3^2\omega _3\mathrm{sin}^2\theta }}+{\displaystyle \frac{2N}{ϵ_0\mathrm{\Omega }_2}}(\sigma _{bc}\mathrm{cos}\theta \sigma _{dc}\mathrm{sin}\theta )].`$ (17)
The solution (6) has been normalized so that it is equal to the solution (3) at $`t=t_2`$, i.e. $`\mathrm{\Psi }(t_2)=\sqrt{\frac{2N\mathrm{}\omega _1}{ϵ_0}}[\sigma _{bc}(t_2)\mathrm{cos}\theta \sigma _{dc}(t_2)\mathrm{sin}\theta ]`$. However, the final form of the polariton is, again for large $`ϵ_2`$, a combination of the signal fields with the coefficients different from $`\mathrm{cos}\theta `$ and $`\mathrm{sin}\theta `$, except in the case of $`d_1=d_3`$ and $`\omega _1=\omega _3`$. In particular for $`\theta =\frac{\pi }{2}`$ one finds that $`\mathrm{\Psi }\sqrt{\frac{\omega _1}{\omega 3}}ϵ_3\frac{d_1d_3}{|d_1d_3|}`$ instead of $`ϵ_3\frac{d_1d_3}{|d_1d_3|}`$. As described in detail in our previous paper , this means that the evolution cannot be fully adiabatic and bright state polaritons (which are later damped) must be invoked. In the case (a) the nonadiabatic element of the evolution is even more conspicuous: only the part of the transformed coherence, namely that proportional to $`\mathrm{cos}\theta `$ is active in the release stage and turns adiabatically, being a shape-preserving solution, into $`\mathrm{cos}\theta ϵ_1`$. The other part of the excitation, namely that proportional to $`\mathrm{sin}\theta `$, ”survives” the release stage inside the medium unless relaxations in the states $`b`$ and $`d`$ are taken into account.

In summary, we have demonstrated that a modification of the atomic coherence due to a stopped light pulse can be used as a new way of changing the properties of the released light. One can in particular release two pulses of different frequencies or polarizations, with their envelopes being regulated in a continuous way. This may serve as a kind of a switch which allows one e.g., to continuously steer the information by sending it to particular channels or to temporarily hide it.

###### Acknowledgements.

This work is a part of a program of the National Laboratory of AMO Physics in Toruń, Poland