text:
# Scattering phase shifts in quasi-one-dimension

## I Introduction

Elastic scattering in one, two and three dimensions is well understood . In a scattering process there are some important physical quantities like scattering amplitude and scattering phase shift. While the scattering intensity is directly related to the scattering amplitude, the scattering phase shifts are also very important physical quantities and the Friedel-sum-rule (FSR) relates them to the density of states (DOS).

At low temperatures, inelastic collisions are greatly suppressed. As a result the phase coherence length of an electron can become a few microns. Mesoscopic systems are defined as systems in which the phase coherence length exceeds the sample size. In such a system elastic scattering is the dominant feature and such mesoscopic samples can be understood as phase coherent elastic scatterers. With the experimental realizations of mesoscopic systems and their possibility of being applied in nano-technology and quantum computing, understanding scattering effects in quasi-one dimensions (Q1D) has also become important at present. This is also essential because the Landauer conductance formula relates the conductance to partial scattering intensities and also one can now probe scattering phase shift directly in an experiment . Mesoscopic samples are normally made up of metals or semiconductors, and the defects in them are generally point defects. Hence we will restrict our analysis to delta function potential impurities.

Recently a new kind of scattering phase shift was discussed in Q1D, in connection with the violation of the parity effect . To explain this phase shift and the violation of parity effect, here we elaborate some portions of Ref.  and explain what are symmetry dictated nodes (SDN) and non-symmetry dictated nodes (NSDN) that can arise in a Q1D system. For example let us consider a rectangular quantum biliard or dot connected to leads by quantum mechanical tunneling as shown in Fig.1. The system has reflection symmetry across the x-axis as well as the y-axis. Also the x-y components separate and spanning nodes (nodes that span across the direction of propagation and shown by dashed line in Fig.1) as well as non-spanning nodes (shown by dotted line in Fig.1) develop in the geometry, dictated by the reflection symmetries. There are various symmetries that give rise to nodes in the wave-function and we call them SDN (for example the antisymmetric property of the many body wave function result in nodes). But if quantum coherence extends to some distance inside the leads then one can model the phase coherent quantum dot as shown in Fig.2.

Note that in this case also reflection-symmetry holds in the x-direction as well as in the y-direction, but x-y components do not separate. So by tuning the boundary condition in y-direction by a gate voltage one can develop nodes that try to develop across y-direction but also act across x-direction and change the phase of the wave function in x-direction by $`\pi `$. We call them NSDN because they do not originate due to the symmetry of the Hamiltonian. There are many configurations of such NSDN , the simplest one was discussed for the stub geometry (shown in Fig.3) in ref. .

Consider for example two finite one dimensional quantum wires of equal length, AB and CD placed perpendicular to each other as shown in Fig.3 by the solid lines. When CD is completely detached from AB as shown in Fig.3(a), then the quantum mechanical wave function in AB and CD in the ground state is shown by the dotted lines. They are basically the ground state wave function in an infinite potential well in one dimension (1D). As is known to us, the ground state wave functions are by symmetry, even parity states without any nodes, except at the boundary. But when CD is attached to AB to give a T-shaped stub structure as shown in Fig.3(b), then CD forms a node at C, which is also the midpoint of AB. The wave function in this case is again shown by dotted lines and the wave function between A and B is no longer an even parity state but an odd parity state. The node at C between A and B does not originate from the symmetry of the Hamiltonian and is not a symmetry dictated node (SDN). It is rather forced by the boundary condition in the y-direction and is a NSDN. An infinitesimal change in the length CD makes this node disappear and then we have no node between A and B. The node at C induces a phase change by $`\pi `$ and when we join A and B together to form a ring-stub system, we also get persistent currents without parity effect, as parity of the persistent currents is sensitive to the number of nodes in the wave function .

Fano resonances  are a very general feature of Q1D  systems in the presence of defects and the Fano resonances are characterized by a zero-pole pair of the transmission amplitude in the complex energy plane. When semi infinite leads are attached to A and B in Fig.3(b) then the transmission amplitude of the stub structure also has zero-pole pair and one gets Fano resonances . At the energy corresponding to the pole a charge gets trapped by the scatterer and there is also an energy when there is a zero transmission across the scatterer. At the energy corresponding to the zero, scattering phase shift discontinuously changes by $`\pi `$ due to the NSDN .

It seems at present that this phase due to NSDN is necessary to understand the experimental results of Refs. . Initial analysis of the experiments in terms of Friedel-sum-rule revealed the shortcomings of applying Friedel-sum-rule to the quantum dot. The phase change due to NSDN can explain the experimental data was first proposed in Ref. and later discussed in Refs..

The Friedel-sum-rule (FSR) can be stated as

$$\theta (E_2)\theta (E_1)\pi N(E_2,E_1).$$
(1)
In 1D, 2D and 3D, the equality is known to be approximate, and is almost exact in the WKB regime where generally transport occurs. Here $`N(E_2,E_1)`$ is the variation in the number of states in the energy interval $`[E_1,E_2]`$ due to the scatterer and

$$\theta =\frac{1}{2}\xi _i=\frac{1}{2i}ln(det[S]),$$
(2)
$`S`$ being an n$`\times `$n scattering matrix and $`e^{i\xi _i}`$, i=1,2,….n are the n eigenvalues of the unitary matrix $`S`$. In differential form the FSR can also be stated as

$$\frac{\theta }{E}=\frac{1}{2i}\frac{}{E}ln(det[S])\pi (\rho (E)\rho _0(E)),$$
(3)
where $`(\rho (E)\rho _0(E))`$ is the variation of the DOS or the difference in the DOS due to the presence of the potential. $`\rho (E)`$ and $`\rho _0(E)`$ can be found by integrating the local DOS $`\rho (x,y,z,E)`$ and $`\rho _0(x,y,z,E))`$ which are related to the electron probability. However, if FSR is to be useful in mesoscopic systems where scattering phase shifts can be accurately measured, and where we are always interested to study a finite region of space, then we must see to what extent it can give the local DOS. In fact in 1D, 2D or 3D it does so extremely well and so we ask the question that how good it is in Q1D. We shall show that in Q1D, FSR will fail to give the local DOS as well as the global DOS at all energy regimes.

For a symmetric scatterer in a strictly 1D system when transmission and reflection amplitudes are denoted by $`t`$ and $`r`$, respectively,

$`S=\left(\begin{array}{cc}r& t\\ t& r\end{array}\right).`$In this case, one can show

$$\frac{\theta }{E}=\frac{arg(t)}{E},$$
(4)
which means sum of the phases of the eigen values of S is equal to the phase of some particular matrix element of S.

But now we know that in systems that are not strictly 1D, one can have zero-pole pairs and then Eq.(4) is not valid because of the $`\pi `$ phase shifts induced by NSDN. Note that for the system in Fig.3(b), although the scattering matrix is $`2\times 2`$, the point C is connected to 3 directions and is hence a Q1D system. A quantum wire with a finite width and only one propagating channel  also is a Q1D system with a $`2\times 2`$ scattering matrix. For such systems, that are not strictly one dimensional but has a $`2\times 2`$ scattering matrix

$`{\displaystyle \frac{\theta }{E}}{\displaystyle \frac{arg(t)}{E}}`$ (5)
$`\text{but}{\displaystyle \frac{\theta }{E}}\pi (\rho (E)\rho _0(E)),`$ (6)
i.e., when we go from 1D (with 2x2 S matrix) to Q1D (also with 2x2 S matrix) then Eq. (3) holds but Eq. (4) does not hold. This analysis was presented by Lee  and by Büttiker and Taniguchi. Their analysis is restricted to the system in Fig.3(b) and $`S`$-matrices that are 2$`\times `$2. Eq.(3) is not violated in the presence of NSDN and $`\pi `$ phase slips because $`det[S]=r^2t^2`$ and if $`arg(t)`$ changes by $`\pi `$ then $`[arg(t^2)]`$ changes by $`2\pi `$ or $`0`$ and hence $`det[S]`$ is unaffected by the $`\pi `$ phase slips. Thus for such a single propagating channel this phase shift is well understood by now. It has been emphasized that the multi-channel case also needs to be studied , specially since scattering phase shifts can now be probed experimentally in the single channel case  as well as in the multi-channel case , but no such study has been reported so far.

In particular, the FSR is very important in condensed matter Physics, because from the scattering phase shift (that can be determined experimentally) one can know the local DOS inside a disordered sample without knowing its internal details. Even theoretically, except in very simple situations, the wave function inside the scatterer has infinite degrees of freedom and exact calculation of local DOS from the exact wave function inside the scatterer may be non-trivial. While a good estimate of local DOS from the S-matrix, which is on its own a very useful quantity, can greatly reduce the complexities. In this work we will study the $`n`$ channel scattering problem in a Q1D quantum wire, with special emphasis on FSR and Wigner delay time . Refs.  and  parameterize the $`S`$-matrix in a particular way (there are in fact many different ways of parameterizing the $`S`$-matrix ) in which the scattering matrix elements become independent of energy. We will show that this energy dependence, that are not important in 1D play a very crucial role in Q1D multi-channel scattering. Hence in section II we will generalize the work of Refs. and  for real energy dependent 2$`\times `$2 scattering matrices. The $`n`$ channel case will be analyzed in section III and IV. In section V we will show some novel phase shifts at critical energies where $`S`$ matrix changes dimensions. Section VI is devoted to conclusions.

## II Scattering in one-dimension and negative values of $`d\theta /dE`$

In Fig.4 we consider a potential that is described in details in the figure caption. The quantum mechanical wave function or the solution to the Schrödinger equation in different regions is also shown and explained in the figure and its caption. We will always normalize the incoming wave-function such that its amplitude is 1. Griffiths boundary conditions for this system gives the following equations  (we use $`2m=1`$ and $`\mathrm{}=1`$).

$$1+r=a+b,$$
(7)
$$ae^{ikl}+be^{ikl}=t,$$
(8)
$$ik(1r)ik(ab)=V(1+r)\text{and}$$
(9)
$$ik(ae^{ikl}be^{ikl})ikt=V(ae^{ikl}+be^{ikl}).$$
(10)
We will first analyze this system in detail and generalize the results of Refs.  further by considering realistic energy dependent $`r`$ and $`t`$, that will later help us to accentuate the new features that can be observed in a multi-channel disordered quantum wire.

First of all let us calculate the local DOS and global DOS to see how much it agrees with $`d\theta /dE`$. Although the basic facts discussed in this section is known in the Greens function formalism, to the best of our knowledge, quantitative disagreement (or agreement) has not been shown so far. Using quantum mechanical expression for the local DOS integrated over the region II in Fig. 4, i.e.,

$$\rho _R=\frac{2}{hv}_0^lae^{ikx}+be^{ikx}^2𝑑x,$$

it is easy to show from Eq. 3 that (the Eq. below is consistent with Ref.  and some hints on its derivation is given in Ref )

$$\frac{d\theta }{d(kl)}\overline{\rho }+\frac{\rho _q}{l}=\rho ^{}\text{(say)},$$
(11)
$$\text{where}\overline{\rho }=a^2+b^2\text{and}$$
(12)
$$\rho _q=_0^l\left(ab^{}e^{2ikx}+ba^{}e^{2ikx}\right)𝑑x.$$
(13)
Here $`\frac{\rho _q}{l}`$ is a term that arises because of quantum mechanical interference and it can be seen that the integrand in Eq.(13) oscillates with $`x`$. For $`|\frac{E}{V}|>1`$ (this is the WKB regime when the electron does not feel the potential strongly and is almost entirely transmitted) $`\frac{\rho _q}{l}`$ is negligibly small. This is shown in Fig.5, where we plot $`\rho ^{^{}}`$ (the dashed curve) and $`\overline{\rho }`$ (the dotted curve). The two curves are almost the same for $`|\frac{E}{V}|>1`$, which means $`\frac{\rho _q}{l}`$, being the difference between the dashed and dotted curves is vanishingly small above this energy. It is known that to get the equality between the LHS and RHS of Eq. 11, it is necessary to drop the term $`\frac{\rho _q}{l}`$ . It is also known that this deviation arises because we are considering the local DOS rather than the global DOS. If we consider the global DOS, i.e.,

$$\rho (E)=\frac{2}{hv}_{\mathrm{}}^{\mathrm{}}\psi ^{}(x)\psi (x)𝑑x$$

where $`\psi (x)`$ is the quantum mechanical wavefunction at x, then instead of Eq. (11) we get

$$\frac{d\theta }{d(kl)}Lt_l\mathrm{}(|a|^2+|b|^2).$$

The equality is still approximate although almost exact in the WKB regime because the RHS is positive definite, while it is well known that the LHS can become negative at low energy (non-WKB regime) as will be demonstrated below. To get the equality it is necessary to neglect something else as we shall soon see. This second thing that we shall drop is however not due to the fact that we are considering the local DOS. It is an inherent approximation of the FSR even when we are considering the global DOS and hence also when we are considering the local DOS. The $`\rho _q/l`$ term or the interference term inside the scatterer does not arise in the case when $`l0`$ as in the case to be considered in section IV. All the deviation to be observed there is due to this inherent weakness of the FSR. We shall also see below that this inherent weakness of the FSR is negligible in 1D, 2D or 3D but becomes very formidable in quasi 1D.

One can prove that

$`\overline{\rho }=a^2+b^2={\displaystyle \frac{1r^{^{}}^4}{1r_{}^{^{}}{}_{}{}^{2}\tau ^2^2}},`$for any energy dependent reflection amplitude $`r^{^{}}`$ of one of the two identical scatterers in Fig.4, where, $`\tau =e^{ikl}`$. Hence as indicated by Eq. 11, it would be interesting if we can obtain a good estimate of $`\overline{\rho }`$ or $`\rho ^{}`$ from $`\theta `$. In the appendix I it is explained that if $`\frac{dr^{}}{dE}=\frac{dt^{}}{dE}0`$ (which means the scatterers are non-dispersive and that only happens at high energy in 1D, 2D and 3D) then, $`\frac{d\theta }{d(kl)}`$ reduces to the expression $`\frac{1r^{^{}}^4}{1r_{}^{^{}}{}_{}{}^{2}\tau ^2^2}`$, and then therefore, $`\frac{d\theta }{d(kl)}=\overline{\rho }`$. It is shown in section III of Ref.  (see Eq. 6 and 7 therein), that to relate $`\frac{d\theta }{dE}`$ to the global DOS, one has to neglect the energy dependence of the self energy, that depends on the coupling of the system to the leads, i.e., $`r^{}`$ and $`t^{}`$. Thus our results are consistent with that. The exact $`\frac{d\theta }{d(kl)}`$ is shown in Fig.5 by the solid curve. Note that in the relevant energy regime ($`|\frac{E}{V}|>1`$), the solid curve is very close to the dashed and dotted curves, which means FSR works very well for the local DOS as well as for the global DOS. But for $`|\frac{E}{V}|<1`$, $`\frac{d\theta }{d(kl)}`$ deviates from $`\overline{\rho }`$. But since transport effects in weak localization or diffusive or ballistic regime occur at Fermi energies, that is normally higher in semiconductors as well as metals in comparison to the energy where the two curves deviate substantially from each other, Friedel sum rule is often useful in condensed matter to obtain a good estimate of local DOS as well as global DOS.

$`d\theta /dE`$ is also well known as Wigner delay time . In the stationary phase approximation, it gives the time spent by the scattered particle at the impurity site. In the low energy regime, where dispersion becomes significant, the stationary phase approximation is not valid and $`d\theta /dE`$ can become negative and does not give a meaningful particle delay time. In this regime $`d\theta /dE`$ becomes negative as the phase velocity becomes larger than the group velocity and even larger than the velocity of light, and although such super-luminous particles can be detected experimentally they cannot carry any signal or information. In Fig.6 we show the negative behavior of $`d\theta /d(kl)`$. We find that as the strength of the impurities is varied, $`d\theta /d(kl)`$ can become more or less negative (see Fig.6), maximizing at $`Vl^2=2.1`$ for the symmetric delta potentials. The energy regime, where $`d\theta /d(kl)`$ can be negative remains the same for all $`V`$ and always $`|\frac{E}{V}|<1`$. We have checked for all these values of $`V`$ that apart from this insignificant energy range, FSR works very well. FSR has a close counterpart in quantum mechanics called Levinson’s theorem. It is known that Levinson’s theorem also breaks down in the presence of zero energy bound states  that can be degenerate with scattering states.

## III  Wigner delay time in quasi-one-dimensions

In Fig.7 we consider a quasi-one-dimensional quantum wire with an attractive impurity at (0, $`y_i`$), having electrons confined along the y-direction but free to move along the x-direction. While the states far away from the impurity are good momentum states, the impurity can mix the different modes and in this region of mode mixing, the wave function is $`\psi (x,y)=\mathrm{\Sigma }_nc_n(x)\chi _n(y)`$, where $`\chi _n(y)`$ are the transverse wave functions in the absence of the impurity and $`c_n(x)`$ are position dependent coefficients that has to be determined by mode matching. The confining potential in the y-direction or the transverse direction is taken to be hard wall. Thus the transverse wave-function is of the form $`\chi _n(y)=Sin\frac{n\pi }{W}(y+\frac{W}{2})`$. For a given width $`W`$ of the quantum wire one can choose the energy range of the incident electron such that only two modes are propagating, although, all the other modes (infinite in number, showing that the internal wave function can have infinite degrees of freedom, which makes it difficult to calculate the exact local DOS from the internal wave function) will be present but as evanescent modes. For example, if the energy of the electron be $`E`$ then for propagation in the n-th transverse mode (in short we will refer this as n-th mode) the wave-function is of the form

$$Sin\frac{n\pi }{W}(y+\frac{W}{2})e^{ik_nx}$$
(14)
where $`k_n=\sqrt{EE_n}`$, $`E_n`$ being $`\frac{n^2\pi ^2}{W^2}`$ and $`n=1,2,3,\mathrm{}`$ $`\mathrm{}`$. Here we have used $`\mathrm{}=2m=1`$. To have the n-th mode to be propagating it is necessary that $`k_n^2>0`$ or

$$n<\frac{W}{\pi }\sqrt{E}.$$
(15)
Thus we can choose the energy range where there will be two propagating modes, i.e., $`n=1`$ and $`n=2`$ satisfy condition (15). The rest of the modes ($`n>\frac{W}{\pi }\sqrt{E}`$) will be evanescent, whose wave functions are of the form

$$Sin\frac{n\pi }{W}(y+\frac{W}{2})e^{\kappa _nx},$$
(16)
where $`\kappa _n=\sqrt{E_nE}`$. These evanescent modes just renormalize the scattering matrix elements and drop out of the problem. The transmission amplitude from m-th incident mode to n-th scattered mode is given by

$`t_{mn}^{^{}}={\displaystyle \frac{i\mathrm{\Gamma }_{mn}}{2d\sqrt{k_mk_n}}},`$ (17)
$`\text{where}d=1+{\displaystyle \stackrel{e}{}}{\displaystyle \frac{\mathrm{\Gamma }_{nn}}{2\kappa _n}}+i{\displaystyle \stackrel{p}{}}{\displaystyle \frac{\mathrm{\Gamma }_{nn}}{2k_n}}.`$ (18)
Here $`^e`$ denotes the sum over all evanescent modes and $`^p`$ denotes the sum over all propagating modes. Eq.(17) holds only for inter-subband transmission amplitudes and all reflection amplitudes. When we say reflection amplitude we mean the following. For an electron incident from the left, all outgoing channels to the left are reflection channels. According to this convention $`t_{11}^{^{}}`$, $`t_{22}^{^{}}`$, $`t_{21}^{^{}}`$ and $`t_{12}^{^{}}`$ are reflection amplitudes. Inter-subband transmission amplitude are then obviously $`t_{14}^{^{}}`$ and $`t_{23}^{^{}}`$. The intra-subband transmission amplitudes $`t_{13}^{^{}}`$ and $`t_{24}^{^{}}`$ (according to numbering of channels explained in Fig.7) are given by

$$t_{13}^{^{}}=1+t_{11}^{^{}}\text{and}t_{24}^{^{}}=1+t_{22}^{^{}}.$$
(19)
Here $`\mathrm{\Gamma }_{nm}`$ is the strength of coupling between the n-th mode and the m-th mode. If we take the impurity to be a delta function potential i.e., $`V(x,y)=\gamma \delta (x)\delta (yy_i)`$, and the confining potential in the y-direction to be hard wall ($`V=\mathrm{}`$ for $`\frac{W}{2}y\frac{W}{2}`$, and $`0`$ everywhere else except the impurity site $`\times `$ ) (see Fig.7) then

$`\mathrm{\Gamma }_{nm}=\gamma Sin{\displaystyle \frac{n\pi }{W}}(y_i+{\displaystyle \frac{W}{2}})Sin{\displaystyle \frac{m\pi }{W}}(y_i+{\displaystyle \frac{W}{2}}).`$

Apart from the two propagating modes we consider two evanescent modes and truncate the infinite series of evanescent modes (note that although the series is strongly converging, the reason for truncating is different, stronger and explained in more detail after Fig. 12) in Eq. (18) and so Eq. (18) becomes

$$1+\frac{\mathrm{\Gamma }_{33}}{2\kappa _3}+\frac{\mathrm{\Gamma }_{44}}{2\kappa _4}+i\left(\frac{\mathrm{\Gamma }_{11}}{2k_1}+\frac{\mathrm{\Gamma }_{22}}{2k_2}\right)=d_2\text{(say)}.$$
(20)
The lowest evanescent mode (putting $`n`$=3 in Eq. (16) gives its wave function) has even parity in the transverse direction. For a negative impurity potential $`i.e.`$, $`\gamma <0`$, it also has a bound state at $`E=E_{3b}`$, where $`E_{3b}`$ is given by the solution of

$$1+\frac{\mathrm{\Gamma }_{33}}{2\kappa _3}+\frac{\mathrm{\Gamma }_{44}}{2\kappa _4}=0.$$
(21)
Since $`E_{3b}<\frac{9\pi ^2}{W^2}`$, $`E_{3b}`$ can be degenerate with scattering states (the $`n`$=1 and $`n`$=2 modes are the scattering states). The higher evanescent mode (putting $`n`$=4 in Eq. (16) gives its wave function) has odd parity in the transverse direction and this too has a bound state at $`E=E_{4b}`$, where $`E_{4b}`$ is given by the solution of

$$1+\frac{\mathrm{\Gamma }_{44}}{2\kappa _4}=0.$$
(22)
Once again depending on $`\gamma `$, $`E_{4b}`$ can be degenerate with the scattering states. The effect of including more evanescent modes is just to renormalize the strength of the impurity potential and does not give anything new .

The scattering matrix in this case is

$$S=\left[\begin{array}{cccc}t_{11}^{^{}}\hfill & t_{12}^{^{}}\hfill & t_{13}^{^{}}& t_{14}^{^{}}\hfill \\ t_{21}^{^{}}\hfill & t_{22}^{^{}}\hfill & t_{23}^{^{}}& t_{24}^{^{}}\hfill \\ t_{31}^{^{}}\hfill & t_{32}^{^{}}\hfill & t_{33}^{^{}}& t_{34}^{^{}}\hfill \\ t_{41}^{^{}}\hfill & t_{42}^{^{}}\hfill & t_{43}^{^{}}& t_{44}^{^{}}\hfill \end{array}\right]=\left[\begin{array}{cc}r_{2c}& t_{2c}\\ \widehat{t}_{2c}& \widehat{r}_{2c}\end{array}\right],$$
(23)
$`\text{where}r_{2c}=\left[\begin{array}{cc}t_{11}^{^{}}& t_{12}^{^{}}\\ t_{21}^{^{}}& t_{22}^{^{}}\end{array}\right]`$

$`\text{and}t_{2c}=\left[\begin{array}{cc}t_{13}^{^{}}& t_{14}^{^{}}\\ t_{23}^{^{}}& t_{24}^{^{}}\end{array}\right].`$$`\widehat{t}_{2c}=t_{2c}`$ due to time reversal symmetry and $`\widehat{r}_{2c}=r_{2c}`$ for a symmetric scatterer as that considered here. Now once again due to micro-reversebility $`t_{12}^{^{}}=t_{21}^{^{}}`$. Also $`t_{12}^{^{}}=t_{14}^{^{}}`$ because in both $`t_{12}^{^{}}`$ and $`t_{14}^{^{}}`$ the density of states in the input as well as the output channel is the same, and also the incident channel momenta and the outgoing channel momenta are the same in the transverse as well as in the propagating direction. Also $`t_{23}^{^{}}=t_{41}^{^{}}`$ because transmission amplitude should be independent of the position of the observer i.e., whether the observer is looking into the plane of the paper or out of the plane of the paper. Thus among the 16 matrix elements in Eq. (23) we are left with only 5 that are distinct. They are $`t_{11}^{^{}},t_{12}^{^{}},t_{22}^{^{}},t_{13}^{^{}}`$ and $`t_{24}^{^{}}`$.
From Eq. (17) and (19),

$`t_{11}^{^{}}={\displaystyle \frac{i\mathrm{\Gamma }_{11}}{2d_2k_1}},`$ (24)
$`t_{12}^{^{}}={\displaystyle \frac{i\mathrm{\Gamma }_{12}}{2d_2\sqrt{k_1k_2}}},`$ (25)
$`t_{22}^{^{}}={\displaystyle \frac{i\mathrm{\Gamma }_{22}}{2d_2k_2}},`$ (26)
$`t_{13}^{^{}}={\displaystyle \frac{1+\frac{\mathrm{\Gamma }_{33}}{2\kappa _3}+\frac{\mathrm{\Gamma }_{44}}{2\kappa _4}+i\frac{\mathrm{\Gamma }_{22}}{2k_2}}{d_2}}\text{and}`$ (27)
$`t_{24}^{^{}}={\displaystyle \frac{1+\frac{\mathrm{\Gamma }_{33}}{2\kappa _3}+\frac{\mathrm{\Gamma }_{44}}{2\kappa _4}+i\frac{\mathrm{\Gamma }_{11}}{2k_1}}{d_2}}.`$ (28)
Knowing these matrix elements, the scattering matrix is completely known and $`\theta `$ can also be calculated.

We find some further relationships between the scattering phase shifts as follows. First of all

$$arg(t_{11}^{^{}})=arg(t_{22}^{^{}})=tan^1\frac{Re(d)}{Im(d)}.$$
(29)
Secondly, when $`\frac{4\pi ^2}{W^2}<E_{3b}<\frac{9\pi ^2}{W^2}`$, i.e., the bound state of the 3rd subband lies in the energy range where one can have two propagating subbands, then the bound state $`E_{3b}`$ drastically changes the scattering matrix elements in that energy range. So in this energy range $`\frac{4\pi ^2}{W^2}`$ to $`\frac{9\pi ^2}{W^2}`$ we find

$$arg(t_{12}^{^{}})\frac{\pi }{2}=\theta +\pi .$$
(30)
Here negative sign is to be taken when $`E_{3b}`$ lies in this energy range. Otherwise the positive sign has to be taken. $`\theta `$ is to be calculated from Eq. (2) using Eq. (23). Thirdly we find

$$arg(t_{11}^{^{}})\pm \pi =arg(t_{12}^{^{}}).$$
(31)
Note that in contrast to Eq. (30) here the choice of $`\pm `$ sign is arbitrary. However consistent with this choice is the following

$$arg(t_{11}^{^{}})\pm \frac{\pi }{2}=\theta +\pi ,$$
(32)
where once again + sign is to be taken when $`E_{3b}`$ is present in this energy range and - sign is to be taken when absent.

We thus find very simple analytical expressions for $`\theta `$ in the sense that one need not calculate it from a $`4\times 4`$ scattering matrix but can calculate it from the argument of a single matrix element like $`t_{11}^{^{}}`$ or $`t_{12}^{^{}}`$ or $`t_{22}^{^{}}`$. These relations are analogous to Eq. (4) in section I obtained for purely one dimensional case, i.e., one need not calculate $`\theta `$ from 2$`\times `$2 matrix but one can find it from the argument of a single matrix element.

In Fig. 8, we plot only the distinct arguments of the scattering amplitudes versus energy of the incident electron. We find that all of them show negative slopes over a very large range of energy and as already discussed, such negative slopes give rise to fundamental questions in quantum mechanics . Now in Q1D we find that this negative slope is not restricted to low energy but can occur at any arbitrary energy. Notice for example, $`arg(t_{13}^{})`$ and $`arg(t_{24}^{})`$ show larger negative slopes at the highest possible energies for two channel propagation. The rest of this section will be devoted to understanding these negative slopes that at first sight looks very different in nature and character in the three curves in Fig. 8, and also to understanding what will happen when there are more than two propagating modes. We will address the FSR in the next section.

It is to be noted that among all these scattering matrix elements $`t_{11}^{^{}}`$ and $`t_{13}^{^{}}`$ exist in the single channel regime (i.e., $`\pi ^2<EW^2<4\pi ^2`$) where $`t_{11}^{^{}}`$ is the reflection amplitude and $`t_{13}^{^{}}`$ is the transmission amplitude. The phase of $`t_{13}^{}`$ in the single channel regime is known to change discontinuously by $`\pi `$ when $`t_{13}^{}`$ is 0, i.e., $`t_{13}^{}`$ has a zero in real energy. In the two channel regime if we write from simplifying Eq. (27)

$$t_{13}^{^{}}=\frac{k_2(2\kappa _3+g_3)+i\kappa _3g_2}{k_2(2\kappa _3+g_3)+i\kappa _3(g_2+\alpha g_1)},$$
(33)
where $`\alpha =\frac{k_2}{k_1}`$ and $`g_s=\frac{2\kappa _4}{\mathrm{\Gamma }_{44}+2\kappa _4}\mathrm{\Gamma }_{ss}`$; with $`s`$=1,2,3. then interestingly, we see that it has a zero in complex energy and not in real energy.

If we modify the Breit-Wigner line shape formula of 1D to include complex zeroes and write

$$t_{mbw}(E)=A\frac{EE_0+i\mathrm{\Gamma }_0}{EE_p+i\mathrm{\Gamma }_p},$$
(34)
where A is a normalization factor, then just as $`\mathrm{\Gamma }_p`$ gives the scale over which $`arg[t_{mbw}(E)]`$ increase at $`E=E_p`$, $`\mathrm{\Gamma }_0`$ gives a scale over which $`arg[t_{mbw}(E)]`$ decrease at $`E=E_0`$ where $`|t_{mbw}(E)|^2`$ also shows a minimum at $`E=E_0`$ (but not zero). One can check this very easily (let us say, when $`E_0`$=2, $`E_p`$=1 and $`\mathrm{\Gamma }_0=\mathrm{\Gamma }_p=0.5`$) and so we do not demonstrate it here. Now from Eq. (33) we see that at an energy which satisfies the condition

$$2\kappa _3+g_3=0,$$
(35)
the real part of the numerator in Eq. (33) is zero. Condition (35) is the same as the condition (21) for a bound state $`E_{3b}`$ coming from the 3rd subband that is degenerate with scattering states. So, around this energy where Eq.(35) is satisfied (lets say at $`E=E_{3b}E_0`$) $`arg(t_{13}^{})`$ will undergo a drop over an energy scale determined by the imaginary part, $`\kappa _3g_2`$, i.e., $`\mathrm{\Gamma }_0\kappa _3g_2`$.

It can be seen in Fig.9 that $`|t_{13}^{^{}}|^2`$ (dotted curve) shows a narrow minimum around an energy E$`W^284`$ (which is the solution of Eq. (35) or Eq. (21)) and at this energy $`arg(t_{13}^{^{}})`$ shows a very sharp drop over a narrow energy range determined by $`\kappa _3g_2`$. Hence by decreasing/increasing this quantity $`\kappa _3g_2`$ we can make the phase drop sharper/broader. $`g_2`$ can be made smaller in two ways, first by decreasing $`\gamma `$ and second by taking the impurity closer to a node in the transverse wave function. The plot for a decreased value of $`\gamma `$ is shown in Fig.10 and it confirms this.

Note that the quantity $`\kappa _3g_2`$ is actually energy dependent. But in Fig.9 and Fig.10 $`\kappa _3g_2`$ is so small that the drop occurs over a scale in which $`\kappa _3g_2`$ is roughly constant. For larger values of $`\kappa _3g_2`$, the phase drop will be determined by a complex competition between $`\kappa _3`$ and $`g_2`$. This is shown in Fig.11. First of all the scale of the phase drop becomes so large that any sensitivity to the position of the bound state can not be seen. Secondly, $`\kappa _3g_2`$ can not be taken to be a constant over this large scale and the enhancement of the negative slope for $`EW^2>79`$ is a signature of the fact that here $`\kappa _30`$ and so $`\kappa _3g_20`$ as $`EW^2`$ increases.

Similarly if we rewrite Eq.(28) as

$`t_{24}^{^{}}={\displaystyle \frac{k_1(2\kappa _3+g_3)+i\kappa _3g_1}{k_1(2\kappa _3+g_3)+i\kappa _3(g_1+\beta g_2)}},`$where $`\beta =\frac{k_1}{k_2}`$; then it is clear that the behavior of $`arg(t_{24}^{^{}})`$ will be qualitatively the same. It is indeed found in Fig.8 that the behavior of $`arg(t_{24}^{^{}})`$ is similar to that of $`arg(t_{13}^{^{}})`$.

$`\theta =\frac{1}{2i}ln[det[S]]`$ is shown in Fig.12 as a function of energy, for different values of $`\gamma `$. The minimum in $`\theta `$ follows the $`E_{3b}`$ and so the energy range where the slope of $`\theta `$ versus E is negative is determined by the $`E_{3b}`$. Note that when $`E_{3b}`$ goes out of this energy range the $`\theta `$ versus E has a positive slope everywhere. So in Fig.12, the negative slope arises whenever a bound state $`E_{3b}`$ is degenerate with the scattering states ($`n`$=1 and $`n`$=2), and non-monotonously scatter and disperse the scattering states. For weaker impurities in Q1D, the negative slope occur at higher energies and also are steeper as demonstrated in Fig.12. This is in contrast to what happens in 1D and demonstarted in Fig. 6, that the energy where the negative slopes occur is always for $`E/V<1`$.

We have used two evanescent modes in our calculations because one can include as many evanescent modes without changing the nature of the negative slopes as long as the positions of the bound states $`E_{3b}`$ and $`E_{4b}`$ remain the same. One can check this that with four evanescent modes and $`\gamma =6.46584`$, the negative slopes are the same as in Fig. 8, which means that the third and the fourth evanescent modes just renormalizes $`\gamma `$ from -6.46584 to -10. The exact renormalization takes place according to a formula $`\gamma ^h=\gamma ^n/d`$, where $`\gamma ^h`$ is the $`\gamma `$ value used here and $`\gamma ^n`$ is the renormalized value of $`\gamma `$ when we use $`m`$ evanescent modes instead of two. $`d=(1+\mathrm{\Gamma }_{55}^n/(2\kappa _5)+\mathrm{\Gamma }_{66}^n/(2\kappa _6)\mathrm{}.\mathrm{\Gamma }_{mm}^n/(2\kappa _m)`$, where $`\mathrm{\Gamma }_{mm}^n==\gamma ^nSin\frac{m\pi }{W}(y_i+\frac{W}{2})^2`$. Solving this one can find the renormalized value of $`\gamma `$ i.e., $`\gamma ^n`$ that keep the minimum of any of the curves for the scattering phase shifts considered here unchanged. It is worthwhile mentioning that at the band edges (i.e., $`E39`$ and 89, in the figures considered in this section), the value of any curve is independent of the number of evanescent modes, as all the modes get decoupled there. In other words, number of evanescent modes considered does not change the nature of the negative slopes. Only the positions of the bound states are important.

In order to generalize to arbitrary number of propagating channels we change our notations slightly. For electrons incident from the left/right, we call the scattered channels towards the left/right as reflection channels (the rest being transmission channels) and change our notation to

$`\stackrel{~}{r}_{11}=t_{11}^{^{}},`$ (36)
$`\stackrel{~}{r}_{22}=t_{22}^{^{}},`$ (37)
$`\stackrel{~}{r}_{12}=t_{12}^{^{}}.`$ (38)
Thus all possible reflection channels are distinct. All the intra-subband transmission channels are also distinct and they are denoted as $`\stackrel{~}{t}_{11}`$ and $`\stackrel{~}{t}_{22}`$ where

$`\stackrel{~}{t}_{11}=t_{13}^{^{}}`$ (39)
$`\text{and}\stackrel{~}{t}_{22}=t_{24}^{^{}}.`$ (40)
All other scattering matrix elements are equal to one of these 5 elements. In this notation the only difference is that the lowest channel ($`n`$=1) on the transmission side is marked 1 instead of 3. So when we say $`\stackrel{~}{t}_{11}`$ we mean transmission amplitude from the $`n`$=1 channel on the left to the $`n`$=1 channel on the right. We find from Eqs. (29), (30), (31) and (32)

$$\frac{d}{dE}arg(\stackrel{~}{r}_{mn})=\frac{d}{dE}(\frac{1}{2i}ln[det[S]])$$
(41)
We find the above relation to be true for any number of propagating modes. So $`m`$ and $`n`$ can take any integer value less than or equal to $`p`$, where $`p`$ is the total number of propagating modes. For two propagating modes $`p`$=2, for three propagating modes $`p`$=3 and so on. So Eq.(41) is analogous to the 1D case given in Eq.(4). That is when the dimension of the matrix $`S`$ becomes very large, then it is sufficient to consider the argument of a single matrix element in order to calculate the complicated quantity on the RHS of Eq. 41. In the energy regime where there are two propagating channels, the negative slopes in $`\theta `$ versus incident energy curves are determined by $`E_{3b}`$, and when there are 3 propagating channels then the negative slopes are determined by $`E_{4b}`$ and so on.

The scattering phase shifts of transmission channels i.e. $`arg(\stackrel{~}{t}_{mn})`$, where again $`m`$ and $`n`$ can take all possible integer values less than or equal to $`p`$, show sharp or gentle phase drops when the scattering states are degenerate with a bound state, depending on the value of the imaginary part in the numerator of $`\stackrel{~}{t}_{mn}`$. In the single channel regime the imaginary part in the numerator is zero and phase drops take the limiting value when the phase drops are absolutely discontinuous by $`\pi `$. Just as the discontinuous phase drop in single channel case do not affect $`\theta `$ in any way, the phase drops of the $`arg(\stackrel{~}{t}_{mn})`$ also do not affect $`\theta `$ in any way and $`\theta `$ behaves similarly as $`arg(\stackrel{~}{r}_{mn})`$.

## IV density of states and Friedel sum rule in quasi-one-dimension

The local DOS is given by the following expression

$$\rho _R=_R𝑑x_{\frac{W}{2}}^{\frac{W}{2}}𝑑y\underset{m,k_m}{}\delta (EE_{m,k_m})\psi _{m,k_m}(x,y)^2$$
(42)
Here $`E`$ is the incident energy and R is the integration region where modes are mixed. $`m`$ and $`k_m`$ are the two quantum numbers that define an incident electron wavefunction, $`A_me^{ik_mx}Sin\frac{m\pi }{W}(y+\frac{W}{2})`$ whose energy is $`E_{m,k_m}`$, where we have taken that the electron is incident from the left i.e., $`x<0`$. $`\psi _{m,k_m}(x,y)`$ is the wavefunction in the region of mode mixing and $`\psi _{m,k_m}(x,y)=_nc_n^{(m)}(x,k_n)Sin\frac{n\pi }{W}(y+\frac{W}{2})`$. Here $`c_n^{(m)}(x,k_n)=C_ne^{ik_nx}`$ for $`n=1`$ and $`n=2`$ and $`c_n^{(m)}(x,k_n)=C_ne^{\kappa _nx}`$ for $`n>2`$; $`x`$ being greater than or equal to 0. The coefficients $`C_n`$ can be determined by using the mode matching technique. The mode matching has been done in details by Bagwell . Here the delta function potential is taken to be extending from $`ϵ`$ to $`+ϵ`$ which has to be set to be tending to 0 in the end. $`\rho _{0R}`$ can be determined by replacing $`\psi _{n,k_n}(x,y)`$ by the plane wave states in absence of the scatterer and doing the integration again.

Thus we find that for any non-zero incident energy

$$(\rho \rho _0)_R=\frac{2}{hv_1}[|t_{13}|^2+|t_{14}|^2+|t_{15}|^2+\mathrm{}]$$

$$+\frac{2}{hv_2}[|t_{23}|^2+|t_{24}|^2+|t_{25}|^2+\mathrm{}].$$
(43)
Here $`v_1=\frac{\mathrm{}k_1}{m},v_2=\frac{\mathrm{}k_2}{m}`$ and $`t_{mn}=\frac{C_n}{A_m}`$. $`t_{mn}`$ can be obtained by solving the matrix Eqs. given in Ref.  (see Eq. 23 therein and we have used the same notation i.e., $`t_{mn}`$ here is the same as $`t_{mn}`$ in Eq. (23) of Ref. ). As can be seen in Fig.12 that $`\frac{d\theta }{dE}`$ is negative over a very large energy range while $`(\rho \rho _0)_R`$ as given by Eq.43 is positive.

Note that if we calculate the global DOS by taking the integration region to be from $`\mathrm{}`$ to $`\mathrm{}`$ instead of just the region R where the modes are mixed then Eq. 43 remain the same. One will get some extra integrals that are indefinite integrals but using the current conservation condition it can be analytically proved (see appendix II) that they cancel each other. This means the contribution to $`\rho (E)`$ and that to $`\rho _0(E)`$ coming from outside the region R cancel each other in the absence of a term like $`\rho _q/l`$ as in section II. Thus in this case (the proof is given in appendix II)

$$\rho (E)\rho _0(E)=(\rho \rho _0)_R$$

and both of them deviate identically from $`\frac{1}{\pi }\frac{d\theta }{dE}`$ due to strong dispersion, at any arbitrary energy.

Since the negative slopes are due to the bound states supported by the negative delta function potential, one may ask what happens for a positive delta function potential that does not support any bound state. This situation is discussed below and it also elaborates the uniqueness of the Q1D, with respect to the violation of Friedel sum rule and shows that large violation can occur in the extreme WKB limit. Fig 13 shows the energy dependence of some important scattering probabilities. As in 1D the scattering probabilities are strongly energy dependent for $`EW^2<50`$. For $`EW^2>50`$ all the curves vary slowly with energy (non-dispersive scattering). However, one scattering probability $`|\stackrel{~}{t}_{12}|^2`$ is also very strongly energy dependent at the highest energy of $`EW^289`$. ($`k_m^2/V>>1`$). In this regime scattering is almost entirely intra-channel as can be seen that $`|\stackrel{~}{t}_{11}|^2`$ and $`|\stackrel{~}{t}_{22}|^2`$ are both almost unity. All other scattering probabilities like reflection probabilities and inter-channel scattering probabilities like $`|\stackrel{~}{t}_{12}|^2`$ are extremely small. There are several other scattering matrix elements that are identical to $`\stackrel{~}{t}_{12}`$. Although being small, they can have strong energy dependence, or their energy derivative can be very large. As can be seen in Fig. 13 that the dashed curve bends down with a steep slope at $`EW^289`$. It is found in Fig. 14 that there is a strong violation of Friedel sum rule in the non-WKB regime as in 1D and also in the extreme WKB regime ($`EW^289`$ where $`k^2/V>>1`$ and scattering is almost completely intra-channel) quite unlike that in 1D. It can be seen in Fig. 13 that in the mid-energy range, the scattering probabilities are not very energy dependent and the Friedel sum rule holds in this mid-energy range as can be seen in Fig. 14. But at regimes where the scattering leads to strong dispersion, FSR is violated. As the strength of the positive delta function potential is increased, this regime where the scattering probabilities are not strongly energy dependent becomes narrower and Friedel sum rule is violated at all energies.

## V  phase behavior at critical energies

Very interesting phase behaviors can be seen at energies where the $`S`$-matrix changes dimension. For example for $`E\frac{4\pi ^2}{W^2}`$ there is only one propagating mode and the $`S`$-matrix is 2$`\times `$2. But for $`E>\frac{4\pi ^2}{W^2}`$, there are two propagating modes and the $`S`$-matrix is 4$`\times `$4. The matrix element $`t_{11}^{^{}}`$ exists on either side of the energy $`\frac{4\pi ^2}{W^2}`$ and in Fig.15 we show the behavior of $`arg(t_{11}^{^{}})`$ in the energy range that includes $`EW^2=4\pi ^2`$. Note that it exhibits a discontinuous phase drop by $`\frac{\pi }{2}`$ at $`EW^2=4\pi ^2`$. So far only discontinuous phase drops of $`\pi `$ has been observed but never $`\frac{\pi }{2}`$. From the properties of a 2$`\times `$2 $`S`$-matrix it follows that if there is a discontinuous phase change then it can only be of $`\pi `$ . So had the $`S`$-matrix been 2$`\times `$2 on either side of $`EW^2=4\pi ^2`$ the phase drop would have been $`\pi `$. But since the $`S`$-matrix is 2$`\times `$2 only on one side, including $`EW^2=4\pi ^2`$, i.e., $`E\frac{4\pi ^2}{W^2}`$, the phase drop is also one half of $`\pi `$. $`|t_{11}^{^{}}|^2`$ also has a zero at $`EW^2=4\pi ^2`$ for all possible choice of parameters , and this zero is associated with a $`\frac{\pi }{2}`$ phase jump instead of a $`\pi `$ phase jump.

Next we take a repulsive $`\delta `$ function potential. It is known  that at critical energies like $`EW^2=4\pi ^2`$, $`|t_{13}^{^{}}|^2`$ shows discontinuities. Here $`|t_{13}^{^{}}|^2`$ does not have a zero but exhibits a discontinuous jump. At these points $`arg(t_{13}^{^{}})`$ also shows non-analytic behavior as demonstrated in Fig.16. In this case $`\frac{d}{dE}arg(t_{13}^{^{}})`$ is discontinuous.

## VI Conclusions

In a multichannel quantum wire with attractive impurities, negative slopes in the scattering phase shift versus incident energy curves can occur at all possible energies. For weaker defects it happens at higher energies and the negative slopes are more pronounced. Such negative slopes mean super luminescence  that can be observed experimentally. Although such a super luminescent particle will not give any information about the particle delay or information delay, they are of interest because they demonstrate fundamental principles in quantum mechanics. Hence Eq. (41) derived in this paper may be of use to experimentalists and theoreticians. Calculations of variation of DOS due to the presence of the scatterer show that FSR can be violated at any arbitrary energy. That is violation of FSR is not restricted to non-WKB regime (the regime where anyway transport does not occur) as in 1D, 2D or 3D. Rather the violation is related to strong energy dependence of the scattering matrix elements (dispersive behavior). While in 1D, 2D and 3D the strong energy dependence of scattering occur only in the non-WKB regime, in Q1D there is no systematics. A Q1D system in the extreme WKB limit can also exhibit strongly energy dependent scattering and there is no definable regime where the FSR will work and hence it may not work even in regimes where transport occurs. For attractive impurities, weaker the impurity, stronger the violation of FSR while for repulsive impurities, stronger the impurity, stronger the violation of FSR. We also show that the discontinuous phase drops in the single channel case have a counterpart in the multichannel case wherein the drops can be continuous and we propose a line shape formula for them in Eq. 34. When there is a third channel of escape for the electron, apart from the channel along which it is incident and the channel where its scattering phase shift is measured, the phase drop becomes continuous. However, these phase drops do not affect $`\frac{1}{2i}lnDet[S]`$ and hence Friedel sum rule. Finally, we discuss some novel scattering phase shifts at energies where the $`S`$ matrix changes dimension.

## VII ACKNOWLEDGMENTS

We thank Dr. A. M. Jayannavar for useful discussions.

## VIII APPENDIX I

Considering the symmetric scattering potential in Fig.4, the scattering matrix $`S`$ of the structure can be found by cascading the scattering matrices of different parts, i.e.,

$`S=\left(\begin{array}{cc}r& t\\ t& r\end{array}\right)=S_1S_2S_3,`$

$`\text{where}S_1=S_3=\left(\begin{array}{cc}r^{}& t^{}\\ t^{}& r^{}\end{array}\right)`$

$`\text{and}S_2=\left(\begin{array}{cc}0& \tau \\ \tau & 0\end{array}\right).`$Here $`\tau =e^{i\varphi }`$, $`\varphi =kl`$ and $`k=\sqrt{\frac{2m}{\mathrm{}^2}E}`$. $`S_2`$ is the scattering matrix for the free region II of length $`l`$ between the two scatterers. $`r^{}`$ & $`t^{}`$ are the reflection & transmission amplitudes due to one of the two potentials when isolated.

After cascading these three matrices the resultant scattering matrix of our system becomes

$`S=\left(\begin{array}{cc}r^{}+{\displaystyle \frac{t^2\tau ^2r^{}}{1r^2\tau ^2}}& {\displaystyle \frac{t^2\tau }{1r^2\tau ^2}}\\ {\displaystyle \frac{t^2\tau }{1r^2\tau ^2}}& r^{}+{\displaystyle \frac{t^2\tau ^2r^{}}{1r^2\tau ^2}}\end{array}\right).`$And so,

$`det[S]=\left(r^{}+{\displaystyle \frac{t^2\tau ^2r^{}}{1r^2\tau ^2}}\right)^2\left({\displaystyle \frac{t^2\tau }{1r^2\tau ^2}}\right)^2\text{(i)}`$

$`={\displaystyle \frac{1}{(1r^2\tau ^2)^2}}(M+N),\text{(ii)}`$

$`\text{where}M=r^2(1r^2\tau ^2)^2\text{(iii)}`$

$`\text{and}N=(1r^2\tau ^2)(2t^2r^2\tau ^2t^4\tau ^2).\text{(iv)}`$From (i),

$`{\displaystyle \frac{det[S]}{\varphi }}=2\left(A{\displaystyle \frac{A}{\varphi }}B{\displaystyle \frac{B}{\varphi }}\right),\text{(v)}`$

$`\text{where}A=r^{}+{\displaystyle \frac{t^2\tau ^2r^{}}{1r^2\tau ^2}}\text{(vi)}`$

$`\text{and}B={\displaystyle \frac{t^2\tau }{1r^2\tau ^2}}.\text{(vii)}`$Using (v), (vi) and (vii) we find

$$\frac{det[S]}{\varphi }=$$

$`2\left[A\left(2r^{}B+{\displaystyle \frac{2r^3\tau ^2B}{1r^2\tau ^2}}\right)B\left({\displaystyle \frac{t^2+r^2t^2\tau ^2}{(1r^2\tau ^2)^2}}\right)\right]{\displaystyle \frac{\tau }{\varphi }}`$

$`+2\left[A\left(1+B\tau +{\displaystyle \frac{2r^2\tau ^3B}{1r^2\tau ^2}}\right)B\left({\displaystyle \frac{2r^{}\tau ^2B}{1r^2\tau ^2}}\right)\right]{\displaystyle \frac{r^{}}{\varphi }}`$

$`+2[2r^{}\tau ABB\left({\displaystyle \frac{2\tau t^{}}{1r^2\tau ^2}}\right)]{\displaystyle \frac{t^{}}{\varphi }}.\text{(viii)}`$

We will neglect the last two terms in comparison to the 1st term in Eq.(viii) but we will retain the energy dependence of $`r^{}`$ and $`t^{}`$ and the 1st term in Eq.(viii) is actually varying with energy very strongly. Hence our results correspond to real potentials and we do not parameterize the $`S`$ matrix in a special way. We apply this result to the case of double delta function potential in Fig.5 and illustrate the significance of the last two terms in comparison with the 1st one. We stress that this calculation in this appendix holds even if the $`\delta `$-function potential is replaced by the square-well or any arbitrary potential. Thus in the regime where $`\frac{r^{}}{\varphi }0`$ and $`\frac{t^{}}{\varphi }0`$ and using $`\frac{\tau }{\varphi }=i\tau `$,

$$\frac{det[S]}{\varphi }$$

$`=2\left[A\left(2r^{}B+{\displaystyle \frac{2r^3\tau ^2B}{1r^2\tau ^2}}\right)B\left({\displaystyle \frac{t^2+r^2t^2\tau ^2}{(1r^2\tau ^2)^2}}\right)\right]i\tau .`$At this point we substitute the values of $`A`$ and $`B`$ from (vi) and (vii) to get

$`{\displaystyle \frac{det[S]}{\varphi }}=2i{\displaystyle \frac{1}{(1r^2\tau ^2)^3}}N.\text{(ix)}`$

$`\text{Now},\theta ={\displaystyle \frac{1}{2i}}ln(det[S])\text{(x)}`$

$`\text{and so}{\displaystyle \frac{\theta }{E}}={\displaystyle \frac{\theta }{\varphi }}{\displaystyle \frac{\varphi }{E}}.`$From (x), (ii) and (ix)

$`{\displaystyle \frac{\theta }{\varphi }}={\displaystyle \frac{1}{2i}}{\displaystyle \frac{1}{det[S]}}{\displaystyle \frac{det[S]}{\varphi }}`$

$`={\displaystyle \frac{1}{1r^2\tau ^2}}{\displaystyle \frac{1}{\frac{M}{N}+1}}.`$Multiplying the numerator and denominator by $`\frac{(1r^2\tau ^2)^{}}{1r^{^{}}^4}`$, we get

$`{\displaystyle \frac{\theta }{\varphi }}={\displaystyle \frac{1r^{}^4}{1r^2\tau ^2^2}}{\displaystyle \frac{(1r^2\tau ^2)^{}}{1r^{}^4}}{\displaystyle \frac{1}{\frac{M}{N}+1}}`$

$`={\displaystyle \frac{1r^{}^4}{1r^2\tau ^2^2}}Q,\text{(xi)}`$

$`\text{where,}Q={\displaystyle \frac{(1r^2\tau ^2)^{}}{1r^{}^4}}{\displaystyle \frac{1}{\frac{M}{N}+1}}\text{(xii)}`$

$`=\left[{\displaystyle \frac{(1r_{}^{}{}_{}{}^{2}\tau _{}^{}{}_{}{}^{2})1+r^{}^4\tau ^4}{1r^{}^4}}+1\right]{\displaystyle \frac{1}{\frac{M}{N}+1}},`$

$`\text{as}|\tau |^4=1.`$As $`1|r^{}|^2=|t^{}|^2`$

$`Q=\left[{\displaystyle \frac{r_{}^{}{}_{}{}^{2}\tau _{}^{}{}_{}{}^{2}+r^{}^4\tau ^4}{(1+r^{}^2)|t^{}|^2}}+1\right]{\displaystyle \frac{1}{\frac{M}{N}+1}}.`$Now substituting the values of $`M`$ and $`N`$ from (iii) and (iv)

$`Q={\displaystyle \frac{\frac{r_{}^{}{}_{}{}^{2}\tau _{}^{}{}_{}{}^{2}(1r^2\tau ^2)+t^{}^2(1+r^{}^2)}{t^{}^2(1+r^{}^2)}}{\frac{r^2(1r^2\tau ^2)t^2\tau ^2(t^22r^2)}{t^2\tau ^2(t^22r^2)}}}.`$
Using, $`r^{}=r^{}e^{i\theta _r}`$ and $`t^{}=t^{}e^{i\theta _t}`$,

$$Q=[|r^{}|^2|t^{}|^2|\tau |^4e^{2i\theta _t}e^{2i(\theta _t\theta _r)}|r^{}|^4|t^{}|^2|\tau |^4\tau ^2e^{4i\theta _t}$$

$$|t^{}|^4\tau ^2e^{4i\theta _t}|t^{}|^4|r^{}|^2\tau ^2e^{4i\theta _t}$$

$$2|r^{}|^4|\tau |^4e^{2i\theta _t}+2|r^{}|^6|\tau |^4\tau ^2e^{2i(\theta _t+\theta _r)}$$

$$+2|r^{}|^2|t^{}|^2\tau ^2e^{2i(\theta _t+\theta _r)}+2|r^{}|^4||t^{}|^2\tau ^2e^{2i(\theta _t+\theta _r)}]/D,$$

$`\text{where},D=(1+r^{}^2)(t^{}^4\tau ^2e^{4i\theta _t}r^{}^4\tau ^2e^{4i\theta _r}`$

$`+2t^{}^2r^{}^2\tau ^2e^{2i(\theta _r+\theta _t)}+r^{}^2e^{2i\theta _r}).\text{(xiii)}`$
Now it follows from unitarity that $`e^{2i(\theta _t\theta _r)}=e^{i\pi }=1`$ and $`|\tau |^4=1`$ and so

$`Q=[|r^{}|^2|t^{}|^2e^{2i\theta _t}|r^{}|^4|t^{}|^2\tau ^2e^{4i\theta _t}|t^{}|^4\tau ^2e^{4i\theta _t}`$

$`|t^{}|^4|r^{}|^2\tau ^2e^{4i\theta _t}2|r^{}|^4e^{2i\theta _t}+2|r^{}|^6\tau ^2e^{2i(\theta _t+\theta _r)}`$

$`+2|r^{}|^2|t^{}|^2\tau ^2e^{2i(\theta _t+\theta _r)}+2|r^{}|^4|t^{}|^2\tau ^2e^{2i(\theta _t+\theta _r)}]/D`$

$`=[e^{2i\theta _t}|r^{}|^2\{|t^{}|^2+2|r^{}|^2\}(|t^{}|^2+|r^{}|^2\{|r^{}|^2+|t^{}|^2\})`$

$`(e^{4i\theta _t}|t^{}|^2\tau ^2e^{2i(\theta _t+\theta _r)}2|r^{}|^2\tau ^2)]/D.`$Now inside the $`\{\}`$ brackets if we use the fact that $`r^{}^2+t^{}^2=1`$, then

$`Q=[e^{2i\theta _t}|r^{}|^2(1+|r^{}|^2)e^{4i\theta _t}|t^{}|^2\tau ^2+2e^{2i(\theta _t+\theta _r)}|r^{}|^2\tau ^2]/D.`$Multiplying numerator and denominator above by $`e^{2i(\theta _r+\theta _t)}`$ and putting $`e^{2i(\theta _t\theta _r)}=e^{i\pi }=1`$, we get

$`Q={\displaystyle \frac{e^{2i\theta _r}|r^{}|^2(1+|r^{}|^2)+\tau ^2(|t^{}|^2+|r^{}|^2)+|r^{}|^2\tau ^2}{D^{}}},`$

$`\text{where,}D^{}=De^{2i(\theta _r+\theta _t)}.\text{(xiv)}`$Again using $`|r^{}|^2+|t^{}|^2=1`$,

$`Q={\displaystyle \frac{(1+|r^{}|^2)(\tau ^2|r^{}|^2e^{2i\theta _r})}{D^{}}}.\text{(xv)}`$Now from (xiv) and (xiii)

$`D^{}=(1+|r^{}|^2)[|r^{}|^2e^{2i\theta _t}|r^{}|^4\tau ^2e^{2i(\theta _t\theta _r)}`$

$`|t^{}|^4\tau ^2e^{2i(\theta _t\theta _r)}+2|t^{}|^2|r^{}|^2\tau ^2].`$As $`\theta _t\theta _r=\frac{\pi }{2}`$,

$`D^{}=(1+|r^{}|^2)[|r^{}|^2e^{2i\theta _t}+\tau ^2(|r^{}|^2+|t^{}|^2)^2]`$where of course $`r^{}^2+t^{}^2=1`$. Substituting $`D^{}`$ in Eq.(xv),

$`Q={\displaystyle \frac{\tau ^2r^{}^2e^{2i\theta _r}}{\tau ^2+r^{}^2e^{2i\theta _t}}}.`$Multiplying numerator and denominator of $`Q`$ by $`e^{2i\theta _t}`$ and using $`e^{2i(\theta _t\theta _r)}=e^{i\pi }=1`$, we get from (xi)

$`{\displaystyle \frac{\theta }{\varphi }}={\displaystyle \frac{1r^{}^4}{1r^2\tau ^2^2}}`$

## IX APPENDIX II

The global density of states is given by

$`\rho (E)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x{\displaystyle _{\frac{W}{2}}^{\frac{W}{2}}}𝑑y`$ (45)
$`{\displaystyle \underset{m,k_m}{}}\left|\psi _{m,k_m}(x,y)\right|^2\delta (EE_{m,k_m})\text{(i)}`$
where $`\psi _{m,k_m}(x,y)=_nc_n^{(m)}(x)\chi _n(y)`$ and $`E_{m,k_m}`$ is the energy of an electron in the leads.
$`E_{m,k_m}={\displaystyle \frac{m^2\pi ^2\mathrm{}^2}{2m_eW^2}}+{\displaystyle \frac{\mathrm{}^2k_m^2}{2m_e}},\text{where}m=\pm 1,\pm 2,`$ as there are two propagating modes in the leads.
As $`\chi _n(y)`$’s form an orthonormal set,

$$\rho (E)=\underset{m,k_m}{}\delta (EE_{m,k_m})_{\mathrm{}}^{\mathrm{}}𝑑x\underset{n}{}\left|c_n^{(m)}(x)\right|^2$$

$`\rho (E)`$ $`=`$ $`{\displaystyle \frac{2}{hv_1}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x{\displaystyle \underset{n}{}}\left|c_n^{(1)}(x)\right|^2`$ (46)
$`+`$ $`{\displaystyle \frac{2}{hv_2}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x{\displaystyle \underset{n}{}}\left|c_n^{(2)}(x)\right|^2\text{(ii)}`$ (47)
Here, $`v_1={\displaystyle \frac{\mathrm{}k_1}{m_e}}`$ and $`v_2={\displaystyle \frac{\mathrm{}k_2}{m_e}}`$.
Now,

$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x{\displaystyle \underset{n}{}}\left|c_n^{(1)}(x)\right|^2`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}\left|c_1^{(1)}(x)\right|^2𝑑x`$ (48)
$`+`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}\left|c_2^{(1)}(x)\right|^2𝑑x`$ (49)
$`+`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}\left|c_3^{(1)}(x)\right|^2𝑑x+\mathrm{}`$ (50)
$`=`$ $`T1\text{(say)},`$ (51)
where electron is incident in the fundamental mode,

$`c_1^{(1)}(x)`$ $`=`$ $`e^{ik_1x}+\stackrel{~}{r}_{11}e^{ik_1x}\text{for}x<0`$ (52)
$`=`$ $`\stackrel{~}{t}_{11}e^{ik_1x}\text{for}x>0`$ (53)
$`c_2^{(1)}(x)`$ $`=`$ $`\stackrel{~}{r}_{12}e^{ik_2x}\text{for}x<0`$ (54)
$`=`$ $`\stackrel{~}{t}_{12}e^{ik_2x}\text{for}x>0`$ (55)
and for $`n>2`$,

$`c_n^{(1)}(x)`$ $`=`$ $`t_{1n}e^{\kappa _nx}\text{for}x<0`$ (56)
$`=`$ $`t_{1n}e^{\kappa _nx}\text{for}x>0`$ (57)
$`T1`$ $`=`$ $`{\displaystyle _{\mathrm{}}^0}𝑑x\left[1+|\stackrel{~}{r}_{11}|^2+2|\stackrel{~}{r}_{11}|Cos(2k_1x+\eta _1)\right]`$ (58)
$`+`$ $`{\displaystyle _0^{\mathrm{}}}𝑑x|\stackrel{~}{t}_{11}|^2+{\displaystyle _{\mathrm{}}^0}𝑑x|\stackrel{~}{r}_{12}|^2+{\displaystyle _0^{\mathrm{}}}𝑑x|\stackrel{~}{t}_{12}|^2`$ (59)
$`+`$ $`|t_{13}|^2+|t_{14}|^2+\mathrm{}\mathrm{}`$ (60)
Here, $`\eta _1`$ is defined as $`\stackrel{~}{r}_{11}=|\stackrel{~}{r}_{11}|e^{i\eta _1}`$
Similarly, for electron incident in the first excited mode,

$`T2`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x{\displaystyle \underset{n}{}}\left|c_n^{(2)}(x)\right|^2`$ (61)
$`=`$ $`{\displaystyle _{\mathrm{}}^0}𝑑x\left[1+|\stackrel{~}{r}_{22}|^2+2|\stackrel{~}{r}_{22}|Cos(2k_2x+\eta _2)\right]`$ (62)
$`+`$ $`{\displaystyle _0^{\mathrm{}}}𝑑x|\stackrel{~}{t}_{22}|^2+{\displaystyle _{\mathrm{}}^0}𝑑x|\stackrel{~}{r}_{21}|^2+{\displaystyle _0^{\mathrm{}}}𝑑x|\stackrel{~}{t}_{21}|^2`$ (63)
$`+`$ $`|t_{23}|^2+|t_{24}|^2+\mathrm{}\mathrm{},`$ (64)
Here, $`\eta _2`$ is defined as $`\stackrel{~}{r}_{22}=|\stackrel{~}{r}_{22}|e^{i\eta _2}`$

Therefore, (65)
$`\rho (E)`$ $`=`$ $`2[{\displaystyle \frac{1+|\stackrel{~}{r}_{11}|^2}{hv_1}}{\displaystyle _{\mathrm{}}^0}dx+{\displaystyle \frac{1+|\stackrel{~}{r}_{22}|^2}{hv_2}}{\displaystyle _{\mathrm{}}^0}dx`$ (66)
$`+`$ $`{\displaystyle \frac{|\stackrel{~}{r}_{12}|^2}{hv_1}}{\displaystyle _{\mathrm{}}^0}𝑑x+{\displaystyle \frac{|\stackrel{~}{r}_{21}|^2}{hv_2}}{\displaystyle _{\mathrm{}}^0}𝑑x`$ (67)
$`+`$ $`{\displaystyle \frac{2|\stackrel{~}{r}_{11}|}{hv_1}}{\displaystyle _{\mathrm{}}^0}𝑑xCos(2k_1x+\eta _1)`$ (68)
$`+`$ $`{\displaystyle \frac{2|\stackrel{~}{r}_{22}|}{hv_2}}{\displaystyle _{\mathrm{}}^0}𝑑xCos(2k_2x+\eta _2)`$ (69)
$`+`$ $`{\displaystyle \frac{|\stackrel{~}{t}_{11}|^2+|\stackrel{~}{t}_{12}|^2}{hv_1}}{\displaystyle _0^{\mathrm{}}}𝑑x+{\displaystyle \frac{|\stackrel{~}{t}_{22}|^2+|\stackrel{~}{t}_{21}|^2}{hv_2}}{\displaystyle _0^{\mathrm{}}}𝑑x`$ (70)
$`+`$ $`{\displaystyle \frac{1}{hv_1}}\left(|t_{13}|^2+|t_{14}|^2+\mathrm{}\right)`$ (71)
$`+`$ $`{\displaystyle \frac{1}{hv_2}}(|t_{23}|^2+|t_{24}|^2+\mathrm{})]`$ (72)
Due to time reversal symmetry, $`\stackrel{~}{r}_{12}=\stackrel{~}{r}_{21}`$ & $`\stackrel{~}{t}_{12}=\stackrel{~}{t}_{21}`$.
Inside the square bracket $`[]`$,
in the 3rd term we put $`\stackrel{~}{r}_{12}=\stackrel{~}{r}_{21}`$,
in the 4th term we put $`\stackrel{~}{r}_{21}=\stackrel{~}{r}_{12}`$,
in the 7th term we put $`\stackrel{~}{t}_{12}=\stackrel{~}{t}_{21}`$,
in the 8th term we put $`\stackrel{~}{t}_{21}=\stackrel{~}{t}_{12}`$.

$`\text{Therefore},`$ (73)
$`\rho (E)`$ $`=`$ $`2[{\displaystyle \frac{1+|\stackrel{~}{r}_{11}|^2+|\stackrel{~}{r}_{21}|^2}{hv_1}}{\displaystyle _{\mathrm{}}^0}dx`$ (74)
$`+`$ $`{\displaystyle \frac{1+|\stackrel{~}{r}_{12}|^2+|\stackrel{~}{r}_{22}|^2}{hv_2}}{\displaystyle _{\mathrm{}}^0}𝑑x`$ (75)
$`+`$ $`{\displaystyle \frac{|\stackrel{~}{t}_{11}|^2+|\stackrel{~}{t}_{21}|^2}{hv_1}}{\displaystyle _0^{\mathrm{}}}𝑑x+{\displaystyle \frac{|\stackrel{~}{t}_{12}|^2+|\stackrel{~}{t}_{22}|^2}{hv_2}}{\displaystyle _0^{\mathrm{}}}𝑑x`$ (76)
$`+`$ $`{\displaystyle \frac{2|\stackrel{~}{r}_{11}|}{hv_1}}{\displaystyle _{\mathrm{}}^0}𝑑xCos(2k_1x+\eta _1)`$ (77)
$`+`$ $`{\displaystyle \frac{2|\stackrel{~}{r}_{22}|}{hv_2}}{\displaystyle _{\mathrm{}}^0}𝑑xCos(2k_2x+\eta _2)`$ (78)
$`+`$ $`{\displaystyle \frac{1}{hv_1}}\left(|t_{13}|^2+|t_{14}|^2+\mathrm{}\right)`$ (79)
$`+`$ $`{\displaystyle \frac{1}{hv_2}}(|t_{23}|^2+|t_{24}|^2+\mathrm{})]`$ (80)
Now adding and subtracting the following terms inside $`[]`$,
$`{\displaystyle \frac{|\stackrel{~}{t}_{11}|^2}{hv_1}}{\displaystyle _{\mathrm{}}^0}𝑑x`$, $`{\displaystyle \frac{|\stackrel{~}{t}_{21}|^2}{hv_1}}{\displaystyle _{\mathrm{}}^0}𝑑x`$, $`{\displaystyle \frac{|\stackrel{~}{t}_{22}|^2}{hv_2}}{\displaystyle _{\mathrm{}}^0}𝑑x`$, $`{\displaystyle \frac{|\stackrel{~}{t}_{12}|^2}{hv_2}}{\displaystyle _{\mathrm{}}^0}𝑑x`$,

$`\text{we get},`$ (81)
$`\rho (E)`$ $`=`$ $`2[{\displaystyle \frac{1+|\stackrel{~}{r}_{11}|^2+|\stackrel{~}{r}_{21}|^2+|\stackrel{~}{t}_{11}|^2+|\stackrel{~}{t}_{21}|^2}{hv_1}}{\displaystyle _{\mathrm{}}^0}dx`$ (82)
$`+`$ $`{\displaystyle \frac{1+|\stackrel{~}{r}_{12}|^2+|\stackrel{~}{r}_{22}|^2+|\stackrel{~}{t}_{12}|^2+|\stackrel{~}{t}_{21}|^2}{hv_2}}{\displaystyle _{\mathrm{}}^0}𝑑x`$ (83)
$`+`$ $`{\displaystyle \frac{2|\stackrel{~}{r}_{11}|}{hv_1}}{\displaystyle _{\mathrm{}}^0}𝑑xCos(2k_1x+\eta _1)`$ (84)
$`+`$ $`{\displaystyle \frac{2|\stackrel{~}{r}_{22}|}{hv_2}}{\displaystyle _{\mathrm{}}^0}𝑑xCos(2k_2x+\eta _2)`$ (85)
$`+`$ $`{\displaystyle \frac{1}{hv_1}}\left(|t_{13}|^2+|t_{14}|^2+\mathrm{}\right)`$ (86)
$`+`$ $`{\displaystyle \frac{1}{hv_2}}(|t_{23}|^2+|t_{24}|^2+\mathrm{})]`$ (87)
Now $`|\stackrel{~}{r}_{11}|^2+|\stackrel{~}{r}_{21}|^2+|\stackrel{~}{t}_{11}|^2+|\stackrel{~}{t}_{21}|^2=1`$

and $`|\stackrel{~}{r}_{12}|^2+|\stackrel{~}{r}_{22}|^2+|\stackrel{~}{t}_{12}|^2+|\stackrel{~}{t}_{21}|^2=1`$

$`\text{Thus},`$ (88)
$`\rho (E)`$ $`=`$ $`{\displaystyle \frac{2}{hv_1}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x+{\displaystyle \frac{2}{hv_2}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x`$ (89)
$`+`$ $`{\displaystyle \frac{2|\stackrel{~}{r}_{11}|}{hv_1}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑xCos(2k_1x+\eta _1)`$ (90)
$`+`$ $`{\displaystyle \frac{2|\stackrel{~}{r}_{22}|}{hv_2}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑xCos(2k_2x+\eta _2)`$ (91)
$`+`$ $`{\displaystyle \frac{2}{hv_1}}\left(|t_{13}|^2+|t_{14}|^2+\mathrm{}\right)`$ (92)
$`+`$ $`{\displaystyle \frac{2}{hv_2}}(t_{23}|^2+|t_{24}|^2+\mathrm{})\text{(iii)}`$ (93)
Now $`{\displaystyle \frac{2}{hv_1}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x+{\displaystyle \frac{2}{hv_2}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x=\rho _0(E)`$ i.e. DOS in the absence of scatterer. Also according to some text books , for plane wave states the current conservation condition

$$\frac{}{t}_\mathrm{\Omega }\psi ^{}\psi 𝑑\tau +_S\stackrel{}{J}.\widehat{n}ds=0$$

can be satisfied if and only if the wave function vanishes at $`\pm \mathrm{}`$ for all energy. Thus the 3rd and 4th term of Eq.(iii) of this section is $`0`$. Besides, neglecting these oscillatory terms in the leads is very standard in the derivation of Friedel sum rule , where it is assumed that the carrier concentration is so high in the leads that the leads are non-polarizable. However, assuming non-polarizable leads alone does not eradicate all energy dependence of the self energy as thought in Ref.  (see Eqs. 2 and 6 therein, the parameter $`t_\alpha `$ can be strongly energy dependent, which is dispersive behavior). So we get

$`\rho (E)`$ $`=`$ $`\rho _0(E)+{\displaystyle \frac{2}{hv_1}}\left(|t_{13}|^2+|t_{14}|^2+\mathrm{}\right)`$ (94)
$`+`$ $`{\displaystyle \frac{2}{hv_2}}(|t_{23}|^2+|^2+\mathrm{})`$ (95)
$`\text{Thus}\rho (E)\rho _0(E)`$ $`=`$ $`{\displaystyle \frac{2}{hv_1}}\left(|t_{13}|^2+|t_{14}|^2+\mathrm{}\right)`$ (96)
$`+`$ $`{\displaystyle \frac{2}{hv_2}}(|t_{23}|^2+|^2+\mathrm{})`$ (97)
$`=`$ $`\left(\rho \rho _0\right)_R`$ (98)


format: text

text:
# Submillimeter mapping and analysis of cold dust condensations in the Orion M42 star forming complex

## 1 Introduction

At about 470 pc, the Orion Nebula is one of the nearest massive star formation regions. It is part of the giant molecular complex of Orion, that has been revealed by the 12CO maps of Kutner et al. (1977). The CO map of Maddalena et al. (1986) has shown the global distribution of gas in the Orion-Monoceros complex. This dense region of the interstellar medium, extending over approximatively 30 degrees from southeast to northwest, is composed of several distinct areas. The most northern part is the gaseous region distributed on a circle around $`\lambda `$ Orionis. To the south, the Orion Giant Molecular Clouds A and B extend over 15 degrees from north to south. To the southeast are located Monoceros and the Southern Filament. This set of giant molecular clouds is located 150 pc below the Galactic plane, and might have been formed by the collision of a cloud falling from the southern Galactic hemisphere onto the galactic disk (Franco et al.  (1988)). These molecular clouds are associated with large HII regions, the Orion Ia, Ib, Ic OB star associations, located on a 15 degree long ridge on the side away from the Galactic plane. The Orion Nebula is part of the Orion A giant molecular cloud, and its best studied region. Behind this H II region is the OMC-1 molecular cloud, which corresponds to a CO emission peak. This cloud is heated by the ultraviolet radiation of the OB stars of the Trapezium cluster (see for example Hillenbrand  (1997)). One can distinguish inside of OMC-1, two intense infrared sources: the ponctual object of Becklin-Neugebauer (Becklin & Neugebauer  (1967)) and, 10 arcseconds to the south, the Kleinmann-Low Nebula (Kleinmann & Low  (1967)) that has an angular diameter of about 30 arcseconds. It is thought that BN is a young massive star (25 M) with an important mass loss. The first far infrared observation towards M42 has been performed by Low & Aumann (1970) in the spectral range from 30 $`\mu m`$ to 1000 $`\mu m`$. Then Harper (1974) mapped M42 at 90 $`\mu m`$ and showed that the infrared emission originates from a more extended region. This area has been since intensively studied, from optical radiations to radio wavelengths (see for example the review of Genzel & Stutzki  (1989)). North from OMC-1, one can clearly see an extended ”integral-shaped” filament (ISF hereafter) with two denser regions called OMC-2 and OMC-3 (see Bally et al.  (1987)). These cores have been studied both in molecular lines (recent studies of Dutrey et al.  (1993), Castets & Langer  (1995), Nagahama et al.  (1998)) and continuum (recent works of Chini et al.  (1997), Lis et al.  (1998), Johnstone & Bally  (1999)). The 1300 $`\mu m`$ maps of OMC-2 and OMC-3 of Chini et al. (1997) have shown the precise filamentary structure of these regions, studying several sources inside. On the dense ridge, the three cores OMC-1, OMC-2, OMC-3 have active star formation and seem to have been formed at roughly the same time. However there are important differences between them. OMC-1 has massive star formation and strong interactions between the gas and the infrared sources, increasing the gas temperature to $``$ 70 K, whereas OMC-2 is less luminous, contains low mass young stars and less energetic outflows than OMC-1. OMC-3 is thought to have a mass comparable to OMC-2 ($``$ 100 M), but has less energetic outflows and lower gas temperature, suggesting that OMC-3 is less evolved than OMC-2, itself being less evolved than OMC-1 (see Castets & Langer  (1995)). The recent works in submillimeter continuum (Lis et al.  (1998), Johnstone & Bally  (1999)) have high angular resolution, and thus have clearly shown the complex structure of the Orion Integral-Shaped Filament and discovered small condensations, as also did the $`{}_{}{}^{13}CO`$ work of Nagahama et al. (1998) for the whole Orion A cloud. However the continuum ground-based observations do not have the sensitivity to map the weak emission from the dust away from the dense ridge. A recent far-infrared balloon-borne observation of Orion (Mookerjea et al.  (2000)) has derived the temperature distribution in Orion, the coldest source they found being 15 K.

We have observed this region with the submillimeter balloon-borne telescope ProNaOS, in order to better understand the variations of temperature and spectral index of the dust, and to study the cold phase of the dense interstellar medium. ProNaOS provides simultaneous measurements within four bands covering the wavelength range 200 $`\mu m`$ \- 800 $`\mu m`$, allowing to constrain both the temperature and the spectral index of the dust, in order to derive column densities and masses. The high sensitivity per beam allows us to study the extended emission of the dust, in and a few parsecs away from the dense ridge. A first ProNaOS mapping of OMC-1 and its surroundings has been analyzed by Ristorcelli et al. (1998). They discovered a cold (12.5 K) condensation and showed for four sources the variability of temperature and spectral index in the Orion Nebula complex. We present in this paper a detailed analysis of the submillimeter maps obtained towards Orion with the second flight of ProNaOS in september 1996, at Fort-Sumner, New Mexico. We have observed a much larger region than in Ristorcelli et al. including the integral-shaped filament and other cold clouds, two being seen for the first time in continuum emission. Section 2 deals with the observation performed, section 3 describes the maps obtained, and section 4 presents a detailed analysis of these maps.

## 2 Instruments and observations

### 2.1 Observations with ProNaOS

ProNaOS (PROgramme NAtional d’Observations Submillimétriques) is a French balloon-borne submillimeter experiment, with a 2 m diameter telescope (Buisson & Durand  (1990)). The main characteristics of the instrument are given in Table 1. The focal plane instrument SPM (Système Photométrique Multibande, see Lamarre et al.  (1994)) is composed of a wobbling mirror, providing a beam switching on the sky with an amplitude of about 6’ at 19.5 Hz, and four bolometers cooled at 0.3 K. They measure the submillimeter flux in the spectral ranges 180-240 $`\mu m`$, 240-340 $`\mu m`$, 340-540 $`\mu m`$ and 540-1200 $`\mu m`$, with sensitivity to low brightness gradients of about 1 MJy/sr in band 4. Thanks to an arrangement with dicroic filters, the measurement is performed simultaneously for all four bands on the same sky pixel. An off-axis star tracker and a 3-axis gyroscope allow to control the telescope pointing. Before a flight, the telescope and the whole payload are aligned and controlled. Two internal black bodies are calibrated with an absolute reference. They provide in-flight calibrations with an accuracy of about 5 % (relative), which has been checked against Saturn. The measurements on Saturn allow us to derive the beam shape up to a radius of six arcminutes from the axis. We use the integral of these beams to compare the Saturn measurement to the ground based calibration which is performed on extended blackbodies filling the beam. Regarding the ratio of the main central beam to the total beam, the beam efficiencies are from 0.78 in band 1 to 0.96 in band 4. This can be expressed as a mean wave-front error of about 15 $`\mu m`$. More details about the instrument can be found in Ristorcelli et al. (1998).

The observing procedure is an altazimut scanning of the beam on the sky. The time-ordered data suffer thus from two main artifacts: the measured flux is differential, and the rotations of the telescope around the direction of observation distort the mapping pattern from a rectangular grid.

### 2.2 Data processing: a new method for ProNaOS data

The first processing applied to the data is correction from map distortion, taking into account the pointing of the telescope, including fine pointing errors due to the swinging of the gondola. Then usually the deconvolution method applied to ProNaOS data was EKH-like (Emerson et al.  (1979)) with a scan by scan filtering in Fourier space (see Sales et al.  (1991)). We have developed another method, based on direct linear inversion on the whole map, using a Wiener matrix. The Wiener filter (Wiener  (1949)) is the linear method that minimizes the reconstruction error. The restored sky map is computed as a linear expression of the data:

$`x_0=Wy`$

$`x_0`$ being the reconstructed sky vector (i.e. the map), W the inversion matrix and y the signal vector. Assuming that the instrument is linear, the observation is:

$`y=Ax+n`$

where x is the true sky vector, A the observation matrix and n the noise vector. For W we use a Wiener matrix (see Tegmark  (1997)):

$`W=[S^1+A^tN^1A]^1A^tN^1`$

where S is the sky covariance matrix and N the noise covariance matrix. The convolution matrix A (point-spread function) takes into account the size and profile of the two beams, and the observing mode. In the case of the Orion region presented here, we have gathered two observation sequences. The first one is centered 10’ north from BN/KL, and covers 50’ by 40’. The second is offset to the west to map the edge of the molecular cloud, and covers about 30’ by 20’. In such data, one of the main problems encountered to properly reconstruct the maps, is the great contrast between the intense areas like OMC-1 and the weak areas, which are not less interesting. We have processed separatly the second sequence, in order to decrease the contrast in these data and reconstruct precisely the intensities of the second sequence. The size of A is about 2000\*3000, which allows us to use a classical bi-conjugate gradient algorithm to invert the main term in W. The correlation matrices are estimated from the data, and can be refined with iterations and tests.

We show in Fig. 1 the map of the first sequence obtained with the EKH-like method, to be compared in Fig. 2 to the optimal maps obtained with the Wiener inversion method. Negative artifacts are clearly visible along the scans of the most intense area (along the elevation of OMC-1) with the EKH-like method, whereas they are significantly removed with the Wiener inversion method.

### 2.3 Observations with DiaBolo

Observations of Cloud 2 have been performed with the DiaBolo photometer at the focus of the IRAM 30 meter telescope in december 6/7 1996. DiaBolo is a two channel (1200 $`\mu m`$ and 2100 $`\mu m`$) photometer which uses 100 mK bolometers. A detailed description of the instrumental setup can be found in Benoît et al. (2000). The configuration for the observation of Cloud 2 used one bolometer per wavelength channel with a beam size of about 30 arcseconds. We used the wobbling secondary mirror to subtract the background emission from the sky and the telescope. The wobbling amplitude was set to 120 arcseconds for a frequency of 0.8 Hz. We mapped a 6 by 6 arcminute area centered on the cloud by scanning the telescope beam in azimut and elevation. The scanning speed was 6 arcsec per second, and the step between two lines was 15 arcseconds. The maps were deglitched by filtering the signal of each scan using a temporal filter. After base line subtraction and rebinning of each scan, the dual beam maps were deconvolved from the beam switching using home-made IDL routines which follow the EKH method (Emerson et al.  (1979)).

## 3 Results

We present in Fig. 2 the images obtained in the two extreme photometric bands of ProNaOS-SPM. The noise level is about 4 MJy/sr rms in band 1 and 0.8 MJy/sr in band 4. However, due to the calibration uncertainty, the flux accuracy is not better than 5 % (1 $`\sigma `$) relative between bands (10-20 % absolute).

The brightest area observed is the molecular cloud OMC-1, which appears on our maps as a very intense region of thermal emission of the dust. Particularly in the area of BN/KL, the emission is very strong (it is the maximum of the submillimeter emission in this region): it reaches 49000 MJy/sr in band 1 (200 $`\mu m`$).

North from this intense central core is clearly visible a large molecular cloud extending over about 3 pc, called the integral-shaped filament, which presents an average spectral intensity at 200 $`\mu m`$ of about 2300 MJy/sr.

One can see west from the central core a weaker condensation that we call Cloud 1, which has an intensity of 570 MJy/sr in band 1. A weaker emission can be seen around this cloud, which is linked to the Cloud 2 northwards.

The Clouds 2, 3, 4 can be seen on the west of the maps in Fig. 2, from east to west. The Cloud 2 is a weak intensity condensation discovered during the first flight of ProNaOS (Ristorcelli et al.  (1998)). Its 200 $`\mu m`$ intensity is 47 MJy/sr, but its maximum is closer to band 2 (260 $`\mu m`$), in which the spectral intensity reaches 50 MJy/sr.

The Clouds 3 and 4 are two extremely weak condensations. The 260 $`\mu m`$ intensity of Cloud 3 is 32 MJy/sr, and the 200 $`\mu m`$ intensity of Cloud 4 is 35 MJy/sr. These cold clouds, as well as Cloud 2, cannot be seen with IRAS 100 $`\mu m`$ maps of the same region.

## 4 Analysis

### 4.1 Dust temperatures and spectral indexes

#### 4.1.1 Derivation

To compare the emission obtained in the four submillimeter bands, and particularly to deduce emission spectra of the sources identified, it is necessary to degrade the angular resolution in all channels to the same beam. This is why we smooth the images in bands 1, 2, 3 with an adequate profile in order to obtain for each band the resolution of the fourth band (3.5’). The corresponding beams have been determined from the observation performed on the planet Saturn, considered as a reference point source.

The ProNaOS integrated fluxes in the four bands are given in Table 2. The integrated fluxes from other data (IRAS: http://www.ipac.caltech.edu/ipac/iras/iras.html, Harper  (1974), Chini et al.  (1984), Mezger et al.  (1990), DiaBolo: Benoît et al.  (2000)) are given in Table 3.

We model the emission of the grains with a modified blackbody law:

$`I_{\nu fit}(\lambda ,T,n)=C.B_\nu (\lambda ,T).\lambda ^\beta `$

where $`\lambda `$ is the wavelength, C a constant, T the temperature of the grains, $`\beta `$ the spectral index and $`B_\nu `$ the Planck function.

The three parameters C, T and $`\beta `$, are adjusted with a least square fit. Assuming optically thin emission (as did Ristorcelli et al.  (1998)), we can derive the optical depth $`\tau _\nu =\frac{I_\nu }{B_\nu }=C.\lambda ^\beta `$.

For some regions: Integral-Shaped Filament (south) and ISF (north), we include the IRAS 100 $`\mu m`$ data in the fit. For the other regions, we cannot use these data, because they are either saturated (OMC-1), or drowned in the halo of OMC-1 (Clouds 1,2,3,4). For the fit of OMC-1, we included the 90 $`\mu m`$ data of Harper (1974), the 1000 $`\mu m`$ data of Chini et al. (1984), and the 1300 $`\mu m`$ data of Mezger et al. (1990). These data have been averaged over the same beam diameter as ProNaOS, i.e. 3.5’.

For Cloud 2, we added the data obtained with the millimeter instrument DiaBolo at 1200 $`\mu m`$ and 2100 $`\mu m`$. Three maps were repeated on Cloud 2, for a total observing time of about 1 hour. The values given in Table 3 have been obtained by integrating in a 3.5 arcmin circle centered on the cloud.

We give spectra of some regions in Fig. 3. The best parameters for each region selected are presented in Table 4.

#### 4.1.2 Variations of the temperature and the spectral index

The dust temperature is highly variable all around this molecular complex: from 12 K to 70 K; and the spectral index changes much too: from 1.1 to 2.2, as was also found for the four sources analyzed by Ristorcelli et al. (1998). The clouds separated from the denser regions (Clouds 3,4) are cold and have high indexes, but we can also see that some cold condensations (Cloud 1, 17 K, Cloud 2, 11.8 K) with high indexes can be close to the central region OMC-1 (1.5 pc). For the two sources in common (OMC-1 and Cloud 2), the parameters derived in this paper are well consistent with Ristorcelli et al. The spectral index derived for the OMC-1 cloud is 1.13 with quite small error bars, which is due to the high number of data points we have in the Rayleigh-Jeans tail of the dust emission, the temperature being around 70 K (see Ristorcelli et al.  (1998) for a precise discussion on this region). For such a region, the submillimeter spectrum we observe could come from a mixture of several dust temperatures. It has been shown by Ristorcelli et al. (1998) that the cold dust column density required to explain the fit temperature (83 K) with $`\beta =2`$ was 40 times more massive than the warm component. They concluded that the spectral index of the dominant component had to be around 1. The ISF extended emission is quite cold (22 K - south, 25 K - north) with spectral indexes of about 1.7. In this temperature range, the ProNaOS data points combined with the IRAS 100 $`\mu m`$ point give quite small error bars. We have simulated a temperature distribution with $`\beta =2`$, with a bell function around 17 K having a fwhm of about 9 K. This leads to fit the spectrum with a temperature of 21 K and a spectral index of 1.8. Thus the temperature and spectral index range we observe in the ISF could partly be explained by a distribution of temperatures with a colder central temperature, and a higher index. The clouds 1 and 2 have the highest spectral index we found in this region (2.2), Cloud 2 being the coldest (11.8 K), with error bars considerably reduced, thanks to the DiaBolo millimeter data, with respect to the parameters given in Ristorcelli et al. The Cloud 3 is another cold cloud (13.3 K) with a high index (1.98), however the error bar on the spectral index is large. The Cloud 4 has ill-defined parameters, but seems also to have comparable values of the temperature and the spectral index.

Recently, high angular resolution submillimeter maps of the integral-shaped filament have been obtained from the ground using the CSO at 350 microns (Lis et al.  (1998)) and the JCMT at 450 and 850 microns (Johnstone & Bally  (1999)). However, these maps do not include the fainter regions evidenced here (Clouds 1-4) nor the diffuse emission around the ISF that we show in this article. Comparing their maps to the ProNaOS ones, it is clear that the extended emission of the dust around the brightest sources is missed by these ground-based experiments. Combining their 350 $`\mu m`$ map and the IRAM 1.3 mm map of Chini et al. (1997), Lis et al. derived a map of the dust emissivity index over the brightest region surrounding OMC-1, including Orion-S and the Orion bar. Under the assumption of a single dust temperature (55 K) over this region, they derived emissivity indexes in the range 1.3 to 2.7, the lower values being observed toward the Orion bar (about 1.7) and IRC2 (1.8) while higher values (typically 2.3) appear north of OMC-1. The index values they derived toward the brightest region of the map (IRC2 and Orion-S) is in the range 1.8-2.0 when averaged within the beam size of our observations. Even if we take into account that the temperature they assumed is too low compared to the one we measure (66 K), these values remain clearly higher than the dust emissivity index value we derived for this region ($`\beta `$=1.13). This discrepancy may originate from their use of two distinct datasets with different systematic calibration errors.

### 4.2 Temperature - index inverse correlation

We have computed the 68 % (1 $`\sigma `$) confidence likelihood contours of the (T,$`\beta `$) values for each cloud (see Fig. 4) -N.B.: the likelihood is proportional to $`e^{\frac{\chi ^2}{2}}`$.

The correlation coefficient between T and $`\beta `$ has been computed globally for the studied regions. We have therefore seven couples (T,$`\beta `$) which correlation we analyse by:

$`C=\frac{(T_i\overline{T}).(\beta _i\overline{\beta })}{\sqrt{{\scriptscriptstyle (T_i\overline{T})^2}.{\scriptscriptstyle (\beta _i\overline{\beta })^2}}}`$

A first analysis of the plot in Fig. 4 shows an inverse correlation between the temperature and the spectral index, which means that the regions of high index have the lowest temperature, and the regions of low index have the highest temperature. The correlation coefficient calculated is -0.92. This correlation is, for a part, dominated by the warm area of OMC-1, but the correlation coefficient is still -0.83 when this area is not included.

One can see in Fig. 4 that the likelihood contours are banana-shaped, that is mainly due for OMC-1 to the insensitivity to temperature variations in the Rayleigh-Jeans tail, and for the cold clouds to the relatively bad sensitivity to spectral index variations, as we do not have millimeter data for most of them. However the fit procedure itself could induce an amount of inverse correlation between the temperature and the spectral index. To prove that this correlation is for the most part not due to the fit, we performed simulations. Repeated fits on noisy simulated data with no intrinsic (T,$`\beta `$) correlation show that the fitting procedure itself induces a correlation coefficient lesser in absolute than -0.4. This artifact of the fit is thus insufficient to explain the correlation found on real data.

Our conclusion is that this inverse correlation found between the temperature and the spectral index has to be an intrinsic physical property of the grains.

Laboratory experiments (see Agladze et al.  (1996) and Mennella et al.  (1998)) showed this effect for temperatures down to 10 K, but it had never been systematically shown until now with observations. Agladze et al. (1996) measured absorption spectra of crystalline and amorphous grains between 0.7 and 2.9 mm wavelength. They deduced an inverse correlation between the power-law index $`\beta `$ and the temperature in the temperature range 10-25 K, and attributed it to two level tunnelling processes. The measures of Agladze et al. are for us insufficient to justify our observation in the submillimeter spectral range, because absorption can be very different than in the millimeter range. Mennella et al. (1998) measured the absorption coefficient of cosmic dust analogue grains, crystalline and amorphous, between 20 $`\mu m`$ and 2 mm wavelength, in the temperature range 24-295 K. They deduced an inverse correlation between T and $`\beta `$, and attributed it to two phonon difference processes.

### 4.3 Column densities

#### 4.3.1 Derivation

To calculate the opacity of the grains, the masses and densities of the studied regions, we need to estimate the column density in each pixel of the maps, i.e. the mass of gas and dust along the line of sight. For this we use the dust 100 $`\mu m`$ opacity from Désert et al. (1990). It allows us to estimate the column density from the spectral intensity. We consider only the thermal emission of the big grains, which dominate widely in this spectral range, and we assume that the spectral index does not change in the ProNaOS spectral range. We take into account the variability of the dust spectral index of the different regions. This gives us a simple self-consistent model that allows us to estimate the column density $`N_H`$, as a function of the spectral intensity and the spectral index. We adopt the value of the opacity $`\kappa _{100\mu m}`$ = 0.361 $`cm^2/g`$ (per gram of total medium: gas and dust) at 100 $`\mu m`$, because this value is well constrained from IRAS data. This is calculated from the extinction curves given in Désert et al. The opacity $`\kappa `$ is defined as $`\frac{\tau }{N_H.m_H}`$, thus the relation between the gas column density and the fit parameters is:

$`N_H=\frac{C.\lambda ^\beta }{\kappa .m_H}`$

where $`m_H`$ is the proton mass and $`\kappa `$ the gas and dust opacity.

Then for the Désert et al. value of $`\kappa _{100\mu m}`$ we have:

$`N_H`$ = 1.67 $`10^{24}`$ . C . (100 $`\mu m)^\beta `$

with $`N_H`$ in protons/cm2 and C in $`\mu m^\beta `$.

We make another simple model using the 100 $`\mu m`$ opacity of the Ossenkopf & Henning (1994) model, which is specific to protostellar cores. The opacity introduced then is $`\kappa _{100\mu m}`$ = 1 $`cm^2/g`$. The proportionality relation is:

$`N_H`$ = 6.02 $`10^{23}`$ . C . (100 $`\mu m)^\beta `$

with $`N_H`$ in protons/cm2 and C in $`\mu m^\beta `$.

#### 4.3.2 Analysis

We compare our column density results to those derived from observations of the rotational transition $`J=10`$ of $`{}_{}{}^{13}CO`$, made by Nagahama et al. (1998). We transform the observed $`W_{{}_{}{}^{13}CO}`$ into the column density in cm-2 by using the correspondence established by Bally et al. (1991) in the Orion Nebula. This correspondence has been calculated using comparisons between the IRAS far-infrared dust emission and the CO lines in Orion. Ristorcelli et al. (1998) have taken for the BN/KL region the coefficient: $`N_H/W_{{}_{}{}^{13}CO}`$ = 2 1021 cm-2/(K km/s), assuming an excitation temperature of about 80 K. We have used the same value of this coefficient for estimating the column density of OMC-1 from the $`{}_{}{}^{13}CO`$ data. For the cold and rather cold clouds we study in this article, we choose to adopt $`N_H/W_{{}_{}{}^{13}CO}`$ = 1021 cm-2/(K km/s), as it is mentioned by Bally et al. (1991) for an excitation temperature of about 30 K. Indeed the $`{}_{}{}^{13}CO`$ emission towards the cold clouds we observe ($`<`$ 20 K) is not likely to trace very well the clump itself, but rather the surrounding gas, which has probably a quite higher temperature. Recent CO observations of Ristorcelli et al. (2001, in preparation) of Cloud 2 seem to support this, as the $`{}_{}{}^{13}CO`$ emission of this cloud is more diffuse than the $`C^{18}O`$ emission. This method gives a simple way to estimate the gas column density from these $`{}_{}{}^{13}CO`$ data. As we can see in Table 5, for OMC-1, ISF-N and ISF-S, there is a rather good agreement between NH estimated from ProNaOS data by the model from Désert et al. and NH from 13CO data. Indeed, the $`\frac{N_Hdust}{N_H{}_{}{}^{13}CO}`$ ratio is 1.3 for OMC-1, 0.9 for ISF south and 1.5 for ISF north. Of course the uncertainties on both values are large and difficult to estimate, but these results prove that these values found are quite robust. The cold clouds (1,2,3,4) have 13CO column densities in rather good agreement with the ProNaOS Ossenkopf & Henning estimations: $`\frac{N_Hdust}{N_H{}_{}{}^{13}CO}`$ is 0.7 for Cloud 1, 1.2 for Cloud 2, 1.7 for Cloud 3 and 2.4 for Cloud 4. This rather good agreement between the $`{}_{}{}^{13}CO`$ estimated column densities and the O&H estimated ones may be explained by the fact that some cold condensations are the site of formation of molecular ice mantles on the grains, and of coagulation of grains. These processes are taken into account by the model of Ossenkopf & Henning. However, the 13CO estimations for NH seem under-estimated for Clouds 3 and 4, compared to the ProNaOS O&H values ($`\frac{N_Hdust}{N_H{}_{}{}^{13}CO}`$ ratios of 1.7 and 2.4). If we suppose a higher excitation temperature than 30 K for these two clouds, it is possible to match the ProNaOS O&H values, for instance assuming T$``$ 80 K leads to $`10^{21}`$ $`cm^2`$ for the 13CO NH value of Cloud 4, but this warm temperature is not likely for this region. Another possibility may be that the O$`\&`$H opacity value still over-estimates the column density for these clouds, that may mean that they are the place of even more grain special effects than those taken into account by Ossenkopf & Henning.

We also show in Table 5 the mass of each region. This is obtained simply by multiplying the estimated gas column density by the surface of the region, assuming a distance of 470 pc.

For the clouds separated from the molecular complex, we may reasonably ask if they can collapse due to their mass and therefore become protostellar objects. As an indicator of the stability of the clouds, we give in Table 5 the Jeans mass. We use for this the expression derived from the equality of gravitational and thermal energy (Larson  (1969)):

$`M_J`$ = 10-18 . $`\frac{D_{cm}}{2}`$ . T

with MJ in solar masses and $`D_{cm}`$ the diameter in cm of the cloud.

T refers in this equation to the gas temperature, whereas we used the dust temperature found. In the dense interstellar medium, the two temperatures are weakly different (Hollenbach  (1988), Tielens & Hollenbach  (1985)). We assume by this way that the mass is essentially the clump’s one, and that the gas temperature of the clump is about the same as the dust one, but it may be different from the temperature of the supposed surrounding gas (traced by the 13CO).

The comparison of the cloud mass to the Jeans mass is here to be taken as an indicator of the cloud’s instability. There are of course large uncertainties in both the clump’s measured mass and the non thermal sources of internal energy, such as magnetic field and turbulence, which can help to balance the gravitational energy. On the other hand, external pressure may lead to stabilization of otherwise evaporating clumps. In our sample of four cold clouds, there is a trend that the clumps closer to the active region (Clouds 1 and 2) are potentially the more unstable.

## 5 Conclusion and summary

Our study shows a large distribution of temperatures and spectral indexes in and around a dense and active molecular complex, the M42 Orion Nebula. The temperature varies from 12 K to 70 K, and the spectral index from 1.1 to 2.2. The finding of two new cold clouds (Clouds 3 and 4) confirms that the existence of cold condensations in such regions is not unusual. However the extended cold clumps are located in the outskirts of the active star forming area. They may be the sites for future star formation.

The statistical analysis of the temperature and spectral index spatial distribution shows an evidence for an inverse correlation between these two parameters. This effect is not well explained yet, especially in the submillimeter spectral range for cold grains ($`<`$ 20 K). It has been shown to occur in the laboratory for warm grains by Mennella et al. (1998), and for cold grains in the millimeter by Agladze et al. (1996).

We estimated the column densities and masses of the observed regions, by simply modelling the thermal emission of the grains from Désert et al. (1990) or Ossenkopf & Henning (1994). There is a good agreement between the $`{}_{}{}^{13}CO`$ column densities and those derived from our submillimeter measurement. This demonstrates the robustness of dust opacity values in the grain models. The submillimeter wide band spectro-imaging is thus a natural way to derive masses in the interstellar medium.

Finally we see a trend that the closer to the complex the cold clouds are, the more unstable they are. The history of star formation around OMC-1 shows that there have already been 3 to 4 successive bursts of star formation in this region with the embedded cluster responsible for the BN/KL object being the latest. The clouds that we observe close to the active region may thus be the seeds of the next generation of stars. This is of course quite speculative and should be sustained by more observations, particularly of the possible embedded protostars.

## 6 Acknowledgements

We wish to thank Drs Nagahama, Mizuno, Ogawa and Fukui, for kindly having provided us their CO maps. We are indebted to the French space agency Centre National d’Études Spatiales (CNES), which supported the ProNaOS project. We are very grateful to the ProNaOS technical teams at the CNRS and CNES, and to the NASA-NSBF balloon-launching facilities group of Fort Sumner (New Mexico).