text:
# 1 Introduction

## 1 Introduction

Turbulence is the state of vortical fluid flow where velocity, pressure, and other flow field properties vary in time and space sharply and irregularly, and can be assumed to be random. The experimental investigation of individual realizations of such flows is impossible because the results are irreproducible: Experiments repeated under identical external conditions produce a different outcome. Therefore experimental investigations of turbulent flows can only provide their average properties.

Less trivial and not always recognized is the following: What is of greatest interest in these experiments are intermediate asymptotic states of wider classes of flows, i.e., coherent, self-consistent fragments common to many different flows. Other measurements reflect special properties of a set-up, which cannot be reproduced in other experiments. It may sometimes be useful to investigate a particular device (e.g. an atomizer) for an immediate practical purpose, but one should be cautious in transferring the results to different flows.

Typical examples of intermediate-asymptotic flows are shear flows, where the flow is homogeneous in the direction of mean velocity which depends only on the coordinate perpendicular to the mean flow. A well known example of shear flows occurs in smooth cylindrical circular pipes far from the entrance and outlet. Another example is the zero-pressure gradient boundary layer above a smooth flat plate far from its tip. In spite of their apparent simplicity experiments with such flows require high experimental culture and are expensive, and therefore relatively rare. When they are successful, like, for instance, the experiments of Nikuradze \[1\] with flows in smooth pipes in the range of Reynolds numbers up to $`3.2410^6`$, performed 70 years ago under the guidance of L.Prandtl, they become milestones in turbulence studies. They are used to check theories based on special hypotheses valid for special classes of flows. Not enough is known now about the solutions of Navier-Stokes equations to avoid such hypotheses.

The responsibility of experimentalists who perform such experiments and process the data is therefore very high. They are to be very careful in their conclusions. We showed, for instance \[2,3\] that the experiments of Princeton group (A.Smits, M.V.Zagarola) presented in the thesis of M.V.Zagarola \[4\], which attempted to increase the range of Reynolds numbers achieved by Nikuradze by an order of magnitude, had a flaw. Starting from $`Re=10^6`$, well below the upper boundary of Nikuradze experiments, their data were influenced by roughness — insufficient polishing of the pipe walls.

In the present work we analyze the experiments with zero-pressure-gradient boundary layers presented in the thesis of J.M.Österlund \[5\]. Like the earlier theses of M.V.Zagarola \[4\] and M.H.Hites \[6\], the thesis of J.M.Österlund presents the results of long-time work of a group headed by a senior scientist (in this case A.V.Johansson), using a complicated, expensive and unique facility. The experimental data and their interpretation presented in this thesis might be accepted by some readers, as in the case of the thesis by Zagarola, without precautions. This was the motivation of our analysis of this thesis.

The number of runs (70 measurements of mean velocity profiles) reported in the thesis \[5\] is larger than in previously reported series, although the range of Reynolds numbers covered is not as wide; less than in the previous thesis of M.H.Hites \[6\]: $`2500<Re_\theta <27,500`$. In the thesis \[5\] the authors make very definite statements: they claim that their experiments confirm the classical two-layer theory, in particular the Reynolds number-independent universal logarithmic law, and, exhibiting no Reynolds number-dependence, disagree with the alternative theory based on the Reynolds number-dependent scaling law.

J.M.Österlund presented the data of 70 mean velocity measurements on the Internet (www.mesh.kth.se/$``$jens/zpg/ ). We present here the results of the processing of all these data. We demonstrate that, properly processed, these data lead to the opposite conclusion: they confirm the Reynolds number-dependent scaling law and disagree with the conclusion of Reynolds number-independence.

## 2 Background

According to the classical two-layer theory of wall-bounded turbulent shear flows at large Reynolds numbers, the distribution of average velocity $`u`$ across the flow in the basic intermediate region adjacent to the viscous sublayer is represented in the form of universal (Reynolds number-independent) von Kármán-Prandtl logarithmic law

$$\varphi =\frac{1}{\kappa }\mathrm{ln}\eta +C.$$
(2.1)
Here we use classical Nikuradze-Schlichting et al. notations:

$$\varphi =\frac{u}{u_{}},u_{}=\left(\frac{\tau }{\rho }\right)^{\frac{1}{2}},\eta =\frac{u_{}y}{\nu },$$
(2.2)
where $`\tau `$ is the shear stress at the wall, $`y`$ the distance from the wall; $`\nu `$ and $`\rho `$ are the fluid’s kinematic viscosity and density. In the thesis more modern notations are used: $`\overline{U}^+`$ instead of $`\varphi `$, $`y^+`$ instead of $`\eta `$. So, the von Kármán-Prandtl universal law (2.1) is represented in the thesis in the form

$$\overline{U}^+=\frac{1}{\kappa }\mathrm{ln}(y^+)+B.$$
(2.3)
Here $`\kappa `$ (von Kármán constant) and $`B`$, according to the logic of the derivation, should be universal constants identical in all high quality experiments. It is known from the literature however that various experiments give substantially different values of these constants. Nikuradze \[1\] determined $`\kappa =0.417`$, $`B=5.84`$, Monin and Yaglom \[7\] give the values $`\kappa =0.40`$, $`B=5.1`$; Schlichting \[8\] gives the values $`\kappa =0.40`$, $`B=5.5`$. The difference is substantial, and for many years doubts have accumulated on the validity of the universal logarithmic law.

In our papers (see e.g. \[9\],\[10\]) the derivation of the universal logarithmic law was reconsidered. It was shown that one of the basic assumptions is not quite correct, and on the basis of an alternative assumption, a different “scaling” (power) law was proposed:

$$\overline{U}^+=(C_1\mathrm{ln}Re+C_2)(y^+)^{c/\mathrm{ln}Re}$$
(2.4)
where the constants $`C_1,C_2,c`$ should be universal and Reynolds number-independent. Comparison with the experimental data of Nikuradze has given the following values of the constants:

$$C_1=\frac{1}{\sqrt{3}},C_2=\frac{5}{2},c=\frac{3}{2},$$
(2.5)
so that the law (2.4) is presented in the form

$$\varphi =\left(\frac{\sqrt{3}+5\alpha }{2\alpha }\right)\eta ^\alpha ,\alpha =\frac{3}{2\mathrm{ln}Re}$$
(2.6)
or, using the notation of Österlund’s thesis \[5\]

$$\overline{U}^+=\left(\frac{1}{\sqrt{3}}\mathrm{ln}Re+\frac{5}{2}\right)(y^+)^{\frac{3}{2\mathrm{ln}Re}}.$$
(2.7)
Asymptotically, at $`Re\mathrm{}`$, the specific choice of $`Re`$ is of no importance: $`Re`$ can be replaced by $`Re^{}=\text{Const}Re`$, and the asymptotic form of (2.7) will remain the same. However, for practical large, but not too large, values of $`Re`$, its actual expression is significant. It should be remembered that in our comparison with Nikuradze’s data \[1\] we used his definition of the Reynolds number: $`Re=\overline{u}d/\nu `$, where $`\overline{u}`$ is the mean flow velocity (bulk discharge rate divided by cross-section area), and $`d`$ is the diameter of the pipe. Furthermore, it was recognized from the beginning \[11\] that the law (2.7) is asymptotic (in the parameter $`1/\mathrm{ln}Re`$). It should be valid at large $`Re`$, but at lesser values higher terms in the expansion of the coefficients could be significant.

The experimental data of the Stockholm group (J.M.Österlund, A.V.Johansson) are presented in the form of graphs in the $`\mathrm{ln}y^+`$, $`\overline{U}^+`$ plane suggested by the universal logarithmic law (2.3) (see Figure 5 on page 43 and Figures 13, 14 on pages 152–153 of the thesis \[5\]). They became available on the Internet (www.mesh.kth.se/$``$jens/zpg/ ) in complete digital form. The data are presented in our table, parameterized by the authors by the parameter $`Re_\theta =U\theta /\nu `$. Here $`U`$ is the free stream velocity, $`\theta `$ is the momentum thickness, a quantity measurable a posteriori, after the run is performed.

## 3 Processing of the mean velocity data

The first question to be answered was as follows: are the mean velocity data presented on the Internet compatible with some scaling law, not necessarily the law (2.7). Therefore the data were plotted in the double-logarithmic coordinates $`(\mathrm{lg}y^+,\mathrm{lg}\overline{U}^+)`$. The result was instructive: the data outside the viscous sublayer form a characteristic shape of a broken line (see Figures 1–70). This shape is similar to the shape obtained for the experiments of the first group according to our classification \[12\] where there was no influence neither of the external turbulence of the free stream nor of roughness. The Stockholm authors place the lower boundary of the intermediate region at $`y^+=200`$. We found this value generally too high, and the standard value $`y^+=70÷100`$ seems to be more appropriate, however we marked the line $`y^+=200`$ on all Figures 1–70.

Thus, all experiments revealed two straight lines forming a broken line in the $`\mathrm{lg}y^+,\mathrm{lg}U^+`$ plane. These straight lines correspond to the scaling laws

$$\text{(I)}\overline{U^+}=A(y^+)^\alpha $$
(3.1)
(in the intermediate region adjacent to the viscous sublayer), and

$$\text{(II)}\overline{U^+}=B(y^+)^\beta $$
(3.2)
(in the intermediate region adjacent to the free stream). The coefficients $`A,\alpha ,B,\beta `$ were obtained by standard statistical processing of data (see the table). The coefficients $`A,\alpha `$ of the straight line (3.1) representing the scaling law for the mean velocity distribution in the basic intermediate region adjacent to the viscous sublayer are obviously Reynolds-number-dependent: For us this was not unexpected, because previous processing of all available experimental data for a much wider range of Reynolds number led us to the same conclusion (see paper \[12\] which was known to the Stockholm group and referenced by them). Therefore we conclude that the validity of some Reynolds-number-dependent scaling law for mean velocity distribution is unquestionably confirmed by the experiments of the Stockholm group as well.

Note that because the Reynolds number range covered by the Stockholm group was not large, substantially less than the range covered by the other groups, in particular the Chicago group \[6\], there would be a danger that they would not notice the Reynolds-number-dependence because the governing parameter is $`\mathrm{ln}Re`$, not $`Re`$ itself. This is not the case: The $`Re`$-dependence of the Stockholm experimental data is sufficiently strong to be revealed by proper processing.

By the way, the authors could notice that their values $`\kappa =0.38`$ and $`\beta =4.1`$ are substantially less than those presented in the literature as standard. However, if the law is universal Reynolds-number-independent, these parameters should be identical for all experiments of sufficient quality!

The argument against the power law used by the authors (see paper \[13\] reproduced in the thesis) is the following. They introduce the “diagnostic function”

$$\mathrm{\Gamma }=\frac{y^+}{\overline{U}^+}\frac{d\overline{U}^+}{dy^+}.$$

Their statement, “The function $`\mathrm{\Gamma }`$ should be a constant in a region governed by a power law” is correct for a fixed Reynolds number. However, this is not true for the “diagnostic function averaged for KTH data”, which is shown in their Figure 6.

We invite the reader to look at any of the Figures 1–70. It is clear that for each run $`\mathrm{\Gamma }`$ is a constant — look at the straight lines in the first intermediate region! However, this constant is different for different runs because the slopes of straight lines is $`Re`$-dependent! Indeed, the slope in the first region decays with growing Reynolds number. It is clear why $`\mathrm{\Gamma }`$ obtained by the authors is decreasing: the runs with larger Reynolds number and smaller slopes contribute more to larger $`y^+`$.

Now, when the validity of some Reynolds-number-dependent scaling law for the experiments of the Stockholm group is unquestionably established, we have to investigate whether this scaling law can be represented in the same form (2.7) obtained for flows in pipes. But what is Re? We cannot assume it arbitrarily to be equal to $`Re_\theta `$.

This effective Reynolds number $`Re`$ should have the form $`Re=U\mathrm{\Lambda }/\nu `$, where $`U`$ is the free stream velocity, $`\nu `$ the kinematic viscosity, and $`\mathrm{\Lambda }`$ a length scale which we cannot à priori identify with the momentum thickness $`\theta `$, as there is no rationale for such identification. So, the basic question is whether one can find for each run such length scale $`\mathrm{\Lambda }`$ so that the scaling law (2.7) will be valid for the mean velocity distribution in the first intermediate region. A priori the very existence of such a length scale is under question, but if it does exist, this means that the law (2.7) is not specific to flows in pipes but can be a general law for wall-bounded shear flows at large Reynolds numbers.

To answer this question one has to take the values $`A`$ and $`\alpha `$ for each run, obtained by standard statistical processing of the experimental data in the first intermediate scaling region, and then calculate two values $`\mathrm{ln}Re_1`$ and $`\mathrm{ln}Re_2`$ by solving two equations suggested by the scaling law (2.7),

$$\frac{1}{\sqrt{3}}\mathrm{ln}Re_1+\frac{5}{2}=A,\frac{3}{2\mathrm{ln}Re_2}=\alpha .$$
(3.3)
If the values $`\mathrm{ln}Re_1`$ and $`\mathrm{ln}Re_2`$ obtained by solving two different equations (3.3) are close, i.e. if they coincide within experimental accuracy, it will mean that the unique length scale $`\mathrm{\Lambda }`$ can be determined so that the experimental scaling law in the region (3.1), whose existence was proved before, coincides with the basic scaling law (2.7). The table shows that indeed these values are close — for all $`Re_\theta >10,000`$, the difference $`\mathrm{ln}Re_1\mathrm{ln}Re_2`$ does not exceed 3%. This allows one to introduce for large Reynolds numbers the effective Reynolds number $`Re`$ according to the relation

$$\mathrm{ln}Re_1=\frac{1}{2}(\mathrm{ln}Re_1+\mathrm{ln}Re_2),\text{or}Re=\sqrt{Re_1Re_2},$$
(3.4)
i.e., the geometric mean of $`Re_1`$ and $`Re_2`$. This Reynolds number allows the definition of the effective length scale $`\mathrm{\Lambda }`$, which plays a similar role in the scaling law for the boundary layer flow as does the pipe diameter for flow in pipes. We remind the reader that the momentum thickness is calculated by integration of the velocity profile obtained experimentally: the calculation of the length scale on the basis of the measured velocity profile is not more complicated. Naturally the ratio of two length scales $`\theta /\mathrm{\Lambda }`$ is different for different runs: both these quantities depend upon the details of flows, in particular they can depend in principle upon the distance between the tip of the plate and the point of observation.

## 4 Checking universality

The scaling law (2.7) can be reduced to a universal form

$$\psi =\frac{1}{\alpha }\mathrm{ln}\left(\frac{2\alpha \overline{U}^+}{\sqrt{3}+5\alpha }\right)=\mathrm{ln}y^+$$
(4.1)
where $`\alpha =\frac{3}{2\mathrm{ln}Re}`$. This formula gives another way to check the applicability of the Reynolds-number-dependent scaling law (2.7) in the intermediate region (3.1). Indeed, according to (4.1), in the coordinates $`\mathrm{ln}y^+,\psi `$, all experimental points should collapse onto the bisectrix of the first quadrant. Figure 71 shows that all data for large Reynolds numbers ($`Re_\theta >15,000`$, 24 runs) presented on the Internet collapse onto the bisectrix with accuracy sufficient to give an additional confirmation to the Reynolds-number-dependent scaling law (2.7). For lesser values of $`Re_\theta `$ a systematic parallel shift is observed (Figures 72–74). Apparently in this case the choice of $`Re`$ according to (3.4) is insufficient because at small Reynolds numbers the higher terms of the expansion could have some influence (see the paper by Radhakrishnan Srinivasan \[13\]).

## 5 Conclusion

The thesis of J.Österlund contains the following statements:

1. (p.22 of thesis), “The classical two layer theory was confirmed and constant values of the slope of the logarithmic overlap region (i.e. von Kármán constant) and additive constants were found and estimated to $`\kappa =0.38`$, $`B=4.1`$, and $`B_1=3.6`$.”
2. (p.29 of thesis), in fact the Introduction to the paper : “Contrary to the conclusions of some earlier publications, careful analysis of the data reveals no significant Reynolds number dependence for the parameters describing the overlap region using the classical logarithmic relation.”

These statements are not correct. Our careful analysis of experimental data presented in the thesis performed in the present work leads to the opposite conclusions.

> $`1^{}.`$ The results contradict the classical two-layer theory. The estimates of the constants obtained by the authors are substantially different from the standard values and this reason alone is enough to reject the assumption of universality, the cornerstone of the classical theory.
>
> $`2^{}.`$ In full agreement with our earlier publications, careful analysis of the data reveals significant Reynolds number dependence for the parameters describing the “overlap” region and confirms the Reynolds-number-dependent scaling law.

The thesis of J.Österlund was not the first investigation of this kind, many similar experimental investigations were performed earlier covering a much larger range of Reynolds number. (In the thesis only one decade of $`Re_\theta `$ was covered: $`2,500<Re_\theta <27,000`$; in previous investigations, in particular in those reflected in the instructive review \[15\], the range $`1,000<Re_\theta <200,000`$ was covered). However, as we showed above, the accuracy of experimental data is sufficient to reveal Reynolds number dependence and correspondence of their data to the Reynolds-number-dependent scaling law (2.7) proposed by us.

The experiments of Österlund et al. allow an additional confirmation of the separation of the basic part of the flow into two self-similar regions (I) and (II) governed by the laws (3.1) and (3.2). It is important that these two self-similar regions cover the whole boundary layer and not a small (1/6!) part of it where the universal law is expected by the authors to be valid. These experiments reveal a weak $`Re`$-number dependence of the parameter $`\beta `$ (Figure 75): it decreased with growing $`Re`$. The data are not sufficient to come to a final decision, but they are in approximate agreement with the correlation

$$\beta =\frac{2}{\mathrm{ln}Re}+0.01.$$
(5.1)
Acknowledgments. This work was supported in part by the Applied Mathematics subprogram of the U.S. Department of Energy under contract DE–AC03–76–SF00098, and in part by the National Science Foudnation under grants DMS 94–16431 and DMS 97–32710.

References

1. Nikuradze, J. 1932. Zur turbulenten Strömung in glatten Rohren. VDI Forschungsheft, no. 356.
2. Barenblatt, G.I., Chorin, A.J., and Prostokishin, V.M. 1997. Scaling laws in fully developed turbulent pipe flow: discussion of experimental data. Proceedings National Academy of Sciences USA 94a, pp.773–776.
3. Barenblatt, G.I., Chorin, A.J. 1998. Scaling of the intermediate region of wall-bounded turbulence: The power law. Physics of Fluids vol. 10, no. 4, pp.1043–1044.
4. Zagarola, M.V. 1996. Mean-flow scaling of turbulent pipe flow. Doctoral Thesis, Princeton University, Princeton, New Jersey.
5. Österlund, J.M. 1999. Experimental studies of zero pressure-gradient turbulent boundary layer flow. Doctoral Thesis, Royal Institute of Technology, Stockholm.
6. Nagib, H. 1997. Scaling of high Reynolds number turbulent boundary layers in the National Diagnostic Facility. Doctoral Thesis, Illinois Institute of Technology, Chicago.
7. Monin, A.S. and Yaglom, A.M. 1971. Statistical Fluid Mechanics, vol. 1, MIT Press, Boston.
8. Schlichting, H. 1968. Boundary Layer Theory, McGraw-Hill, New York.
9. Barenblatt, G.I., Chorin, A.J., and Prostokishin, V.M. 1997. Scaling laws in fully developed turbulent pipe flow. Applied Mechanics Reviews vol. 50, no. 7, pp.413–429.
10. Chorin. A.J. 1998. New perspectives in turbulence. Quarterly of Applied Mathematics, vol. LVI, no. 4, pp.767–785.
11. Barenblatt, G.I., 1993. Scaling laws for fully developed turbulent shear flows. Part 1: Basic hypotheses and analysis. Journal of Fluid Mechanics, vol. 248, pp.513–520.
12. Barenblatt, G.I., Chorin, A.J., and Prostokishin, V.M. 1999. Self-similar intermediate structures in turbulent boundary layers at large Reynolds numbers. UC Berkeley CPAM preprint 775. Journal of Fluid Mechanics, in press.
13. Srinivasan, Radhakrishnan, 1998. The importance of higher-order effects in the Barenblatt-Chorin theory of wall-bounded fully developed turbulent shear flows. Physics of Fluids vol. 10, no. 4, pp.1037–1039.
14. Österlund, J.M., Johansson, A.V., Nagib, H.M. and Hites, M.H. 2000. A note on the overlap region in turbulent boundary layers. Physics of Fluids, vol. 12, no. 1, pp.1–4.
15. Fernholz, H. H. and Finley, P. J., 1996. The incompressible zero-pressure-gradient turbulent boundary layer: an assessment of the data. Progr. Aerospace Sci. 32, pp.245–311.


format: text

text:
# Linear and Nonlinear Experimental Regimes of Stochastic Resonance

## I Introduction.

The renewed interest of the last two decades on stochastic processes modeling different phenomena of physics, chemistry and engineering sciences has led to the discovery of noise-induced phenomena in nonlinear systems away from equilibrium. In these systems a variation of the level of external noise can qualitatively change the response of the system. The paradigmatic example of these noise-induced phenomena is stochastic resonance  (for a recent review see ). Other examples of noise-induced phenomena comprise resonant activation , noise-induced transitions  and noise enhanced stability .

Stochastic resonance (SR) manifests itself as an enhancement of the system response for certain finite values of the noise strength. In particular the signal-to-noise ratio (SNR) shows a maximum as a function of the noise intensity. In other words, a statistical synchronization of the random transitions between the two metastable states of the nonlinear system takes place in the presence of an external weak periodic force and noise. Such a physical system presents a time-scale matching condition, which can be observed by tuning the noise level to such a value that the period of the driving force approximately equals twice the noise-induced escape time. The SR phenomenon appears in a large variety of physical systems and has been observed in different systems, ranging from sets of neurons , to lasers  and to solid-state devices, like SQUIDs and tunnel diodes .

The SR phenomenon is a well investigated phenomenon both from a theoretical and an experimental point of view . However, few studies systematically analyze the SR phenomenon for different values of the frequency and amplitude of the modulating signal . In this article we systematically study the SR phenomenon as a function of the frequency and amplitude of the modulating signal in a physical bistable system based on a tunnel diode. Our experimental set-up allows us to investigate this phenomenon in a range of amplitude and frequency spanning several orders of magnitude. By varying the amplitude and the frequency of the modulating signal we detect both the regime of the SR phenomenon described by the linear response theory and the nonlinear deviation from it. In the linear regime we observe the customary behavior of stochastic resonance whereas in the nonlinear regime we detect a saturation of the power spectral density measured at the frequency of the modulating signal and a depletion of the power spectral density of the noise at the same frequency. When the noise depletion takes place we observe phase and frequency synchronization between the stochastic output and the deterministic input.

The paper is organized as follows. In section II we describe the experimental apparatus and we discuss the stochastic differential equation associated to the electronic circuit based on a tunnel diode. In section III we present our experimental results of the power amplification and SNR as a function of the noise intensity for different values of the amplitude and frequency of the modulating signal. In this section we discuss the detected deviations form the predictions of the linear response theory and we present an evidence of phase and frequency synchronization detected in the nonlinear regime of high values of the amplitude of the modulating signal. In section IV we briefly draw our conclusions.

## II Experimental apparatus and the tunnel diode

The experimental setup used for investigating the SR phenomenon is a bistable electronic system based on a tunnel diode. The physical system is a series of a resistor(tunable to a desired value) and a tunnel diode in parallel to a capacitor. The tunnel diode is a highly doped semiconductor device with a typical current-voltage characteristic showing a region of negative differential resistance, which is due to a tunneling current from the valence band of the n-doped region to the conducting band of the p-doped region. There are two stable states and one unstable state . For details about this experimental set-up see ref. .

A network of general purpose very low-noise wide band operational amplifier is used to sum the driving periodic signal and the noise signal. The noise signal is the output of a commercial digital noise source, whose spectral density is approximately flat up to 20 MHz and whose root mean square voltage $`V_{rms}`$ may be selected within the range from $`0.133`$ to $`5.5`$ V with a 27 mV resolution. At the output of the operational amplifier the statistical properties of the noise are altered by the filtering of the operational amplifier. We measure the noise $`v_n(t)`$ at the output of the operational amplifier and we observe that it is a Gaussian noise characterized by a spectral density which is flat at low frequency ($`f<1.25`$ MHz). By defining the correlation time of the Gaussian noise as the time at which the normalized autocorrelation function assumes the value $`1/e`$, we measure $`\tau _n=120.`$ ns. In our measurements we vary the amplitude and the frequency of the driving periodic signal $`v_s(t)=V_scos(2\pi f_st)`$. Specifically, we vary the amplitude $`V_s`$ from $`0.0067`$ to $`1.00`$ V, and the frequency $`f_s`$ from $`1`$ Hz to $`1`$ MHz. The output voltage across the diode $`v_d`$ is detected by a digital oscilloscope and transferred to a PC. A typical time series has 4096 records. The digitized time series are analyzed on-line by using a fast Fourier transform (FFT) routine. The result of the FFT routine provides the power spectral density of the $`v_d`$ time signal.

We model our electronic circuit by writing down its differential equation. This equation is

$$\frac{dv_d}{dt}=\frac{dU(v_d,t)}{dv_d}+\frac{1}{RC}v_n(t),$$
(1)
which is formally equivalent to a stochastic differential equation describing the position of an overdamped random particle moving in a generalized potential. In this equation $`R`$ is the biasing resistor, $`C`$ is the parallel capacitor (in our case 45 pF) of the circuit and $`v_n(t)`$ is a noise voltage mimicking the presence of a finite temperature in the corresponding overdamped system with a physical particle. The generalized potential $`U(v_d,t)`$ associated to our physical system is

$`U(v_d,t)={\displaystyle \frac{V_bv_d}{RC}}+{\displaystyle \frac{v_d^2}{2RC}}+{\displaystyle \frac{1}{C}}{\displaystyle _0^{v_d}}I(v)𝑑v`$ (2)
$`{\displaystyle \frac{V_sv_d}{RC}}\mathrm{sin}(\omega _st),`$ (3)
here $`V_B`$ is the biasing voltage of the electronic bistable network.

Our circuit presents two control parameters affecting the shape of the associated generalized potential. They are the biasing resistance $`R`$ and the biasing voltage $`V_b`$. We control both of them independently. We perform our experiments by ensuring a symmetric escape from one potential well to the other. This is done by selecting the values of the two control parameters (R=770 $`\mathrm{\Omega }`$ and $`V_b=6.76`$ V) in such a way that we do not detect power spectral density of the output voltage $`v_d`$ above the noise level at even harmonics of the frequency of modulating signal. We also verify that for this choice of control parameters the experimentally measured residence times have approximately the same value in the two potential wells.

The aim of this study is to investigate the SR phenomenon over a wide interval of the frequency $`f_s`$ of the modulating signal. To perform such a task, we have to overcome two experimental conflicting constraints: (i) we are forced to set the time constant of our system $`\tau RC`$ to a low value satisfying the inequality $`f_s1/2\pi \tau `$ and (ii) we need to use a high value of $`\tau `$ to maintain the ratio $`\tau _n/\tau `$ as low as possible to conduct our experiments in the“white noise” limit of $`v_n(t)`$. The best compromise we find is to set $`\tau RC=34.6`$ ns. With this choice $`1/2\pi \tau 4f_s^{max}`$ and $`\tau _n/\tau 3.47`$. In other words we guarantee the investigation of the SR phenomenon over a rather wide range of $`f_s`$ by performing our experiments in a regime of moderately colored noise.

## III Stochastic Resonance for different values of frequency and amplitude of the modulating signal

A bistable system based on a tunnel diode provides a versatile physical system in which stochastic resonance can be investigated . In this study we perform an investigation of the SR phenomenon as a function of a wide range of amplitude and frequency of the modulating signal. The first investigation concerns the power $`P_1`$ of the output signal $`v_d(t)`$ localized around the frequency of the modulating signal. This quantity is obtained by integrating the spectral density $`S(\omega )`$ over the delta-like peak observed at angular frequency $`2\pi f_s`$. The signal “power” $`P_1`$ used in the theory of linear signal processing obtained by integrating the spectral density over the $`f`$ peaks at $`f_s`$ and -$`f_s`$ is

$$P_1=4\pi |M_1|^2,$$
(4)
here $`|M_1|`$ is the magnitude of coefficient of the Fourier series $`<v_d(t)>=_{\mathrm{}}^{\mathrm{}}M_n\mathrm{exp}(in2\pi f_st)`$ taken at the frequency of the modulating signal. The linear response theory for stochastic resonance predicts that $`|M_1|`$ is proportional to $`V_s`$ at fixed value of $`V_n`$ . We test this dependence on a large interval of values of the amplitude signal $`V_s`$. In Fig. 1 we show the measured values of $`|M_1|`$ as a function of $`V_s`$ varying in the interval from 0.0067 to 1.00 V. The measure are done by setting $`f_s=10`$ Hz and $`V_{rms}=1.89`$ V. From the figure is evident that the prediction $`|M_1|V_s`$ obtained by using the linear response theory is valid only within the amplitude interval $`0.017<V_s<0.067`$ V. The deviation observed for the lowest investigated value of the modulating amplitude signal ($`V_s=0.0067`$ V) is probably due to experimental detection problems related with the low value of the signal whereas the deviation observed when $`V_s>0.067`$ is entirely ascribed to a deviation of the physical system from the behavior predicted by the linear response theory. In particular a saturation of the power localized at $`f_s`$ is detected for large values of the amplitude of the modulating signal.

The second investigation concerns the frequency dependence of $`|M_1|`$. Within the framework of the linear response theory, at fixed values of $`V_n`$, $`|M_1|`$ is related to $`V_s`$, $`f_s`$ and $`V_n`$ through the relation

$$|M_1|<v_d^2>\frac{V_s}{V_{rms}^2}\frac{\lambda _{min}}{(\lambda _{min}^2+(2\pi f_s)^2)^{1/2}},$$
(5)
where $`<v_d^2>`$ denotes a stationary mean value of the unperturbed system and $`\lambda _{min}`$ is the smallest non-vanishing eigenvalue of the Fokker-Plank operator of the system without periodic driving . This quantity is an exponential function of the noise amplitude $`V_{rms}`$ under the hypothesis of white noise. We set $`V_s`$=0.067 V to ensure that we are in a region of parameters where the linear response theory may apply and we perform our experiments as a function of $`f_s`$ for various values of $`V_{rms}`$. In Fig. 2 we show the results obtained. The general trend predicted by Eq. (4) is observed. When $`f_s<<\lambda _{min}/2\pi `$ a constant value of $`|M_1|`$ is detected whereas in the opposite regime the value of $`|M_1|`$ is decreasing as a function of the frequency. Concerning the functional form of $`|M_1|`$, we observe that $`|M_1|f_s^{1.3}`$. This is close but not coincident with the behavior expected from the linear response theory $`|M_1|f_s^1`$. The observed deviation might be ascribed to one or more than one of the following possibilities: (i) a distortion introduced by the noise background present in our measurements; (ii) an additional frequency dependence which is present through the term $`<v_d^2>`$ of Eq. (4)and (iii) the colored nature of the noise $`v_n(t)`$.

One key aspect of the SR phenomenon is the statistical synchronization that takes place when the Kramers time $`T_K(V_{rms})`$ between two noise induced inter-well transition is of the order of half period of the periodic forcing. In other words statistical synchronization occurs when $`f_s=1/2T_K(V_{rms})`$. In our measurement, we verify the validity of this description with the following procedure. We set $`V_s=0.067`$ V and we measure the SNR of the output signal $`v_d(t)`$ at the frequency of the modulating signal for six values of $`f_s`$ ranging from 1 to 105 Hz. We use these experimental results to single out for which value of $`V_{rms}V_{rms}^{}(f_s)`$ a maximum of the SNR is detected. This is of course the state of maximal statistical noise induced synchronization. We then compare $`V_{rms}^{}(f_s)`$ with the function $`y(V_{rms})=r_K(V_{rms})/2`$, where the Kramers rate $`r_K(V_{rms})`$ is measured in the absence of a modulating signal. The results are shown in Fig. 3. From the figure it is clear that statistical synchronization is observed for all the investigated frequencies, supporting the traditional interpretation  of the stochastic resonance mechanism over a frequency range of the modulating signal spanning six frequency orders of magnitude. It is worth pointing out that synchronization is observed in spite of the fact that our experiments are performed in a regime of moderately colored noise.

We now investigate the SR phenomenon by studying both the signal power amplification

$$\eta =\frac{P_1}{P_{in}}=4\left[\frac{|M_1|}{V_s}\right]^2$$
(6)
and the signal to noise ratio

$$SNR=10\mathrm{log}_{10}\left[\frac{P_1}{N_1}\right].$$
(7)
The $`SNR`$ is customary given in dB and it is obtained by dividing the output signal power level $`P_1`$ to the noise level signal $`N_1`$. Both quantities are measured at the frequency of the modulating signal.

We investigate both the effect of varying the amplitude and the frequency of the modulating signal on $`\eta `$ and $`SNR`$. We first consider the role of the frequency of the modulating signal. Specifically we investigate the SR phenomenon as a function of $`V_{rms}`$ by keeping $`V_s`$ constant (we choose $`V_s=0.067`$ V) and by varying $`f_s`$ from 1 to $`10^6`$ Hz. For each pair of the control parameters $`V_s`$ and $`f_s`$, we vary $`V_{rms}`$ from 0.67 to 5.33 V. The measured values of the power amplification $`\eta `$ are collected in Fig. 4, where we show $`\eta `$ as a function of $`V_{rms}`$ for 7 different values of $`f_s`$, which are 1, 10, 100, $`10^3`$, $`10^4`$, $`10^5`$ and $`10^6`$ Hz. The classical profile of the SR phenomenon  is observed for the lowest values of $`f_s`$. For higher values of $`f_s`$, $`\eta `$ deviates from the canonical SR profile by lowering and broadening its maximum. These results are in qualitative agreement with the explicit results theoretically obtained for the signal power amplification in a model bistable system .

The next investigation of the power amplification $`\eta `$ concerns the study performed by keeping $`f_s`$ constant whereas $`V_s`$ is varied. We set $`f_s=10`$ Hz and we vary $`V_s`$ from 0.0067 to 1.00 V. For all the selected values of $`V_s`$ we check that the amount of the amplitude of the modulating signal is not sufficient to induce deterministic jumps between the two wells. In Fig. 5 we show the experimental values of $`\eta `$ obtained for 8 different values of $`V_s`$. The selected values are 0.0067, 0.017, 0.033, 0.067, 0.167, 0.333, 0.667 and 1.00 V. In the figure the top line corresponds to $`V_s=0.0067`$ V whereas the bottom line refers to the value $`V_s=1.00`$ V. The profile of $`\eta `$ becomes progressively more sharp around the value $`V_{rms}1.5`$ V for decreasing values of $`V_s`$. This experimental finding is in qualitative agreement with the results of theoretical calculations of Ref. . In the figure, the curves measured for highest values of $`V_s`$ tend to collapse into a unique curve that the theory indicates as the limit behavior predicted by the linear response theory for negligible values of the amplitude of the modulating signal. On the other hand, a difference is detected when one considers the lowest values of $`V_s`$. In these cases we experimentally detect a deviation from the expected limit curve. These deviations are observed in our experiment because for these values of $`V_s`$ (0.0067, 0.017 and 0.033 V) the signal $`P_1`$ becomes of the same level of the noise level $`N_1`$ and it is therefore indistinguishable from it. One can verify quickly the above statements by inspecting Fig. 6 where we present the $`SNR`$ measured under the same conditions of Fig. 5. In Fig. 6 the bottom curve refers to the case $`V_s=0.0067`$ V whereas the top curve is obtained by setting $`V_s=1.00`$ V. From the figure it is evident that for $`V_s`$ equal to 0.0067, 0.017 and 0.033 V, the $`SNR`$ becomes zero within the experimental errors for a wide range of values of $`V_{rms}`$. This effect is reflected into the deviation of $`\eta `$ from the limit curve observed in Fig. 5.

A simultaneous inspection of Figs 5 and 6 show that the results obtained with the highest values of $`V_s`$ are associated with high values of the $`SNR`$ but are at the same time seriously affected by nonlinear distortion. This nonlinear distortion manifests itself in the broadening of the SNR curve. In other words, by using high values of $`V_s`$ it is possible to detect a wide interval of $`V_{rms}`$ where the SR phenomenon occurs, however this interval is not well described in terms of linear response theory. On the other hand by using low values of $`V_s`$ one observe experimentally SR on a more limited interval of $`V_{rms}`$ but the experimental results are in this interval well described by a linear response theory. Hence from an experimental point of view the more straightforward investigation of the SR phenomenon requires the selection of a value of the amplitude of the modulating signal which allows the detection of a large but undistorted signal. In our case this condition is attained when $`V_s0.067`$ V.

We also investigate the detailed behavior of the output noise power level $`N_1`$ at the frequency of the modulating signal. In Fig. 7 we show $`N_1`$ as a function of $`V_{rms}`$ for several values of $`V_s`$ ranging from $`0.0067`$ to $`1.00`$ V. In these investigation $`f_s`$ is kept constant at the value of 10 Hz. We observe that the noise power $`N_1`$ slightly decreases in an interval of values of $`V_{rms}`$ by increasing $`V_s`$. The noise level $`N_1`$ sharply increases at the onset of the SR phenomenon, reaches a maximum and then decreases. Depending on the value of $`V_s`$, the noise level may decrease monotonically to the asymptotic value observed for high values of $`V_{rms}`$ or reach a minimum value and then increases until reaching the same asymptotic value. In other words, we detect a dip in the noise level for a finite value of $`V_{rms}`$ for high values of $`V_s`$. The dip is shown in the inset of Fig. 7 for the measurements done by setting $`f_s=10`$ Hz. The noise dip is more pronounced for high values of the signal amplitude. A similar behavior is also observed for values of $`f_s`$ satisfying the condition $`f_s<1`$ kHz. By taking into account the results previously obtained concerning the deviation form the behavior predicted in terms of linear response theory, we conclude that this behavior is belonging to the nonlinear response regime of stochastic resonance. In this regime, sometimes called weak-noise limit, the amplitude of the periodic signal can not be regarded as weak with respect to the noise intensity ($`V_sV_{rms}`$) and the linear-response theory or the perturbation theory is no longer valid .

In the nonlinear response regime we experimentally detect the phenomenon of phase and frequency locking. Specifically by increasing the value of the amplitude of the modulating signal one observes jumps between the two stable states occurring at phases which are progressively more synchronized with the phases of the modulating signal. Moreover a locking between the period of the output signal and the period of the input signal is also observed for given values of the noise amplitude. We address this last phenomenon as frequency locking. An example of phases and frequency locking is shown in Fig. 8 where we show the digitized time series of $`v_d(t)`$ and $`v_s(t)`$ recorded by setting $`f_s=10`$ Hz, $`V_s=0.667`$ V and $`V_{rms}=1.33,2.00`$ and $`4.67`$ V for the top, middle up and middle down time series of Fig. 8. The bottom time series of Fig. 8 shows the time evolution of the modulating signal. All the time series are digitized synchronously. By inspecting Fig. 8 one notes that for low levels of noise amplitude ($`V_{rms}=1.33`$ V, top time series) the system jumps randomly from one state to the other but the jumps are statistically synchronized in phases. In fact they occur preferentially at times $`t=nT/2`$, where $`T`$ is the period of the modulating signal and $`n`$ is an integer. When the noise amplitude is increased ($`V_{rms}=2.00`$ V, middle top time series) we still observe a phase synchronization, but in this case jumps occurs preferentially at times $`t=nT/2+T/4`$. Moreover for the present value of the noise amplitude jumps occurs with probability almost one at each period. This means that in addition to the phase synchronization we also observe frequency synchronization. It is worth pointing out that the noise amplitude at which we observe phase and frequency synchronization is always detected in the region of the dip observed in the noise level of Fig. 7. In other words the dip of the noise may be interpreted as a manifestation of the fact that phases and frequency locking are simultaneously present. By increasing the noise amplitude ($`V_{rms}=4.67`$ V, middle bottom time series) jumps becomes very frequent inside a single period of the modulation signal so that phase and frequency synchronization is progressively lost. A similar effect has been theoretically considered in the literature  recently.

## IV Conclusions

We report an experimental study of stochastic resonance in a physical system. Our physical system, which is characterized by versatility and high stability, allows us to investigate with high precision the SR phenomenon in a wide range of parameters, such as the frequency and the amplitude of the modulating signal and the noise amplitude. In the experiments presented here, the frequency range is spanning up to seven orders of magnitude whereas the amplitude range spans more than two orders of magnitude.

Theoretical and experimental investigations have been mainly focused on the linear regime of SR. However for a complete description of the SR phenomenon it is also important to investigate the nonlinear response regime of SR. We experimentally investigate the degree of consistence of our experimental results with the results expected in terms of the linear response theory. We find a range of experimental parameters within which the linear response theory describes quite well the investigated dynamics. However, outside these intervals, nonlinear deviations from the prediction of the linear response theory are clearly detected. These deviations primarily manifest them-self (i) in a saturation of the output power spectral density signal and of the signal amplification and (ii) in a non-monotonic behavior of the output noise level associated with a high degree of phase and frequency synchronization.

We wish to thank ASI, INFM and MURST for financial support.


text:
# Single and Paired Point Defects in a 2D Wigner Crystal

## Acknowledgments

PP acknowledges stimulating discussions with Antonio Castro-Neto and partial funding from DMR98-12422. LC acknowledges useful discussions and computational help with Dr. G. Bauer and Dr. B. Militzer and support by Fundacão de Amparo à Pesquisa do Estado de São Paulo (FAPESP). DC is supported by DMR98-02373. Computations were performed at National Computational Science Alliance (NCSA).


text:
# Resonances in Fock Space: Optimization of a SASER device

## Abstract

We model the Fock space for the electronic resonant tunneling through a double barrier including the coherent effects of the electron-phonon interaction. The geometry is optimized to achieve the maximal optical phonon emission required by a SASER (ultrasound emitter) device. PACS numbers: 73.20.Dx, 73.40.Gk, 73.50.Rb

, and

The possibility of generating coherent phonons in a double barrier semiconductor heterostructure was first proposed a few years ago. This is the basis of a SASER device which transforms the electric potential energy in a single vibrational mode of the lattice. This is facilitated by the electronic confinement in a double barrier structure. The phonon emission appears when the energy of the resonant state is one quantum $`\mathrm{}\omega _0`$(LO phonon energy) bellow the energy of the incoming electrons. As in laser devices this is enhanced if the first excited state of the well lies bellow the Fermi energy and becomes overpopulated. According to ref. the emitted LO phonons decay coherently into a pair of $`LO`$ and TA phonons the last the useful ones in a SASER device.

In this paper we want to explore the case in which well’s ground state mediates the decay of the emitter states into the collector’s ones plus a phonon. This feature represents a resonance in the electron-phonon Fock space and is observed as a satellite peak in the current. This resonant condition is tuned directly by the applied voltage and we expect that its optimization could also provide enough emission of primary phonons to allow for SASER operation. We carry out the modeling of the electronic structure and the electron-phonon interaction to get a minimal structure in the Fock space. Thus, the optimization of the phonon emission for different geometries of the device (height and width of the barriers, field intensity) can be discussed in simple terms.

We consider a one-dimensional model for a double barrier including the interaction with LO phonons in the well, neglecting the effects derived from the accumulated charge. This will give results comparable to the 3-D case when $`\epsilon _F`$ is small, thus limiting the number of traversal modes; or in the presence of a high magnetic field perpendicular to the plane of the barriers which quantize these modes in Landau levels.We do not consider the phonon-phonon interaction that leads to the decay of the LO phonons.

The Hamiltonian is a sum of an electronic contribution, a phonon contribution and an electron-phonon interaction term.

$$=_e+_p+_{ep}$$

$$_e=\underset{j}{}E_jc_j^+c_j\underset{j,k}{}V_{j,k}(c_j^+c_k+c_k^+c_j),$$

$$_p=\mathrm{}\omega _0\underset{k[\mathrm{well}]}{}b_k^+b_k,\mathrm{and}_{ep}=V_g\underset{k[\mathrm{well}]}{}c_k^+c_k(b_k^++b_k).$$

where $`c_j^+`$and $`c_j`$ are electron operators on site $`j`$, $`E_j`$ is the site dependent diagonal energy and $`V_{j,k}=V\delta _{j\pm 1,k}`$ are the hopping parameters. We assume that the potential drop $`\mathrm{𝐞𝐕}`$ is linear through the double barrier and limited to it. NL and NR are the number of sites in the left and right barriers and Nw are those in the well, the associated lengths are Li=Ni 2.825Å. There is a single well state in the energy range of interest.

Since the most important interaction between electrons and phonons in polar semiconductors involves longitudinal optical (LO) phonons, only one phonon mode with frequency $`\omega _0`$ is considered. The electron-phonon interaction is limited to the well region and the coupling to the phonons is denoted with $`V_g`$. The model is represented schematically in figure 1.

For simplicity we restrict the problem to the case in which we have either $`0`$ or $`1`$ phonons with no phonons in the well before the scattering process. By modifying $`V_gV_g\sqrt{n+1}`$ this also represents a finite temperature emission $`nn+1.`$ The effective mass is taken to be $`0.067`$ me, the LO phonon frequency $`\mathrm{}\omega _0=36`$ meV and the value of the hopping parameter $`V=7.1018`$eV. $`V_g`$is taken $`10`$meV which gives a typical electron-phonon interaction strength $`g=(V_g/\mathrm{}\omega _0)^20.1`$. The barrier heights are $`300`$meV and the Fermi energy $`\epsilon _F`$ is taken between $`10`$ and $`20`$ meV.

This discrete model is solved exactly using a decimation procedure for the sites in the barriers and the well . The leads are taken into account by adding a proper self-energy. The transmittances are computed from the Green’s functions for the system .

Let us denote with $`T_{0,0}^{RL}`$ ($`T_{0,0}^{LR}`$) and $`T_{1,0}^{RL}`$ ($`T_{1,0}^{RL}`$) the transmission coefficients from left (right) to right (left) where the subscripts $`0`$ and $`1`$ denote the number of phonons in the outgoing (first subscript) and in the incoming (second subscript) channel. The total current is a sum of an elastic current $`I_{\mathrm{el}}`$and an inelastic current $`I_{\mathrm{in}}`$ (with the emission of one phonon during the scattering process). These currents can be calculated from the following expressions

$`I_{el}`$ $`=`$ $`(2e/h){\displaystyle [T_{0,0}^{RL}f_L(\epsilon )T_{0,0}^{LR}f_R(\epsilon )]𝑑\epsilon },`$
$`I_{in}`$ $`=`$ $`(2e/h){\displaystyle [T_{1,0}^{RL}f_L(\epsilon )T_{1,0}^{LR}f_R(\epsilon )]𝑑\epsilon };`$
where $`f_L(\epsilon )`$ and $`f_R(\epsilon )`$ are the Fermi functions for the left and right leads.

For a given configuration of the system a curve of inelastic current vs. applied bias is obtained and its maximum value $`I_{\mathrm{in}}^{\mathrm{max}}`$ can be extracted. Figure 2 shows $`I_{\mathrm{in}}`$-V curves as we change NR. The peaks in these curves correspond to the inelastic contributions to the main peak and to the satellite peak in the total current respectively. This figure also shows that the peaks are shifted to higher voltages as NR is increased. This shows a strong renormalization of the resonant energies due to the electrodes. In figure 3 we present $`I_{\mathrm{in}}^{\mathrm{max}}`$ vs. NR curves for different values of NL. These curves exhibit a maximum for $`I_{\mathrm{in}}^{\mathrm{max}}`$ as a function of NR. The optimal configurations correspond to asymmetric structures with wider right barriers. This can be understood by means of the following argument. Increasing the lifetime of the electrons in the well favors the electron-phonon interaction and thus increases the inelastic current. This can be done by choosing wider (or higher) barriers. In spite of this, as an effect of the asymmetry produced by the applied bias, the lifetime is still controlled mainly by the right barrier. On the other hand increasing the length of the barriers increases the reflectivity of the device diminishing the currents, here it is the left barrier which plays the main role. Then there is a trade off between these two effects that maximizes the phonon emission.

In summary, we have used a simple model to show that the asymmetry in double barrier structures plays an important role in the $`I_{\mathrm{in}}`$-V characteristics and to predict how it can be controlled to optimize LO phonon emission. In particular we show that the optimal configuration corresponds to a collector barrier with a length which doubles that of the emitter.

We acknowledge finantial support from CONICET, SeCyT-U.N.C., AnPCyT and an international grant from Andes-Vitae-Antorchas.


text:
# Simulating the effects of intergalactic grey dust

## 1. Introduction

Recent observations of Type IA Supernovae (SNe) at redshifts up to $`z0.8`$ (Riess et al. (1998); Perlmutter et al. 1999, hereafter P (99)) have made possible classical cosmological tests that require standard candles, such as the magnitude-redshift relation. The most dramatic result is that these SNe appear dimmer (by $`0.2`$ magnitudes) at high redshift than would be predicted in a non-accelerating universe, suggesting at face value that we live in an accelerating universe. However, other explanations are possible, including the one we consider here, namely that distant SNe are dimmer due to extinction by intergalactic dust.

As distant standard candles, SNe are sensitive probes of extinction in the intergalactic medium (IGM). The distribution of matter in the IGM may now be modelled accurately in the context of modern cosmology using hydrodynamic simulations (see e.g., Cen et al. (1994); Hernquist et al. (1996); Davé et al. (1999)). The resulting IGM is not smooth, but rather traces large-scale structure. If such structure contains not only dark matter and gas but also dust, this would result in significant variations in SNe brightnesses due to intervening extinction. Observationally, the distribution of Type Ia SNe magnitudes has a very small dispersion (P (99)). Thus by comparing simulations to the distribution of observed brightnesses, we can set limits on the amount of dust extinction and possibly constrain its spatial distribution with respect to intergalactic gas. In this Letter we present a technique for doing this, and apply it to the SNe observations of P (99).

Intergalactic grey dust has been examined in a series of papers by Aguirre (1999a,b; hereafter A (99)) and Aguirre & Haiman (1999), who develop a scenario in which small grains are preferentially destroyed during ejection from galaxies, polluting the IGM with large dust grains that are effectively grey in the bandpasses of the SNe data. This greyness is necessary in order not to violate tight limits on reddening from SNe data, which imply that galactic-type dust would provide negligible absorption (P (99)). Furthermore, a significant fraction of the dust must reside in the IGM. If the grey dust causing extinction were present only in the ISM of the supernova host galaxy, this would introduce too large a dispersion in observed SNe magnitudes (Riess et al. (1998)). In this study, we assume that grey dust blends smoothly from galaxies into the surrounding IGM. Our analysis is insensitive to the intrinsic properties of the dust, such as its opacity and grain size, since we use the simulations to directly translate the observed SNe magnitude distribution into a dust extinction in magnitudes. It is, however, sensitive to the way in which dust traces the distribution of gas in the IGM, and we will consider several simple but plausible variations of this relation.

In § 2 we describe our simulations of the IGM, and of dust extinction. In § 3 we describe our analysis method and results, including constraints on grey dust afforded by current SN observations. In § 4 we discuss systematic uncertainties, and the implications of our results.

## 2. Simulations of Intergalactic Dust

We employ a hydrodynamic simulation of a $`\mathrm{\Lambda }`$-dominated cold dark matter model, with $`\mathrm{\Omega }_m=0.4`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.6`$, $`\mathrm{\Omega }_b=0.02h^2`$, $`H_0=65\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, and $`\sigma _8=0.8`$. Our simulation volume is $`50h^1`$ Mpc with $`10h^1`$ kpc spatial resolution, having $`144^3`$ dark matter and $`144^3`$ gas particles, and was evolved from $`z=490`$ using Parallel TreeSPH (Davé, Dubinski & Hernquist (1997)).

In order to obtain dust column densities along lines of sight, we consider three different ways that dust may trace intergalactic hydrogen gas:

1. $`\rho _{\mathrm{dust}}\rho _{\mathrm{gas}}`$: Dust traces gas linearly.
2. $`\rho _{\mathrm{dust}}\rho _{\mathrm{metal}}`$: Dust traces metals linearly.
3. $`\rho _{\mathrm{dust}}\rho _{\mathrm{gas}}^2`$: Dust traces gas quadratically.

As our simulation makes no direct prediction for the metallicity of gas, we adopt a heuristic prescription (c.f. Cen & Ostriker (1999)), in the second case above. We assume that the metallicity is $`10^2`$ solar if the gas overdensity is less that 10, solar if the overdensity is greater than 1000, and log-linear in between.

Figure 1: Dust extinction maps in $`2.2^{}\times 2.2^{}`$ patches of sky out to $`z=0.05`$ (left) and $`z=0.5`$ (right), for (a) $`\rho _{\mathrm{dust}}\rho _{\mathrm{gas}}`$ and (b) $`\rho _{\mathrm{dust}}\rho _{\mathrm{gas}}^2`$ The median extinction magnitude to $`z=0.5`$ is set equal to $`0.4`$ for both (a) and (b). The maps correspond to what would be seen against a white background.

To extract dust extinction values from the simulations, we assume that the gas associated with each particle is spread over its SPH smoothing volume (see e.g., Hernquist & Katz (1989)). We perform a numerical integration of gas column density along 5000 rays cast through these volumes, at the same time applying one of the three transformations given above to relate gas to dust densities. To reach the required path lengths out to $`z0.5`$, we follow rays through 26 simulation volumes, each ray entering through a random point on a random face. This yields the column density of dust to $`z=0.5`$ along each line of sight.

In Figure 1 we show extinction maps of $`2.2^{}\times 2.2^{}`$ patches of sky, for (a) $`\rho _{\mathrm{dust}}\rho _{\mathrm{gas}}`$ and (b) $`\rho _{\mathrm{dust}}\rho _{\mathrm{gas}}^2`$. The median extinction to $`z=0.5`$ was set to be equal (to $`0.4`$ mag) for cases (a) and (b). In Figure 2, we show how the mean dust extinction and its dispersion varies with redshift. Here, the dispersion in all three panels was set to be the same small value, equal to the difference in quadrature of the dispersion in SNe magnitudes at high ($`\sigma _{z=0.5}=0.157`$) and low ($`\sigma _{z=0.05}=0.154`$) redshifts observed by P (99), namely 0.03 mag. Figs. 1 and 2 show that the mean extinction is much greater when the dust is more smoothly distributed. We will now quantitatively explore the constraints that can be put on the dust extinction by using the observed distribution of SNe magnitudes.

Figure 2: Mean and dispersion of the dust extinction as a function of redshift. Filled squares show mean, open circles with line show the dispersion. Values have been scaled to a dispersion of 0.03 magnitudes at $`z=0.5`$.

## 3. Comparison with observations

We make use of two characteristics of the observed SN data in our comparison, the change with redshift of the dispersion in SNe magnitudes, and the shape of the histogram of SNe magnitudes. As mentioned above, P (99) found little difference in the dispersions of two samples with $`\overline{z}0.05`$ and $`\overline{z}0.5`$. As they stated, this leaves little room for dispersion due to dust, as this dispersion is expected to increase for longer path lengths (see Fig. 2). In order to quantify this, and the effect of the distribution shape, we generate simulated SNe magnitudes, and compare them to the P (99) data using a maximum likelihood approach.

The observational datasets we use are both taken from P (99) (their Tables 1 and 2), being the high-$`z`$ SNe of the Supernova Cosmology Project, and the low-$`z`$ sample of Calán-Tololo SNe survey (Hamuy et al. (1996)). We use 40 (non-reddened SNe) of the former SNe, between $`z=0.172`$ and $`z=0.83`$ ($`\overline{z}0.5`$) , and 16 of the latter SNe which lie between $`z=0.02`$ and $`z=0.101`$ ($`\overline{z}0.05`$) .

We generate simulated datasets for each dust model described in § 2. For each dust model, we vary two parameters; first, $`M_C`$, a cosmological magnitude shift applied to all simulated SNe at a given $`z`$, normalized so that $`M_C=0`$ at $`z=0.5`$ corresponds to the best fitting cosmology found by P (99) (with $`\mathrm{\Omega }_m=0.28`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.72`$); and second, $`A_V`$, the median V-band magnitude of dust extinction out to $`z=0.5`$111 We use the median in order to be less sensitive to long tails of the distribution (c.f. Figure 3) when renormalizing extinction values.. $`M_C0.2`$ then corresponds to an open model with $`\mathrm{\Omega }_m0.3`$, and $`M_C0.4`$ to an $`\mathrm{\Omega }_m=1`$ model. We generate the simulated datasets as follows:
(1) We renormalize the 5000 simulated lines-of-sight so that the median dust extinction to $`z=0.5`$ equals $`A_V`$.
(2) We add the cosmological shift, $`M_C`$.
(3) We broaden the magnitude distribution, which involves convolving it separately with a Gaussian of width given by the observational error of each SN, and then adding these distributions together. When doing this, we include an “intrinsic” SN dispersion of $`\sigma _{\mathrm{int}}=0.17`$ mag (P (99)). Varying this (by $`\pm 0.1`$ mag) makes little difference to the results (see §4).
(4) We truncate the distribution at a magnitude difference $`\mathrm{\Delta }_M=1`$, to roughly account for the fact that SNe along lines of sight passing through high extinction regions would not make it into these magnitude-limited samples.

Figure 3: The PDF of magnitude differences, $`\mathrm{\Delta }_M`$, for 5000 simulated SNe. The shaded region is included in the calculation of the dispersion; non-shaded region indicates heavily obscured SNe that are not included in the PDF normalization. The histogram shows 40 of the high-$`z`$ SNe of P (99).

We derive the probability distribution function (PDF) of SNe magnitude differences $`\mathrm{\Delta }_M`$ predicted by the simulation, so that the predicted number of SNe between $`\mathrm{\Delta }_M`$ and $`\mathrm{\Delta }_M+d\mathrm{\Delta }_M`$ is $`NP(\mathrm{\Delta }_M)d\mathrm{\Delta }_M`$, where $`N`$ is the number of observed SNe. Figure 3 shows histograms of the PDF of $`\mathrm{\Delta }_M`$, for $`A_V=0.4`$ mag and $`M_C=0.4`$. We also plot the observational data of P (99). In effect, for this plot, we have brightened the simulated SNe to mimic an $`\mathrm{\Omega }_m=1`$ model and then dimmed them with dust. We can see that in panel (a), where the dust is fairly smoothly distributed, the simulated PDF is not too different from the observations. In the other panels, which have clumpier dust, there is a skewness not seen in the observed data, due to a long tail of high extinction lines-of-sight.

For each set of simulated lines of sight, we form the relative likelihood of drawing all the observed SNe magnitudes, $`=_i^NP(\mathrm{\Delta }_{M,i})`$, where $`\mathrm{\Delta }_{M,i}`$ is the magnitude difference of SN $`i`$. We define the quantity $`S=2\mathrm{ln}`$ and assume that $`S`$ follows a $`\chi ^2`$ distribution in order to derive confidence limits on the parameters $`M_C`$ and $`A_v`$. In the present analysis, we use results at two different redshifts ($`z=0.05`$ and $`z=0.5`$), and combine the two by adding the values of $`S`$. It is this step, combining the likelihoods at two different redshifts, which constrains the amount of additional dispersion (or the change in the shape of the magnitude distribution) due to grey dust between low and high redshifts. Contours of $`\mathrm{\Delta }S=2.3,6.2`$, and $`11.8`$ (the difference in $`S`$ from its minimum value), representing $`1,2`$ and $`3\sigma `$ intervals of joint confidence, are plotted in Figure 4.

Figure 4: 1, 2 and 3$`\sigma `$ contours for the cosmological magnitude shift, $`M_C`$, with respect to a model with $`\mathrm{\Omega }_\mathrm{\Lambda }=0.72`$ and $`\mathrm{\Omega }_M=0.28`$, and the median extinction due to grey dust ($`A_V`$).

Figure 4 shows that the smoothest distribution of dust easily accomodates enough extinction to reconcile a non-accelerating ($`M_C=0.2`$) or flat ($`M_C=0.4`$) universe, as there is a strong degeneracy between the cosmological magnitude shift and dust extinction. The other cases, where dust traces metals and $`\rho _{\mathrm{gas}}^2`$, are coincidentally quite similar. Such a distribution of dust would be mildly disfavored in a non-accelerating open universe and ruled out at $`99\%`$ confidence in an Einstein-de-Sitter universe. The constraints arise mostly due to the shape of the distribution; the fact that the observed dispersions are similar at $`z=0.05`$ and $`z=0.5`$ is relatively unimportant, only making a noticeable difference in panel (a) where the shape is similar to the observed distribution.

## 4. Discussion

We have found that current SNe datasets appear to have some power to constrain grey dust models. We find it somewhat unlikely (at the $`1.52`$ level) that sufficient grey dust could be distributed in the relatively clumpy fashion expected for the metal distribution (c.f., Gnedin & Ostriker (1997), Davé et al. (1998), Cen & Ostriker (1999)) to reconcile the P (99) data with a non-accelerating universe. This suggests that any substantial grey dust component must be largely segregated from the galaxies where it was formed. In the specific grey dust model of A (99), sputtering does destroy dust more effectively in denser regions, but once in the IGM the large grains have long lifetimes, so that dust is still likely to trace the metals. The question of how the dust and metals are distributed can be answered self-consistently by modelling the relevant physical processes (dust and metal ejection from galaxies) directly in the simulations (e.g., Aguirre et al., in preparation).

There are many alternative explanations for SNe appearing dimmer in the past. For example, there may be metallicity effects in the host galaxy (Höflich et al. (1999)), intrinsic evolution (Riess et al. (1999)), observational selection differences between CCDs used for distant SNe and photographic plates used for nearby samples (Howell, Wang & Wheeler (1999)), or time evolution of the gravitational constant (Amendola, Corasaniti & Occhionero (1999); Garcia-Berro et al. (1999)). Such effects would alter the interpretation of our parameter $`M_C`$. Gravitational lensing magnification (Metcalf (1999)) would also change the magnitude dispersion at higher redshifts.

The study we have presented is reasonably general, in that our analysis is not dependent on the (unknown) microscopic properties of this hypothetical intergalactic dust. Still, there are some possible systematic uncertainties, which we now consider. We assume (as do P (99)) that the observed dispersion in SNe magnitudes in excess of the estimated observational errors is an intrinsic property of SNe and does not vary with redshift. One could envision scenarios in which the intrinsic dispersion is lower at high redshift, thus allowing more dispersion from dust. However, as mentioned earlier, most of the statistical power of our analysis comes from the shape of the distributions, which is not significantly affected by a change in the intrinsic dispersion. Also, we find that when we decrease the intrinsic SN dispersion ($`\sigma _{\mathrm{int}}`$), models with more dust fit slightly better at high-$`z`$, but the low-$`z`$ fit becomes worse, a trade-off which means that the overall results hardly change. The distribution shape depends on our assumption that the intrinsic dispersion has a Gaussian distribution in magnitudes, whereas a distribution skewed to fainter magnitudes could weaken our constraints. With a larger sample of low-z SNe we might be in a position to test this, by using the distribution of low-z SN magnitudes to make a simulated high-z sample with the correct intrinsic distribution shape. Also, we have made a simplifying approximation in our simulated high-z samples, by using dust extinction magnitudes that result from integrating the dust contribution from $`z=0`$ to exactly $`z=0.5`$. The real SNe are at a range of redshifts; our constraints are effectively conservative because of this. We have also assumed that the total mass of dust does not change with redshift from $`z=0.5`$ to $`z=0`$; dust increasing with time would strengthen our constraints, while a decrease seems implausible. Finally, the truncation of our simulated magnitude distribution at $`+1`$ mag is a rather approximate procedure. We have tried changing the cutoff, and find that with none (an unrealistic case) the constraints become stronger, as we might expect. With a lower cutoff, the right side of the observed distribution is not reproduced. One additional degree of freedom would involve changing the functional form of the cutoff. In the future, as observations improve, we plan to simulate such observational selection effects more carefully.

On the simulation side, our relatively low mass resolution ($`m_{\mathrm{gas}}=8.5\times 10^8M_{}`$) means that we miss fluctuations in the extinction that occur on smaller mass scales. If we had higher resolution, this should have the effect of making our dust constraints tighter. We find that the dispersion in projected magnitudes largely depends on small scale fluctuations. We tested this by splitting the simulation into small sub-volumes and then shuffling them to remove large scale correlations (see e.g., Bhavsar & Ling (1988)) before projection, and we obtained similar results. Changing the assumed cosmological model could have some effects, although it is difficult to conceive of models which have less power on small scales while fitting other constraints (see e.g., White & Croft (2000)).

The models of A (99) have grains that do produce some reddening, typically in the infrared. A complementary approach to ours is therefore to use information from different color bands. Such an approach was used by P (99), for example, to constrain normal galactic-type dust. Recent near-IR observations by Riess et al. (2000) have made reddening constraints tight even for non-standard large dust grains. Our approach does not make any use of color information, and so constrains the most extreme scenario in which the dust is totally “grey”.

In summary, we have used cosmological hydrodynamic simulations to explore how intergalactic grey dust could affect observations of high redshift supernovae, and how supernova data can constrain grey dust extinction. We conclude that only a fairly smooth distribution of dust could readily mimic the effect of an accelerating Universe. Such a distribution would be strange, as the dust would be strongly biased away from metal producing regions. More realistic dust distributions are mildly disfavored, but upcoming samples of SNe (at current redshifts) should enable us to put tighter limits as they more precisely determining the shape of the brightness distribution. At higher redshifts ($`z1`$), grey dust predicts that SNe should show increased dimming, while decreased dimming would be strong evidence for a cosmological constant.

We thank Anthony Aguirre, Bruce Draine, and Bob Kirshner for useful discussions, and David Weinberg for helpful comments on the manuscript.


text:
# Time Dependent Theory for Random Lasers

The interplay of localization and amplification is an old and interesting topic in physics research. With promising properties, mirror-less random laser systems are widely studied both experimentally and theoretically. Recently new observations  of laser-like emission were reported and showed new interesting properties of amplifying media with strong randomness. First, sharp lasing peaks appear when the gain or the length of the system is over a well defined threshold value. Although a drastic spectral narrowing has been previously observed , discrete lasing modes were missing. Second, more peaks appear when the gain or the system size further increases over the threshold. Third, the spectra of the lasing system is direction dependent, not isotropic. To fully explain such an unusual behavior of stimulated emission in random systems with gain, we are in need of new theoretical ideas.

Theoretically, a lot of methods have been used to discuss the properties of such random lasing systems. Based on the time-dependent diffusion equation, earlier work of Letokhov  predicted the possibility of lasing in a random system and Zyuzin  discussed the fluctuation properties near the lasing threshold. Recently, John and Pang  studied the random lasing system by combining the electron number equations of energy levels with the diffusion equation. Such a consideration predicted a reduction in the threshold gain for a laser action due to the increased optical path from diffusion. It also verified the narrowing of the output spectrum when approaching the gain threshold. By using the diffusion approach is not possible to explain the lasing peaks observed in the recent experiments in both semiconductor powders and in organic materials. The diffusive description of photon transport in gain media neglects the phase coherence of the wave, so it gives limited information for the wave propagation in the gain media. Another approach which is based on the time-independent wave equations for the random gain media can go beyond the diffusive description . But as was shown recently , the time-independent method is only useful in determining the lasing threshold. When the gain or the length of system is larger than the threshold value, the time-independent description will give a totally unphysical picture for such a system. To fully understand the random lasing system, we have to deal with time-dependent wave equations in random systems. On the other hand, the laser physics community  has developed some phenomenological theories to deal with gain media which were overlooked by the researchers working on random systems.

In this paper we introduce a model by combining these semi-classical laser theories with Maxwell equations. By incorporating a well-established FDTD ( finite-difference time-domain)  method we calculate the wave propagation in random media with gain. Because this model couples electronic number equations at different levels with field equations, the amplification is nonlinear and saturated, so stable state solutions can be obtained after a long relaxation time. The advantages of this FDTD model are obvious, since one can follow the evolution of the electric field and electron numbers inside the system. From the field distribution inside the system, one can clearly distinguish the localized modes from the extended ones. One can also examine the time dependence of the electric field inside and just outside the system. Then after Fourier transformation, the emission spectra and the modes inside the system can be obtained.

Our system is essentially a one-dimensional simplification of the real experiments . It consists of many dielectric layers of real dielectric constant of fixed thickness, sandwiched between two surfaces, with the spacing between the dielectric layers filled with gain media (such as the solution of dye molecules). The distance between the neighboring dielectric layers is assumed to be a random variable. The overall length of the system is $`L`$.

Our results can be summarized as follows: ($`i`$) As expected for periodic and $`short`$ ($`L<\xi `$, $`\xi `$ is the localization length) random system, an extended mode dominates the field and the spectra. ($`ii`$) For either strong disorder or the long ($`L\xi `$) system, we obtain a low threshold value for lasing. By increasing the length or the gain (higher gain can be achieved by increasing the pumping intensity) more peaks appear in the spectra. By examining the field distribution inside the system, one can clearly see that these lasing peaks are coming from localized modes. ($`iii`$) When the gain or the pumping intensity increases even further, the number of lasing modes do not increase further, but saturate to a constant value, which is proportional to the length of system for a given randomness. And ($`i`$v) the emission spectra are not same for different output directions which show that the emission is not isotropic. These findings are in agreement with recent experiments  and also make new predictions.

The binary layers of the system are made of dielectric materials with dielectric constant of $`\epsilon _1=\epsilon _0`$ and $`\epsilon _2=4\times \epsilon _0`$ respectively. The thickness of the first layer, which simulates the gain medium, is a random variable $`a_n=a_0(1+W\gamma )`$ where $`a_0=300`$nm, $`W`$ is the strength of randomness and $`\gamma `$ is a random value in the range \[-0.5, 0.5\]. The thickness of second layer, which simulates the scatterers, is a constant $`b=180`$nm. In the first layer, there is a four-level electronic material mixed inside. An external mechanism pumps electrons from ground level ($`N_0`$) to third level ($`N_3`$) at certain pumping rate $`P_r`$, which is proportional to the pumping intensity in experiments. After a short lifetime $`\tau _{32}`$, electrons can non-radiative transfer to the second level ($`N_2`$). The second level ($`N_2`$) and the first level ($`N_1`$) are called the upper and the lower lasing levels. Electrons can be transfered from the upper to the lower level by both spontaneous and stimulated emission. At last, electrons can non-radiative transfer from the first level ($`N_1`$) back to the ground level ($`N_0`$). The lifetimes and energies of upper and lower lasing levels are $`\tau _{21}`$, $`E_2`$ and $`\tau _{10}`$, $`E_1`$ respectively. The center frequency of radiation is $`\omega _a=(E_2E_1)/\mathrm{}`$ which is chosen to be equal to $`2\pi \times 610^{14}`$ $`Hz`$ ($`\lambda =499.7`$ $`nm`$). Based on real materials , the parameters $`\tau _{32}`$, $`\tau _{21}`$ and $`\tau _{10}`$ are chosen to be $`1\times 10^{13}`$s, $`1\times 10^9`$s and $`1\times 10^{11}`$s. The total electron density $`N_0^0=N_0+N_1+N_2+N_3`$ and the pump rate $`P_r`$ are the controlled variables according to the experiments .

The time-dependent Maxwell equations are given by $`\times 𝐄=𝐁/𝐭`$ and $`\times 𝐇=\epsilon 𝐄/𝐭+𝐏/𝐭`$, where $`𝐁=\mu 𝐇`$ and $`𝐏`$ is the electric polarization density from which the amplification or gain can be obtained. Following the single electron case, one can show  that the polarization density in the presence of an electric field obeys the following equation of motion:

$$\frac{d^2P(t)}{dt^2}+\mathrm{\Delta }\omega _a\frac{dP(t)}{dt}+\omega _a^2P(t)=\frac{\gamma _r}{\gamma _c}\frac{e^2}{m}\mathrm{\Delta }N(t)E(t)$$
(1)
where $`\mathrm{\Delta }\omega _a=1/\tau _{21}+2/T_2`$ is the full width at half maximum linewidth of the atomic transition. $`T_2`$ is the mean time between dephasing events which is taken to be $`2.18\times 10^{14}`$s, $`\mathrm{\Delta }N(t)=N_1N_2`$ and $`\gamma _r=1/\tau _{21}`$ is the real decay rate of the second level and $`\gamma _c=\frac{e^2}{m}\frac{\omega _{a}^{}{}_{}{}^{2}}{6\pi \epsilon _0c^3}`$ is the classical rate. It is easy to derive  from Eq. (1) that the amplification line shape is Lorentzian and homogeneously broadened. Eq. (1) can be thought as a quantum mechanically correct equation for the induced polarization density $`P(t)`$ in a real atomic system.

The equations giving the number of electrons on every level can be expressed as follows:

$`{\displaystyle \frac{dN_3(t)}{dt}}`$ $`=`$ $`P_rN_0(t){\displaystyle \frac{N_3(t)}{\tau _{32}}}`$ (2)
$`{\displaystyle \frac{dN_2(t)}{dt}}`$ $`=`$ $`{\displaystyle \frac{N_3(t)}{\tau _{32}}}+{\displaystyle \frac{1}{\mathrm{}\omega _a}}E(t){\displaystyle \frac{dP(t)}{dt}}{\displaystyle \frac{N_2(t)}{\tau _{21}}}`$ (3)
$`{\displaystyle \frac{dN_1(t)}{dt}}`$ $`=`$ $`{\displaystyle \frac{N_2(t)}{\tau _{21}}}{\displaystyle \frac{1}{\mathrm{}\omega _a}}E(t){\displaystyle \frac{dP(t)}{dt}}{\displaystyle \frac{N_1(t)}{\tau _{10}}}`$ (5)
$`{\displaystyle \frac{dN_0(t)}{dt}}`$ $`=`$ $`{\displaystyle \frac{N_1(t)}{\tau _{10}}}P_rN_0(t)`$ (6)
where $`\frac{1}{\mathrm{}\omega }E(t)\frac{dP(t)}{dt}`$ is the induced radiation rate from level 2 to level 1 or excitation rate from level 1 to level 2 depending on its sign.

To excite the system, we must introduce sources into the system. To simulate the real laser system, we introduce sources homogeneously distributed in the system to simulate the spontaneous emission. We make sure that the distance between the two sources $`L_s`$ is smaller than the localization length $`\xi `$. Each source generates waves of a Lorentzian frequency distribution centered around $`\omega _a`$, with its amplitude depending on $`N_2`$. In real lasers, the spontaneous emission is the most fundamental noise  but generally submerged in other technical noises which are much larger. In our system, the simulated spontaneous emission is the $`only`$ noise present, and is treated self-consistently. This is the reason for the small background in the emission spectra shown below.

There are two leads, both with width of 3000 nm, at the right and the left sides of the system and at the end of the leads we use the Liao method  to impose an absorbing-boundary conditions (ABC). In the FDTD calculation, discrete time step and space steps are chosen to be $`10^{17}`$s and $`10^9`$m respectively. Based on the previous time steps we can calculate the next time step ($`n+1`$ step) values. First we obtain the $`n+1`$ time step of the electric polarization density $`P`$ by using Eq. (1), then the $`n+1`$ step of the electric and magnetic fields are obtained by Maxwell’s equations and at last the $`n+1`$ step of the electron numbers at each level are calculated by Eq. (2). The initial state is that all electrons are on the ground state, so there is no field, no polarization and no spontaneous emission. Then the electrons are pumped and the system begins to evolve according to equations.

We have performed numerical simulations for periodic and random systems. First, for all the systems, a well defined lasing threshold exists. As expected, when the randomness becomes stronger, the threshold intensity decreases because localization effects make the paths of waves propagating inside the gain medium much longer.

For a periodic or $`short`$ ($`L<\xi `$) random system, generally only one mode dominates the whole system even if the gain increases far above the threshold. This is due to the fact that the first mode can extend in the whole system, and its strong electric field can force almost all the electrons of the upper level $`N_2`$ to jump down to the $`N_1`$ level quickly by stimulated emission. This leaves very few upper electrons for stimulated emission of the other modes. In other words, all the other modes are suppressed by the first lasing mode even though their threshold values are only a little bit smaller than the first one. This phenomenon also exists in common lasers .

For $`long`$ ($`L>>\xi `$) random systems, richer behavior is observed. First we find that all the lasing modes are localized and stable around their localization centers after a long time. Each mode has its own specific frequency and corresponds to a peak in the spectrum inside the system. When the gain increases beyond the threshold, the electric field pattern (see Fig. 1a) shows that more localized lasing modes appear in the system and the spectrum inside the system (see Fig. 1b) gives more sharp peaks just as observed in the experiments . This is clearly seen in Figs 1a and 1b for a 80 cell random system, above threshold. In Figs. 1c and 1d similar values are shown for the 160 cell system. Notice that both the number of localized modes (Fig. 1c) of the field, as well the number of lasing peaks (Fig. 1d) are larger now. The exact position of the lasing peaks depends on the random configuration. Notice that the lasing peaks are much narrower than the experimental ones . This is due to the 1d nature of our model. In the present case only two escaping channels exist, so it’s more difficult for the wave to get out from the system which has a higher quality factor. When the gain is really big, we find the number of lasing modes will not increase any more, so a saturated number $`N_m`$ of lasing modes exists for the long random system. This is clearly seen in Fig. 2, where we plot the number of modes $`N_m`$ vs the pumping rate $`P_r`$. In Fig. 3, we plot the spectral intensity vs the wavelength for different input transitions (or equivalently pumping rates). Notice these results are in qualitative agreement with the experimental results shown in Fig. 2 of the paper of Cao et. al. .

These multi-lasing peaks and the saturated-mode-number phenomena are due to the interplay between localization and amplification. Localization makes the lasing mode strong around its localization center and exponentially small away from its center so that it only suppress the modes in this area by reducing $`N_2`$. When a mode lases, only those modes which are $`far`$ $`enough`$ from this mode can lase afterwards. So more than one mode can appear for a long system and each mode seems to $`repel`$ each other. Because every lasing mode dominates a certain area and is separated from other modes, only a limited number of lasing modes can exist for a finite long system even in the case of large amplification. We therefore expect that the number of surviving lasing modes $`N_m`$ should be proportional to the length of the system $`L`$ when the amplification is very large. Since the ”mode-repulsion” property is coming from the localization of the modes we expect that the average mode length $`L_m=L/N_m`$ should be proportional to localization length $`\xi `$ too. In Fig.4, we plot $`N_m`$ vs the length of the systems $`L`$ when we increase the length from 80 cells to 320 cells and keep all other parameters the same. In Fig.4, we also plot the average mode length $`L_m`$ vs the localization length $`\xi `$ when we change the random strength W for a 320 cell system. The localization lengths are calculated using the transfer-matrix method by averaging 10,000 random configurations. These results confirm that indeed $`N_mL`$ and $`L_m\xi `$. It will be very interesting if these predictions can be checked experimentally.

The emission spectra at the right and left side of the system are quite different. This can be explained from the field patterns shown in Fig. 1a and Fig. 1c. Notice the localized modes are not similar at both sides of the system. This is the reason for this difference in the output spectrum. In Fig. 1d we denote with $`l`$ and $`r`$ the output modes from the left and the right side of our 1d system, respectively. The non-isotropic output spectra of real $`3D`$ experiments  might be explained by assuming that every localized mode has its intrinsic direction, strength and position, and the detected output spectra in experiments at different directions are the overlap of contributions from many modes. So generally they should be different. Most of the modes are not able to escape in our model and this is because of the 1D localization effects and exchange of energy between modes.

In summary, by using a FDTD method we constructed a random four-level lasing model to study the interplay of localization and amplification. Unlike the time-independent models, the present formulation calculates the field evolution beyond the threshold amplification. This model allows us to obtain the field pattern and spectra of localized lasing modes inside the system. For random systems, we can explain the multi-peaks and the non-isotropic properties in the emission spectra, seen experimentally. Our numerical results predict the ”mode-repulsion” property, the lasing-mode saturated number and average modelength. We also observed the exchange of energy between the localized modes which is much different from common lasers and this is essential for further research of mode competition and evolution in random laser. All of these properties are from the interplay of the localization and amplification where new physics phenomena can be found.

Ames Laboratory is operated for the U.S. Department of Energy by Iowa State University under Contract No. W-7405-Eng-82. This work was supported by the director for Energy Research, Office of Basic Energy Sciences.