text:
# Maps of the millimetre sky from the BOOMERanG experiment

## 1. Introduction

BOOMERanG is balloon-borne microwave telescope, sensitive at 90, 150, 240 and 410 GHz, with a resolution of $`10^{}`$. The instrument was equipped with very sensitive spider-web bolometers (Mauskopf et al. 1997). The image of the sky is obtained by slowly scanning the full payload in azimuth ($`\pm 30^o`$) at constant elevation. The scan center constantly tracks the azimuth of the lowest Galactic foreground region, situated in the southern hemisphere, in the Horologium constellation. Every day of the flight the instrument obtains a fully cross-linked map of about $`45^o\times 30^o`$ of the sky. The instrument was flown in a long duration circum-Antarctic flight from Dec.28, 1998 to Jan.8, 1999, and has been described in Piacentini et al. (2002) and in Crill et al. (2003). The main target of the experiment was the detection of anisotropy in the Cosmic Microwave Background. The maps produced by the experiment have been published in de Bernardis et al. (2000), Masi et al. (2001), Netterfield et al. (2002), Ruhl et al. (2003). Due to the limited size and to 1/f noise, these maps do not contain information for angular scales $`5^o`$ and have been filtered accordingly. The maps are calibrated to $`10\%`$ in gain and to $`10\%`$ in beam FWHM. In the 150 GHz map, the signal is well above the noise, and maps taken at different scan speeds and in different locations are perfectly consistent, demonstrating the low level of systematic effects. The level of the noise is of the order of 50 $`\mu K`$ per $`7^{}`$ pixel and of 20 $`\mu K`$ per $`28^{}`$ pixel. At 90, 150 and 240 GHz the rms signal has the spectrum of CMB anisotropy, and does not fit any reasonable spectrum of foreground emission. The temperature fluctuations detected in the high latitude part of the 150 GHz map are remarkably gaussian (Polenta et al. 2002, De Troia et al. 2003). The angular power spectrum of the 150 GHz map in the multipole range $`50<\mathrm{}<1000`$ has been computed in Netterfield et al. (2002) and Ruhl et al. (2003). Three peaks have been detected in the power spectrum, at multipoles $`\mathrm{}`$ 210, 540, 845 (de Bernardis et al., 2002). Interstellar Dust contamination of the 150 GHz power spectrum has been shown to be less than 1% (Masi et al. 2001). These results fit the scenario of acoustic oscillations of the primeval plasma at horizon and subhorizon scales (Sunyaev and Zeldovich 1970, Peebles and Yu 1970). In the framework of adiabatic inflationary structure formation the cosmological parameters have been estimated from the measured power spectrum (Lange et al. 2001, Netterfield et al., 2002, de Bernardis et al. 2002, Ruhl et al. 2003). The three cosmological parameters best constrained by the BOOMERanG data are the curvature parameter $`\mathrm{\Omega }=1.03\pm 0.05`$ (the universe is nearly flat), the spectral index of the primordial density perturbations $`n_s=1.02\pm 0.07`$ (nearly scale invariant), the physical density of baryons $`\mathrm{\Omega }_bh^2=0.023\pm 0.003`$ (consistent with Big Bang Nucleosynthesis).

Here we compare the BOOMERanG maps to the maps recently obtained at similar frequency and resolution by the WMAP satellite (Bennett et al., 2003). Working from the advantage L2 point of the Sun-Earth system, in its first year of operation WMAP has produced full sky maps of the microwave sky at 22, 32, 41, 60, 94 GHz, with resolution of the order of $`15^{}30^{}`$ and noise of $`35\mu K`$ per pixel in 28’ pixels. The maps are precisely calibrated (better than 1%). The power spectra obtained from these maps are fully consistent with the CMB power spectrum measured by BOOMERanG at 150 GHz, and with the adiabatic inflationary scenario: they allow the precise determination of most of the cosmological parameters; these estimates improve the precision of the BOOMERanG ones by a factor 2-3 for the parameters mentioned above; the full-sky coverage and polarization sensitivity allows the determination of previously poorly constrained parameters, like the reionization optical depth (Spergel et al. 2003).

The comparison of the BOOMERanG and WMAP maps is carried out with three targets:

$``$ Compare independent maps of the CMB with similar angular resolution, and confirm the detection of primordial structures in both experiments;

$``$ Improve the BOOMERanG calibration using the precise calibration of WMAP, and demonstrate that 1% gain calibration is possible for the new, polarization-sensitive BOOMERanG-B2K survey;

$``$ Infer the level and properties of foregorund contamination in forthcoming deep surveys of the CMB.

## 2. BOOMERanG vs WMAP

In fig. 1 we compare the maps of BOOMERanG at 90, 150, 240 GHz to the maps of WMAP at 41, 60, 94 GHz in the same high latitude region. The maps are Healpix (Gorski et al. 1998) representations with npix=512 (7’ per pixel). The region selected has the best coverage in the BOOMERanG 150 GHz channels. A $`5^o`$ high-pass has been applied to the WMAP maps to easy the comparison to the BOOMERanG maps, which are intrinsically insensitive to large angular scales.

The morphological agreement of the structures detected is evident. Significant differences are in the level of the noise and in the presence of important contamination by AGNs in the lower frequency channels. Once these sources are masked, the pixel-pixel correlation between the maps is quite good.

A quantitative analysis must take into account all important differences between the two datasets:

$``$ The beams of BOOMERanG and WMAP are different. The BOOMERanG beam at 150 GHz is closely fit by a 11’ gaussian down to 2% of the axial gain, while at lower levels is better fit by a 13’ gaussian. The WMAP beam at 94 GHz is fit by a 13’ gaussian down to 5% of the axial gain, and has wide ”shoulders” at lower levels. This difference is difficult to treat in pixel space, while its treatment is relatively simple in multipoles space. An example of this can be found in Abroe et al. (2003), comparing MAXIMA and WMAP maps of the CMB.

$``$ The BOOMERanG maps do not include structures larger than $`10^o`$, while WMAP maps are accurate at all scales. Once again, this is easier to account for in multipoles space.

$``$ In this pixelization, the noise of the 150 GHz BOOMERanG map is around 50 $`\mu K`$ per pixel, while the noise of the 94 GHz WMAP map is around 180 $`\mu K`$ per pixel.

A full analysis taking into account all these details is described in Hivon et al. (2003): the main result is that the calibration of BOOMERanG-98 at 150 GHz is found to be consistent with the precise calibration of WMAP within 5%. Moreover, $`<1\%`$ calibration of BOOMERanG is found to be reachable with this correlation method.

Given the close consistency of the maps of BOOMERanG and WMAP, the consistency of the angular power spectra measured by the two experiments is almost a trivial consequence. In fig.2 we compare the angular power spectra detected by the two instruments. The gain calibration error ($`<1\%`$ for WMAP and $`10\%`$ for BOOMERanG) is not included in the error bars, which account for random errors only. It is evident that the WMAP experiment has much higher sensitivity than BOOMERanG at multipoles $`\stackrel{<}{}450`$: the WMAP measurement in this range is limited by cosmic variance, and the Power Spectrum can be considered definitive. In the region of the second peak the two experiments have similar sensitivity, while in the region of the third peak the BOOMERanG experiment has better sensitivity due to the smaller beam and the lower noise per pixel. For the determination of the cosmological parameters $`\mathrm{\Omega }`$, $`\mathrm{\Omega }_bh^2`$, $`n_s`$, in the adiabatic inflationary scenario, WMAP takes advantage of the accurate calibration, while the BOOMERanG is still competitive because of the wider multipoles coverage, including the third peak of the spectrum.

These results give only an idea of what can be expected from Planck-HFI (Lamarre et al. 2003), the satellite instrument developed to fully exploit the capabilities of cryogenic bolometers. This will work for about two years from the same deep space location as WMAP, using bolometric detectors at 0.1K, even more sensitive than the BOOMERanG ones.

## 3. The polarization-sensitive BOOMERanG: B2K

After the 1998-99 flight, BOOMERanG has been recovered and upgraded. Additional attitude sensors have been implemented (a day-time star camera developed in Toronto and a pointed star sensor developed in Rome), and the focal plane has been rebuild to accomodate polarization sensitive bolometers (PSB) developed in JPL/Caltech (Jones et al. 2003). The new focal plane is sketched in fig.3.

It has 8 pixels in two rows. The distance between pixels projected in the sky is 0.5o both in azimuth and in elevation. Each pixel of the bottom row is a 145 GHz PSB with two independent bolometers sensing the two orthogonal polarization directions (Jones et al. 2003); each pixel on the top row senses a single polarization direction in two frequencies (240 and 350 GHz). The 150 GHz beam is $`9^{}`$ FWHM, while the 240 and 350 GHz beams are $`6^{}`$ FWHM. More details on the instrument can be found in Montroy et al. (2003). The experiment has been designed to measure the E modes (gradient) in the CMB polarization pattern. This has been detected only by the DASI experiment, but the sensitivity was not enough to constrain cosmological models more than with anisotropy data. WMAP has published a detection of the TE cross correlation from the first year of operation, consistent with the adiabatic model inferred by the anisotropy measurement. B2K should be able to measure the TE and EE power spectra at frequencies higher than the DASI and WMAP ones (30 to 94 GHz), thus facing different polarized foregrounds and nicely complementing them. The B2K payload has been flown from the McMurdo base on Jan.7, 2003, for a total of 11 days of operation in the stratosphere. In fig.4 we plot preliminary maps from B2K obtained from the PSBs at 145 GHz. We plan to re-fly B2K with an upgraded focal plane, to go after the polarized foreground from cirrus dust and AGNs. This information is essential for all the planned B-modes experiments (e.g. BICEP, Dome-C etc.) and is very difficult to measure from ground. The BOOMERanG optics can host an array of $`\stackrel{>}{}100`$ PSB at $`\stackrel{>}{}350`$ GHz, providing a deep, high resolution survey of polarized foreground emission at frequencies close to the ones used for CMB polarization research.

The BOOMERanG experiment is supported in Italy by Agenzia Spaziale Italiana, Programma Nazionale Ricerche in Antartide, Universita’ di Roma La Sapienza; by PPARC in the UK, by NASA, NSF OPP and NERSC in the U.S.A., and by CIAR and NSERC in Canada.

> `Abroe, E., et al., 2003, submitted to Ap.J,` astro-ph/0308355
>
> `Bennett, C., et al., 2003, Ap.J., in press,` astro-ph/0302207
>
> `Crill, B., et al. 2003, Ap.J.S., 148, 527`
>
> `de Bernardis, P., et al., 2000, Nature, 404, 955`
>
> `de Bernardis, P., et al., 2002, Ap.J., 564, 559`
>
> `De Troia, G., et al., 2003, MNRAS, 343, 284`
>
> `Gorski, K. M., Hivon, E. and Wandelt, B.D., 1998, in`
> `Analysis Issues for Large CMB Data Sets, A. J. Banday,`
> `R. K. Sheth and L. Da Costa eds., ESO, Ipskamp, NL,`
> `pgg. 37-42, see also` http://www.eso.org/science/healpix/
>
> `Hivon, E., et al., 2003, in preparation`
>
> `Jones, W.C., et al., 2002,` astro-ph/0209132
>
> `Lamarre, J.M., et al., 2003,` astro-ph/0308075
>
> `Lange, A.E., et al., 2001, Phys.Rev., D63, 042001`
>
> `Masi, S., et al., 2001, et al., Ap.J, 553, L93`
>
> `Mauskopf, P., et al. 1997, Applied Optics, 36, 765`
>
> `Montroy, T., et al., 2003,` astro-ph/0305593
>
> `Netterfield, B., et al., 2002, Ap.J., 571, 604`
>
> `Piacentini, F., et al. 2002, Ap.J.Suppl., 138, 315`
>
> `Polenta, G., et al., 2002, Ap.J., 572, L27`
>
> `Peebles, P.J.E, and Yu, J.T., 1970, Ap.J., 162, 815`
>
> `Ruhl, J., et al., 2003, Ap.J. in press` ( astro-ph/0212229 )
>
> `Spergel, D.N., et al., 2003, Ap.J., in press,` astro-ph/0302209
>
> ```
> Sunyaev, R.A. & Zeldovich, Ya.B., 1970, Astr.SpaceSci., 7, 3
>
> ```


format: text

text:
# Non-linear and quantum optics of a type II OPO containing a birefringent element Part 2 : bright entangled beams generation

## 1 Introduction

The domain of quantum information with intense light beams is becoming more and more studied as practical applications are developed teleportkimble ; grangier ; leuchs . A large number of experiments rely on so-called EPR beams, that is a pair of fields containing quantum correlations on two different non commuting observables epr . Several methods have been used to date. One way to generate them is to mix two independent squeezed beams produced by two independent OPOs below threshold teleportkimble ; eprbachor , or to use third order non linearity in a fiber leuchs . Another configuration consists in using a non degenerate OPO below threshold eprkimble ; eprtaiyuan ; eprpolzik . We propose here a new method where two bright entangled beams are produced directly at the output of a non degenerate OPO pumped above threshold.

Optical Parametric Oscillators are very efficient sources of non classical light taiyuan ; kimblevidesq ; mertz ; mlynek ; zhang ; mescond . In particular, it has been shown experimentally that the intensities of the signal and idler beams exhibit very large quantum correlations taiyuan ; mertz ; mescond . It is also predicted that the phase of the signal and idler should show quantum anticorrelations reynaud . However, the measurement of these phase anticorrelations is not possible in a standard OPO in which frequency degenerate operation occurs only accidentally and is very sensitive to thermal and mechanical drifts. Thus, in practical experiments, the frequencies of the signal and idler beam are different. This implies that the relative phase between the two fields has a rapidly varying component. Furthermore, even if the device is actively stabilized on the frequency degeneracy working point pfister , this relative phase between signal and idler undergoes a diffusion process courtois , similarly to the phase of a laser field. In order to be able to measure experimentally these phase anticorrelations, one can use a phase locked OPO : one possible technique proposed in wong  is to use a birefringent plate inserted inside the OPO cavity and making an angle with the axis of the non-linear crystal which leads to self-phase locking. This phenomenon has been studied in detail in fabre ; part1  at the classical level. It is the purpose of this paper to evaluate the influence of the phase locking mechanism on the quantum properties of the system.

The paper is organized as follows. In section 2, we give the equations describing the quantum fluctuations of the system. We study in section 3 the noise spectra of the amplitude quadratures difference and phase quadratures sum showing quantum noise reduction, *i.e* that these variances are below what would be observed for independent coherent states. In section 4, various entanglement criteria are discussed, namely EPR correlation properties (section 4.1), inseparability criterion (section 4.2) and teleportation fidelity (section 4.3). Finally, section 5 is devoted to the case of non-orthogonal quadratures which enable us to increase dramatically the violation of the inseparability criterion.

## 2 Basic equations

In order to calculate the fluctuations for the signal and idler fields when the OPO is pumped above threshold, we apply the input-output linearization technique described *e.g.* in fabre90 . We first determine the classical stationary solutions of the evolution equations : these solutions are denoted $`A_i^{(s)}`$ where $`i=0,1,2`$ corresponds to the pump, signal and idler fields. We then linearize the evolution equations around these stationary values by setting : $`A_i=A_i^{(s)}+\delta A_i,i=0,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}2}`$. Introducing the input fluctuations in all the modes, we obtain a set of linear coupled equations for the fluctuations of the different fields. By solving these equations in the Fourier domain, one can calculate the individual noise spectra as well as the correlations between the different quadratures.

Let us recall the system under study which has been discussed in detail in part1 . We consider a $`\chi ^{(2)}`$ medium placed inside a ring optical cavity doubly resonant for signal and idler (fig. 1).

Type-II phase-matching results in orthogonally polarized signal and idler fields. The reflection coefficients for signal and idler are assumed to be equal in module : $`r_j=r\mathrm{exp}(i\zeta _j),j=1,2`$. We assume also that the cavity has a high finesse for these modes so that we can write $`r=1\kappa `$ with $`\kappa 1`$ ; the transmission is thus $`t\sqrt{2\kappa }`$. Finally, the cavity is characterized by the round-trip losses for the signal and idler waves (crystal absorption, surface scattering), $`\mu `$, and we define a generalized reflection coefficient : $`r^{}=r(1\mu )1(\mu +\kappa )=1\mu ^{}`$.

Let us denote $`l`$ the crystal length and $`L`$ the free propagation length inside the cavity. The crystal indices of refraction are $`n_1`$ and $`n_2`$ respectively for the signal and idler waves. The waveplate is rotated by an angle $`\rho `$ with respect to the crystal neutral axes. In contrast to part1 , we assume here that this angle is small.

The crystal is birefringent and can be characterized by its mean index $`\overline{n}`$, birefringence $`\theta `$ and effective nonlinear coupling constant $`g`$ :

$`\overline{n}`$ $`=`$ $`{\displaystyle \frac{n_1+n_2}{2}}`$ (1)
$`\theta `$ $`=`$ $`{\displaystyle \frac{\omega _0}{2c}}(n_2n_1)l`$ (2)
$`g`$ $`=`$ $`l\chi ^{(2)}\sqrt{{\displaystyle \frac{\mathrm{}\omega _0\omega _1\omega _2}{2c^3\epsilon _0n_0n_1n_2}}}`$ (3)
where $`\omega _0/2`$ is the common frequency of signal and idler.

The mean round trip phase-shift is $`\delta =\frac{\omega _0}{2c}(ne+\overline{n}l+L)+\overline{\zeta }`$ where $`\overline{\zeta }`$ is the mean mirror birefringence, $`\overline{\zeta }=\frac{\zeta _1+\zeta _2}{2}`$. We assume that the signal and idler modes are close to resonance so that the corresponding round trip phase shifts

$$\delta _{1,2}=\delta \pm \frac{\theta +\zeta _2\zeta _1}{2}\psi $$
(4)
where $`\psi `$ is the birefringence of the waveplate, are close to an integer multiple of $`2\pi `$ :

$$\delta _{1,2}=2\pi p_{1,2}+\mathrm{\Delta }_{1,2}$$
(5)
where $`p_{1,2}`$ are integers and $`\mathrm{\Delta }_{1,2}`$ are small detunings.

In this case, the classical evolution equations given in part1  become :

$$\begin{array}{ccc}\tau \frac{dA_1^{(s)}}{dt}& =& A_1^{(s)}(\mu ^{}+i\mathrm{\Delta }_1)+gA_0^{(s)}A_2^{(s)}+i2\rho e^{i(\theta \psi )}A_2^{(s)}\hfill \\ \tau \frac{dA_2^{(s)}}{dt}& =& A_2^{(s)}(\mu ^{}+i\mathrm{\Delta }_2)+gA_0^{(s)}A_1^{(s)}+i2\rho e^{i(\theta \psi )}A_1^{(s)}\hfill \\ A_0^{(s)}& =& A_0^{in}\frac{g}{2}A_1^{(s)}A_2^{(s)}\hfill \end{array}$$
(6)
where $`A_0^{in}`$ is the input pump amplitude at frequency $`\omega _0`$.

The stationary solution of these equations, corresponding to the phase locked regime, can be found only within a given range of detunings $`(\mathrm{\Delta }_1,\mathrm{\Delta }_2)`$ which has been detailed in part1 .

The evolution equations for the fluctuations then are :

$$\begin{array}{ccc}\tau \frac{d\delta A_1}{dt}& =& \delta A_1(\mu ^{}+i\mathrm{\Delta }_1)+\hfill \\ & & g\left(A_0^{(s)}\delta A_2^{}+A_2^{(s)}\delta A_0\right)+i2\rho e^{i(\theta \psi )}\delta A_2+\hfill \\ & & \sqrt{2\kappa }\delta A_1^{in}+\sqrt{2\mu }\delta B_1^{in}\hfill \\ \tau \frac{d\delta A_2}{dt}& =& \delta A_2(\mu ^{}+i\mathrm{\Delta }_2)+\hfill \\ & & g\left(A_0^{(s)}\delta A_1^{}+A_1^{(s)}\delta A_0\right)+i2\rho e^{i(\theta \psi )}\delta A_1+\hfill \\ & & \sqrt{2\kappa }\delta A_2^{in}+\sqrt{2\mu }\delta B_2^{in}\hfill \\ \delta A_0^{(s)}& =& \delta A_0^{in}\frac{g}{2}\left(A_2^{(s)}\delta A_1+A_1^{(s)}\delta A_2\right)\hfill \end{array}$$
(7)
where $`\delta A_i^{in}`$ and $`\delta B_i^{in}`$ correspond to the vacuum fluctuations entering the cavity due respectively to the coupling mirror and to the losses. From these equations, one can calculate the Fourier components of the fluctuations and thus the spectra of variances and correlations for different quadratures of the signal and idler fields.

## 3 Variances

The amplitude quadrature and phase quadrature fluctuations can be written in the Fourier domain as

$`\delta \stackrel{~}{x}_i(\mathrm{\Omega })`$ $`=`$ $`\delta \stackrel{~}{A}_i(\mathrm{\Omega })e^{i\phi _i}+\delta \stackrel{~}{A}_i^{}(\mathrm{\Omega })e^{i\phi _i},`$ (8)
$`\delta \stackrel{~}{p}_i(\mathrm{\Omega })`$ $`=`$ $`i\left(\delta \stackrel{~}{A}_i(\mathrm{\Omega })e^{i\phi _i}\delta \stackrel{~}{A}_i^{}(\mathrm{\Omega })e^{i\phi _i}\right),`$ (9)
$`i=1,\mathrm{\hspace{0.17em}2}`$ where $`\delta \stackrel{~}{A}_i`$ denotes the Fourier component and $`\phi _i`$ is the phase of the stationary solution $`A_i^{(s)}`$.

To characterize the entanglement, one can calculate the normalized variances of the sum and difference of the amplitude and phase quadratures eprkimble ; reid

$$\begin{array}{ccc}S_x^\pm (\mathrm{\Omega })& =& \frac{1}{2}\left|\delta \stackrel{~}{x}_1(\mathrm{\Omega })\pm \delta \stackrel{~}{x}_2(\mathrm{\Omega })\right|^2\\ S_p^\pm (\mathrm{\Omega })& =& \frac{1}{2}\left|\delta \stackrel{~}{p}_1(\mathrm{\Omega })\pm \delta \stackrel{~}{p}_2(\mathrm{\Omega })\right|^2\end{array}$$
(10)
The fact that $`S_x^{}`$ and $`S_p^+`$ are smaller than one is the signature of simultaneous quantum correlations on both phase and amplitude quadratures. This is the case in a standard OPO above threshold with the restriction already mentioned in the introduction that because of the phase diffusion, $`S_p^+`$ can not be measured easily.

¿From the above calculated fluctuations, one can obtain the variances for any operating point. For the sake of simplicity, we consider here the case of detunings $`\mathrm{\Delta }=\pm 2\rho `$ which correspond to the minimum threshold for which $`\theta \psi =\frac{\pi }{2}`$ and $`A_i^{(s)}=\sqrt{\frac{2\mu ^{}}{g}(\sigma 1)}e^{i\phi _i},i=1,2`$ where $`\sigma `$ is the input pump intensity normalized to the threshold of the standard OPO. We now have :

$`S_x^{}(\mathrm{\Omega })`$ $`=`$ $`{\displaystyle \frac{\left(\mathrm{\Omega }^2\frac{4\rho ^2}{\mu ^2}\right)^2+\frac{4\rho ^2}{\mu ^2}\frac{\kappa }{\mu ^{}}+\frac{\mu }{\mu ^{}}\mathrm{\Omega }^2}{\left(\mathrm{\Omega }^2\frac{4\rho ^2}{\mu ^2}\right)^2+\mathrm{\Omega }^2}}`$ (11)
$`S_p^+(\mathrm{\Omega })`$ $`=`$ $`1{\displaystyle \frac{\frac{\kappa }{\mu ^{}}}{\mathrm{\Omega }^2+\sigma ^2}}`$ (12)
where $`\mathrm{\Omega }=\frac{\omega \tau }{2\mu ^{}}=\frac{\omega }{\mathrm{\Omega }_c}`$ is the frequency normalized to the cavity cut-off frequency, $`\mathrm{\Omega }_c`$.

The expression for the phase sum variance is identical to the one obtained in the case of a standard (non phase locked) OPO : one obtains very large phase anticorrelations for any waveplate angle. The expression for the intensity difference variance is different from the expression obtained in the case of a standard OPO and does not depend on the pump power $`\sigma `$. The intensity difference noise reduction is now degraded at low frequencies by the phase locking : inhibition of phase diffusion prevents $`S_p^{}`$ from going to $`+\mathrm{}`$ as $`\mathrm{\Omega }`$ goes to zero. Thus, in order to verify the corresponding Heisenberg inequality, $`S_p^+`$ cannot go to zero but instead takes a finite value.

The corresponding spectra are shown on fig 2 and 3.

These figures show that the amplitude quadrature difference variance goes to zero as $`\rho `$ goes to zero while the phase quadrature sum goes to zero close to threshold.

An important characteristic of the system is the minimal value obtained for $`S_x^{}`$ as a function of $`\mathrm{\Omega }`$, $`\left(S_x^{}\right)_{min}`$ which is plotted in fig. 4 as a function of $`\rho `$.

This plot shows that one obtains a significant noise reduction as long as $`\rho `$ is small compared to $`\kappa `$.

## 4 Entanglement criteria

### 4.1 EPR correlations criterion

In the Gedanken Experiment discussed in epr , one tries to infer the quantum state of a particle by measurements on a spatially separated particle : the optical equivalent of that scheme is to try to infer the signal beam by measurements on the spatially separated idler beam. In this case, quantum correlations do not suffice to characterize the efficiency of the system and another criterion has to be considered reid . The noise on the measured observable has to be taken into account : a large noise will reduce the precision with which the unmeasured system is known. One is lead to introduce inference errors or conditional variances reid . The conditional variances for amplitude and phase are given by :

$`V_x`$ $`=`$ $`\underset{g}{\mathrm{min}}|\delta \stackrel{~}{x}_1g\delta \stackrel{~}{x}_2|^2=|\delta \stackrel{~}{x}_1C_x\delta \stackrel{~}{x}_2|^2`$ (13)
$`V_p`$ $`=`$ $`\underset{g}{\mathrm{min}}|\delta \stackrel{~}{p}_1g\delta \stackrel{~}{p}_2|^2=|\delta \stackrel{~}{p}_1C_p\delta \stackrel{~}{p}_2|^2`$ (14)
where

$`C_x`$ $`=`$ $`{\displaystyle \frac{\delta \stackrel{~}{x}_1(\mathrm{\Omega })\delta \stackrel{~}{x}_2(\mathrm{\Omega })}{\sqrt{\delta \stackrel{~}{x}_1^2\delta \stackrel{~}{x}_2^2}}}`$ (15)
$`C_p`$ $`=`$ $`{\displaystyle \frac{\delta \stackrel{~}{p}_1(\mathrm{\Omega })\delta \stackrel{~}{p}_2(\mathrm{\Omega })}{\sqrt{\delta \stackrel{~}{p}_1^2\delta \stackrel{~}{p}_2^2}}}`$ (16)
are the normalized correlations. Ref. reid  shows that one must have $`V_xV_p<1`$ in order to have a demonstration of the EPR paradox.

Even with correlations factors close but not equal to unity, the conditional variances can be large if the quantum noise on the individual beams is large as is the case for the amplitude quadrature.

The conditional variances $`V_x`$ and $`V_p`$ are shown on fig. 5 and 6. Due to the large excess noise on the amplitude quadrature at low frequencies, $`V_x`$ is large for low frequencies while $`V_p`$ remains low for all frequencies. This difference between the two conditional variances is one of the main differences between our case and the other devices which have been introduced to generate bright entangled beams eprbachor ; eprkimble ; eprtaiyuan ; eprpolzik : in these cases, one obtains theoretically symmetric behaviors for amplitude and phase which is not the case in our system.

The variation of $`V_xV_p`$ with $`\mathrm{\Omega }`$ is plotted on fig. 7. This figure shows that EPR correlations exist in a broad frequency range even for large values of the waveplate angle. The behavior of the criterion is dominated by that of $`V_p`$ which goes to zero as $`\mathrm{\Omega }`$ goes to zero and is independent of the waveplate angle.

### 4.2 Inseparability criterion

However, this criterion is not always relevant. In the vast majority of quantum information protocols, what matters is the non-separability of the quantum state that is the impossibility to express it as a product of independent sub-states. In order to evaluate the inseparability of two systems, Peres proposed a necessary criterion which applies to finite dimension systems peres  which was latter shown by Horodecki horodecki  to be sufficient for dimensions smaller than $`2\times 2`$ \- or $`2\times 3`$ -. In the particular case of Gaussian states, Duan *et al.* duan  and Simon simon  developed a necessary criterion for inseparability which can be written as :

$$S_x^{}+S_p^+2$$
(17)
where $`S_x^{}`$ and $`S_p^+`$ are the previously introduced variances spectra. Let us note the entanglement produced by the self-phase locked OPO is not in the so-called standard form (duan , equations 10-11) : in this case, the above criterion is a sufficient but not necessary condition for inseparability.

The value of $`S_x^{}+S_p^+`$ as a function of $`\mathrm{\Omega }`$ is plotted in fig 8. At low frequencies, the behavior is dominated by the behavior of $`S_x^{}`$ which is large. As higher frequencies are considered, the figure shows that the system is indeed entangled over a large range of frequencies. When $`\rho \kappa `$, the excess noise on the amplitude quadrature makes the sum larger than two meaning that the system is separable. The curves are plotted for a value $`\sigma =1`$, that is at exact threshold. Increasing $`\sigma `$ does not increase significantly the variance $`S_x^{}`$ and does not modify $`S_p^+`$ and will thus move all the curves toward the inseparability limit while keeping the sum below two. It is the behavior of the amplitude fluctuations which dominates that of the sum of the variances.

### 4.3 Fidelity of teleportation

In teleportation experiments, one tries to transfer quantum states from one point to another with a unity gain. In this case, the successful of the teleportation is characterized by the overlap between the input and output states. This overlap is usually called the fidelity teleportkimble . It can be shown that in order to obtain a fidelity larger than $`\frac{1}{2}`$, quantum correlations must be used, while a fidelity larger than $`\frac{2}{3}`$ denotes the possession by the receiver of the best copy available grangier . In the particular case of gaussian states111and only in this case, one can show that the fidelity is related to the intensity difference and phase sum variances defined in equation 10 by the relation :

$$=\frac{1}{\sqrt{(1+S_x^{})(1+S_p^+)}}$$
(18)
The fidelity as a function of the normalized frequency is shown in fig 9 for different angles of the waveplate.

One sees that for small values of the waveplate angle $`\rho `$, one obtains a fidelity larger that the ”best copy” limit ($`2/3`$) over a large frequency range.

In order to determine the optimal operating parameters, one can plot the values of waveplate angle $`\rho `$, normalized intensity $`\sigma `$ and analysis frequency $`\mathrm{\Omega }`$ for which the fidelity is equal to the ”best-copy” limit, that is $`=\frac{2}{3}`$. The plot consists of two zones : close to threshold, *i.e.* $`\sigma <1.1`$, one obtains a fidelity larger than $`\frac{2}{3}`$ for large values of $`\rho `$ and for $`\mathrm{\Omega }<0.5`$ (dark grey surface) ; for $`\sigma >1.1`$, a large fidelity is observed over a larger range of frequencies but only for very small values of $`\rho `$, typically $`\rho <0.02`$ (light grey surface).

One observes that in order to obtain a fidelity significantly larger than $`\frac{2}{3}`$, one can use very small angles, typically less than $`1^{}`$ (0.015 radians). As mentioned in ref. part1 , this is feasible experimentally, since a small angle does not reduce the range of the crystal temperature and cavity length where the locking is observed but only its total surface : once frequency degeneracy has been obtained, existing active stabilization of the system allows for operation in this zone. However, working very close to threshold will reduce the parameter range where phase locking is observed and may be difficult experimentally.

We have shown in this section that our system produces bright entangled beams which fulfill different entanglement criteria over a large range of frequencies. At low frequencies, the entanglement is observed for $`\rho \kappa `$, that is for a rotation of the waveplate small compared to the mirror transmission.

## 5 Use of optimized quadratures

The limitation of the EPR correlations zone to small angle is due to the fact that the intensity correlations are no longer perfect in a phase locked OPO as was the case in the standard OPO. A simple picture shows the origin of this phenomenon. $`A_1`$ and $`A_2`$ having a fixed phase relation, the output beam is in a fixed polarization state which we denote $`A_{bright}`$. The orthogonal polarization mode has a zero mean value and is denoted $`A_{dark}`$. One can show that they correspond to the polarization eigenmodes of the ”cold” OPO, that is without the pump beam eigenmodes .

One can also show that their variance spectra are

$$\begin{array}{ccc}S_x(A_{dark})& =& |\delta \stackrel{~}{x}_{dark}|^2=\frac{1}{2}|\delta \stackrel{~}{x}_1\delta \stackrel{~}{x}_2|^2=S_x^{}\\ S_p(A_{bright})& =& |\delta \stackrel{~}{p}_{bright}|^2=\frac{1}{2}|\delta \stackrel{~}{p}_1+\delta \stackrel{~}{p}_2|^2=S_p^+\end{array}$$
(19)
The exact expressions of the mean fields $`A_{bright}`$ and $`A_{dark}`$ depend on the choice of the operating point. For the minimum threshold point where we have performed our calculations, we have

$`A_{bright}`$ $`=`$ $`{\displaystyle \frac{A_1+A_2}{\sqrt{2}}}`$ (20)
$`A_{dark}`$ $`=`$ $`{\displaystyle \frac{A_1A_2}{\sqrt{2}}}`$ (21)
which are left and right circularly polarized due to the $`\frac{\pi }{2}`$ phase shift between $`A_1`$ and $`A_2`$.

As we have seen previously that $`S_x^{}`$ and $`S_p^+`$ can be smaller than the corresponding standard quantum limit, $`A_{bright}`$ is phase squeezed while $`A_{dark}`$ is squeezed on the orthogonal quadrature. This is shown in fig. 11 which represents the Fresnel diagram common to the fields $`A_1`$, $`A_2`$, $`A_{bright}`$ and $`A_{dark}`$, the phase reference being the phase of the field $`A_{bright}`$.

The increase in $`S_x^{}`$ at low frequencies corresponds to to the tilt of the corresponding squeezing ellipse (in black) with respect to the standard OPO situation (vertical ellipse in dark grey). Due to the presence of the waveplate, the amplitude quadrature of $`A_{dark}`$ (measurement of the width along the horizontal axis) is not the optimal quadrature : it is thus natural to detect this optimal quadrature instead of the amplitude quadraturejosse . It is possible to dephase222for instance, by using a combination of $`\lambda /2`$ and $`\lambda /4`$ waveplates. appropriately $`A_1`$ with respect to $`A_2`$ so that $`A_{bright}`$ and $`A_{dark}`$ are now squeezed on orthogonal quadratures. In the simple case of the minimum threshold point, the new entangled are simply obtained by a linear superposition of the original fields along the neutral axis of the crystal. In the simple case of the minimum threshold point, these fields are given by :

$`A_1^{}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left[\mathrm{cos}\left({\displaystyle \frac{\theta _{opt}}{2}}\right)A_1i\mathrm{sin}\left({\displaystyle \frac{\theta _{opt}}{2}}\right)A_2\right]`$ (22)
$`A_2^{}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left[\mathrm{sin}\left({\displaystyle \frac{\theta _{opt}}{2}}\right)A_1+i\mathrm{cos}\left({\displaystyle \frac{\theta _{opt}}{2}}\right)A_2\right]`$ (23)
where $`\theta _{opt}`$ is given by :

$$\theta _{opt}=\frac{1}{2}\mathrm{arctan}\left(\frac{2C_{}\sqrt{S_x\left(A_{dark}\right)S_p\left(A_{dark}\right)}}{S_x\left(A_{dark}\right)S_p\left(A_{dark}\right)}\right)$$
(24)
and

$$C_{}=\frac{\delta \stackrel{~}{x}\left(A_{dark}\right)\delta \stackrel{~}{p}\left(A_{dark}\right)}{\sqrt{S_x\left(A_{dark}\right)S_p\left(A_{dark}\right)}}$$
(25)
The inseparability criterion developed by Duan *et al.* duan  is

$$S_x^+\left(A_{bright}^{}+\right)+S_p\left(A_{dark}^{}\right)2$$
(26)
where $`S_{x,p}\left(A_{bright,dark}\right)`$ are defined above and $`A_{bright,dark}^{}`$ are the optimally squeezed states :

$$A_{bright,dark}^{}=\frac{A_1^{}\pm A_2^{}}{\sqrt{2}}$$
(27)
Setting $`\rho `$ and $`\sigma `$, one can plot $`S_p\left(A_{bright}^{}\right)+S_x\left(A_{dark}^{}\right)`$ as a function of the frequency (figure 12). This figure shows that one can obtain a much larger violation of the inseparability criterion using optimized quadratures. For $`\rho =0.01\text{rad}`$, the value obtained is very close to the curve obtained for a standard OPO (where frequency non degeneracy and phase diffusion prevent the measurement). Furthermore, the divergence obtained at low frequencies is completely eliminated. For larger values of the waveplate angle, the improvement is still noticeable.

## 6 Conclusion

We have shown that a type II OPO containing a birefringent element can be a very efficient source of intense entangled beams. This entanglement has been characterized by the use of several criteria, namely intensity difference and phase sum variances, inference errors which denote the possibility to make an EPR-type *Gedanken Experiment*, inseparability criteria which correspond to the quality of the entanglement and fidelity which represent the quality of teleportation with unity gain. These criteria have been shown to be fulfilled for a large range of operating parameters thus showing the quality of the system. The inseparability criterion can be largely improved by the use of optimized quadratures. We are currently setting up an experiment based on this system.

###### Acknowledgements.

Laboratoire Kastler-Brossel, of the École Normale Supérieure and the Université Pierre et Marie Curie, is associated with the Centre National de la Recherche Scientifique.
This work was supported by European Community Project QUICOV IST-1999-13071.
T. Coudreau is also at the Pôle Matériaux et Phénomènes Quantiques FR CNRS 2437, Université Denis Diderot, 2, Place Jussieu, 75251 Paris cedex 05, France.
We thank V. Josse and N. Treps for fruitful discussions.


text:
# Distribution of reflection eigenvalues in many-channel chaotic cavities with absorption

(November 12, 2003 )

## Abstract

The reflection matrix $`R=S^{}S`$, with $`S`$ being the scattering matrix, differs from the unit one, when absorption is finite. Using the random matrix approach, we calculate analytically the distribution function of its eigenvalues in the limit of a large number of propagating modes in the leads attached to a chaotic cavity. The obtained result is independent on the presence of time-reversal symmetry in the system, being valid at finite absorption and arbitrary openness of the system. The particular cases of perfectly and weakly open cavities are considered in detail. An application of our results to the problem of thermal emission from random media is briefly discussed.

Introduction.– When absorption in a cavity is finite, part of the incoming flux gets irreversibly lost in the walls, breaking thus the unitarity of the scattering matrix. The mismatch between incoming and outgoing fluxes in the scattering process can be naturally described by means of the reflection matrix $`R=S_\gamma ^{}S_\gamma `$, where $`S_\gamma `$ denotes the subunitary (at nonzero absorption) scattering matrix to be precisely defined below. The matrix $`R`$ has the positive eigenvalues $`r_c1`$, the so-called reflection eigenvalues. Their statistical properties in chaotic resonance scattering attract much attention presently, both experimentally Doron et al. (1990); Méndez-Sánchez et al. (2003) and theoretically Bruce and Chalker (1996); Beenakker et al. (1996); Beenakker (1998); Kogan et al. (2000); Ramakrishna and Kumar (2000); Savin and Sommers (2003); Fyodorov (2003), using the random-matrix theory approach Verbaarschot et al. (1985); Beenakker (1997).

In the limit of weak absorption, $`R`$ was found Doron et al. (1990); Ramakrishna and Kumar (2000) to be related to the Smith’s time-delay matrix at zero absorption those statistical properties were extensively studied in recent time Lehmann et al. (1995a); Fyodorov and Sommers (1997); Savin et al. (2001); Brouwer et al. (1997); Sommers et al. (2001). Such an analysis was recently generalized by us to the case of arbitrary absorption Savin and Sommers (2003), see also Fyodorov (2003) for the related study. In the particular case of the single-channel cavity, obtained results explain partly the recent experiment Méndez-Sánchez et al. (2003) on the reflection coefficient distribution.

The opposite case of a large number of equivalent channels plays an important role in chaotic scattering, corresponding to the semiclassical limit of matrix models Lewenkopf and Weidenmüller (1991); Lehmann et al. (1995b). In this case $`N`$ strongly overlapping resonances (poles of the $`S`$-matrix) are excited in the open cavity through the scattering channels, whose number $`M`$ scales with $`N`$. A new energy scale appears then in addition to the mean level spacing $`\mathrm{\Delta }`$ of the closed cavity Lehmann et al. (1995b, a): the empty gap between the real axis in the complex energy plane and the cloud of the resonances in its lower half plane. In the physically justified limit $`M/N1`$ Lehmann et al. (1995a) (though both $`M,N\mathrm{}`$) considered below this gap is given by the well-known Weisskopf width $`\mathrm{\Gamma }_W=MT\mathrm{\Delta }/2\pi `$, with $`T1`$ being the transmission coefficient Lewenkopf and Weidenmüller (1991). In the presence of finite absorption the widths of the resonances acquire an additional contribution, the absorption width $`\mathrm{\Gamma }_a`$ Savin and Sommers (2003). The two competing processes, the escape into continuum and dissolution into the walls, are therefore characterized by the ratio $`\mathrm{\Gamma }_a/\mathrm{\Gamma }_W=\gamma /T`$, where the convenient dimensionless parameter $`\gamma 2\pi \mathrm{\Gamma }_a/\mathrm{\Delta }M`$ measures the absorption strength.

In this paper we derive the probability distribution function $`𝒫(r)=M^1_{c=1}^M\overline{\delta (rr_c)}`$ of reflection eigenvalues, with the bar indicating the statistical average, in the limit of the large number of channels. This function has, in particular, an important application for the statistics of thermal emission from random media and is known only at perfect coupling ($`T=1`$) Beenakker (1998). Here, we calculate it at arbitrary $`T`$.

Our starting point is the general relation established in Ref. Savin and Sommers (2003) between $`R`$ and the effective non-Hermitian Hamiltonian $`=H(i/2)VV^{}`$ of the open system, with the Hermitian part $`H`$ standing for the closed counterpart and the amplitudes $`V_n^c`$ describing the coupling between $`N`$ intrinsic and $`M`$ channel states. It is achieved by making use of the equivalence Kogan et al. (2000); Ramakrishna and Kumar (2000); Savin and Sommers (2003) of uniform absorption to the pure imaginary shift $`E_\gamma E+\frac{i}{2}\mathrm{\Gamma }_a`$ of the scattering energy $`E`$, thus giving $`S_\gamma S(E_\gamma )=1iV^{}(E_\gamma )^1V`$. This leads to the following connection between the reflection matrix $`R`$ and the time-delay matrix with absorption $`Q_\gamma `$ ($`\mathrm{}=1`$) Savin and Sommers (2003):

$$R=1\mathrm{\Gamma }_aQ_\gamma ,$$
(1)
where $`Q_\gamma Q(E_\gamma )=V^{}[(E_\gamma )^{}]^1(E_\gamma )^1V`$. The distribution function $`𝒫(r)`$ is therefore directly related as

$$𝒫(r)=\gamma ^1𝒫_\tau \left(\gamma ^1(1r)\right)$$
(2)
to the distribution $`𝒫_\tau (\tau )=M^1_{c=1}^M\overline{\delta (\tau Mq_c/t_H)}`$ of the proper delay times (eigenvalues $`q_c`$ of $`Q_\gamma `$) scaled with $`M`$ and measured in units of the Heisenberg time $`t_H=2\pi /\mathrm{\Delta }`$.

The function (2) is normalized to unity. The average reflection coefficient $`r=M^1\overline{\mathrm{tr}R}`$ is given by

$$r=_0^1𝑑rr𝒫(r)=\frac{T+\gamma (1T)}{T+\gamma }.$$
(3)
This result follows from the connection Savin and Sommers (2003) between $`r`$ and the “norm-leakage” decay function $`P(t)`$ introduced in Savin and Sokolov (1997): $`r=1\gamma [1\mathrm{\Gamma }_a_0^{\mathrm{}}𝑑te^{\mathrm{\Gamma }_at}P(t)]`$. In the limit considered $`P(t)`$ reduces to the simple exponent $`e^{\mathrm{\Gamma }_Wt}`$, since the widths ceased to fluctuate Lehmann et al. (1995b); Savin and Sokolov (1997). Below we will also derive Eq. (3) using a different method.

The saddle-point equation.– The scaling limit studied is far from being trivial and can not be obtained from the known expressions for $`𝒫(r)`$ at finite $`M`$ Savin and Sommers (2003) simply by letting $`M\mathrm{}`$ there. One needs to start from the original definition that determines the seeking distribution through the jump of the resolvent $`G(t)=M^1\overline{\mathrm{tr}(t\mathrm{\Gamma }_aQ_\gamma )^1}1/t+(\gamma /t^2)K(t)`$ along the discontinuity line $`t>0`$ as follows:

$$𝒫(r)=\frac{\gamma }{\pi (1r)^2}\mathrm{Im}K(ti0)|_{t=1r}.$$
(4)
Following Refs. Savin and Sommers (2003); Sommers et al. (2001), one can obtain the exact (at $`N\mathrm{}`$) representation for $`K`$, which is convenient for our purposes,

$$K(t)=\frac{1}{8}\mathrm{str}\left[\left(\mathrm{\Lambda }i\frac{1t/2}{\sqrt{1t}}\mathrm{\Lambda }_1\right)k\widehat{\sigma }\right]_{},$$
(5)
with the shorthand $`(\mathrm{})_{}=𝑑\mu (\widehat{\sigma })(\mathrm{})e^{(M/2)}`$ and

$``$ $`=`$ $`(\gamma /t)\mathrm{str}[(1t/2)\mathrm{\Lambda }\widehat{\sigma }i\sqrt{1t}\mathrm{\Lambda }_1\widehat{\sigma }]`$ (6)
$`\mathrm{str}\mathrm{ln}[1+(2g)^1(\mathrm{\Lambda }\widehat{\sigma }+\widehat{\sigma }\mathrm{\Lambda })],`$
as the integral over the noncompact saddle-point manifold of the $`8\times 8`$ supermatrix $`\widehat{\sigma }`$ subject to the constraint $`\widehat{\sigma }^2=1`$. $`\mathrm{\Lambda },\mathrm{\Lambda }_1`$ appearing above are the supermatrix analog of the Pauli matrices $`\sigma _3,\sigma _1`$ and $`k=+1(1)`$ in the space of commuting (anticommuting) variables. Definitions of the superalgebra as well as explicit parametrization, which depends on whether TRS is preserved or broken, can be found in Verbaarschot et al. (1985); Efetov (1996). At last, the constant $`g=2/T11`$ is related to the transmission coefficient $`T`$.

In the limit $`M\mathrm{}`$ the integration over $`\widehat{\sigma }`$ can also be done in the saddle-point approximation. Function (5) is given by a saddle-point value of $`\widehat{\sigma }`$, which is found, as usual, by equating the variance $`\delta `$ to zero. The structure of the latter equation is determined by that of $``$. A careful analysis shows that, independently on TRS, the supermatrix $`\widehat{\sigma }`$ reduces (and so does the corresponding algebra) in the saddle-point to the $`2\times 2`$ usual matrix $`\sigma =a\sigma _3+b\sigma _1`$ with the imposed constraint $`a^2+b^2=1`$. We find that $`K=aib(1t/2)/\sqrt{1t}`$, whereas the equation $`\delta =0`$ gives $`(\gamma /t)\sqrt{1t}K+ib/(g+a)=0`$. Eliminating in this equation $`a`$ and $`b`$, one arrives finally after some algebra at the following saddle-point equation

$$\frac{1}{t}=\frac{1}{2}\frac{1}{\gamma K}+\frac{(K+g)/2}{\pm \sqrt{(K+g)^2+1+4/\gamma ^2g^2}2/\gamma },$$
(7)
where the sign “$`+`$” must be taken. Independently on the taken choice of the sign, equation (7) can be further equivalently represented as the following forth order equation in $`K`$:

$`{\displaystyle \frac{K^4+2gK^3}{t}}(1{\displaystyle \frac{1}{t}})`$ $`+`$ $`{\displaystyle \frac{K^2}{4}}\left[(g+{\displaystyle \frac{2}{\gamma }})^2{\displaystyle \frac{8g}{t\gamma }}(1{\displaystyle \frac{2}{t}})^2\right]`$ (8)
$`+`$ $`{\displaystyle \frac{K}{\gamma }}(1{\displaystyle \frac{2}{t}}){\displaystyle \frac{1}{\gamma ^2}}=0.`$
In the zero absorption limit, $`\gamma 0`$ with fixed $`t/\gamma =\tau `$, this equation simplifies to the cubic one found earlier Sommers et al. (2001).

The choice of the sign “$`+`$” is justified by two reasons. First, only then the complex solution of Eq. (7) yields the distribution $`𝒫(r)`$ nonzero at positive $`r`$. Second, the function $`K(t)`$, as follows directly from its definition, must have at $`t\mathrm{}`$ the positive limit $`K_{\mathrm{}}=t_H^1\overline{\text{tr }Q_\gamma }<1`$ Savin and Sommers (2003), which is the mean scaled proper delay time $`\tau `$. The equation for the latter follows then readily from (8) as $`[(g+2/\gamma )^21]\tau ^2+4\tau /\gamma 4/\gamma ^2=0`$ that gives $`\tau =T/(\gamma +T)`$, reproducing thus Eq. (3) by virtue of $`r=1\gamma \tau `$.

The higher moments of the distribution can be easily calculated in the same way, exploiting the $`1/t`$–expansion further. For example, substituting $`K=\tau +\gamma \tau ^2/t+O(t^2)`$ in (8) and equating there the $`1/t`$–term to zero, we find $`\tau ^2=\tau ^2[\gamma ^2+2(\gamma +T)]/(\gamma +T)^2`$. That yields the variance of reflection eigenvalues as follows:

$$r^2r^2=\frac{\gamma ^2T^2(2\gamma +2T(1\gamma )T^2)}{(\gamma +T)^4}.$$
(9)
Remarkably, the variance and, therefore, fluctuations of $`r`$ are suppressed in the both limits of weak and strong absorption.

To understand qualitatively the structure of the solutions of the saddle-point equation, it is instructive to consider the behavior of $`1/t`$ as a function of $`K`$ shown schematically on Fig. 1. We readily see that the distribution function $`𝒫(r)`$ is nonzero only in the finite domain $`r_{\mathrm{min}}<r<r_{\mathrm{max}}`$. The value of the upper border $`r_{\mathrm{max}}<1`$ at any $`g`$ and $`\gamma `$. In a vicinity of $`r_{\mathrm{max}}`$ the distribution behaves as $`𝒫(r)\sqrt{r_{\mathrm{max}}r}`$. For the lower border we find $`r_{\mathrm{min}}=0`$ and $`𝒫(r)1/\sqrt{r}`$, when $`(g2/\gamma )^21`$ or absorption being from the interval $`T\gamma T/(1T)`$, and $`r_{\mathrm{min}}>0`$ with $`𝒫(r)\sqrt{rr_{\mathrm{min}}}`$ otherwise borders . Further analytical study is possible in the following particular cases, which we consider in detail now.

Perfect coupling, $`g=T=1`$. In this case, the region between the two vertical dash-dotted lines on Fig. 1 shrinks to the line $`K=1`$. Eliminating the resulting common factor $`(K+1)`$, one gets from Eq. (8) the following cubic equation

$`K^3+K^2+{\displaystyle \frac{(1r)^2}{r\gamma ^2}}\left[K\left({\displaystyle \frac{2\gamma }{1r}}1\gamma \right)+1\right]=0,`$ (10)
written already in the variable $`r=1t`$. We find that the discriminant $`D=[(1+\gamma )/(3\gamma ^2r)]^3(r_+r)(rr_{})(1r)^4`$ of this equation, with

$$r_\pm =\frac{8+20\gamma ^2\gamma ^4\pm \gamma (8+\gamma ^2)^{3/2}}{8(1+\gamma )^3},$$
(11)
is positive at $`r>0`$ (and thereby $`\mathrm{Im}K`$ is nonzero) only in the domain $`\text{max}[0,r_{}]<r<r_+`$. This sets explicitly the exact borders for the searched distribution. $`𝒫(r)`$ can be found from (10) by applying Cardan’s formulae, reproducing exactly the result of Ref. Beenakker (1998) obtained by a different method. When absorption is weak, $`\gamma <1`$, the behavior of $`𝒫(r)`$ can be well approximated by the following simple interpolation expression interpol :

$$𝒫_{\gamma <1}(r)=\frac{C_\gamma }{2\pi }\frac{\sqrt{(r_+r)(rr_{})}}{(1r)^2\sqrt{r}},$$
(12)
with $`C_\gamma =1+\frac{3}{2}\gamma +O(\gamma ^2)`$ being the normalization constant. In the opposite case of strong absorption, the distribution is found to be close to

$$𝒫_{\gamma >1}(r)=\frac{2\sqrt{1r_+}}{\pi r_+}\frac{\sqrt{r_+/r1}}{(1r)^2}.$$
(13)
The limiting distribution (12) or (13) becomes asymptotically exact as $`\gamma `$ diminishes or grows, respectively. All these features are illustrated on Fig. 2a.

Special coupling, $`g^2=1+4/\gamma ^2`$. Under this condition Eq. (7) reduces to a quadratic equation in $`K`$, giving readily

$$𝒫_g(r)=\frac{2\sqrt{1r_g}}{\pi r_g}\frac{\sqrt{r_g/r1}}{(1r)^2},$$
(14)
with $`r_g=2/(1+\sqrt{1+\gamma ^2/4})=2/(1+g/\sqrt{g^21})`$.

Weak coupling, $`g2/T1`$. One can find approximate expressions for $`𝒫(r)`$, making use of the scaling $`KO(1/g)`$ in the exact equations. When absorption is weak we may expand the square root in (7) and keep only the main contribution $`(2g/\gamma )(1+2gK)^1`$ instead of the last term there. That leads again to a quadratic equation in $`K`$, giving readily the behavior $`𝒫(r)(\gamma /2g\pi ^2)^{1/2}\sqrt{r_{\mathrm{max}}r}/(1r)^2`$ near the upper border $`r_{\mathrm{max}}1\gamma /8g`$. For $`\gamma `$ from the small interval $`[T,T/(1T)]`$, when $`r_{\mathrm{min}}=0`$, $`𝒫(r)`$ is given in the leading approximation by Eq. (14), with $`r_g`$ replaced by $`r_{\mathrm{max}}`$, because the applicability condition for (14) is now satisfied up to $`O(g^2)`$. At other values of small $`\gamma `$ there is no reliable approximation for the lower border available and $`r_{\mathrm{min}}>0`$ is to be found from the general equation borders . At last, for very strong absorption and arbitrary non-perfect coupling, $`T1`$, we can use scaling $`KO(1/\sqrt{\gamma })`$ in (7). That allows us to replace there the last term with $`\frac{1}{2}(K+g)/(1+gK2/\gamma )`$, yielding finally

$$𝒫_{\gamma 1}(r)=\frac{C}{2\pi }\frac{\sqrt{(r_0+\mathrm{\Delta }rr)(rr_0+\mathrm{\Delta }r)}}{(1r)^2(1T+r)},$$
(15)
with the normalization constant $`C\gamma 2(1T)`$ and

$$r_01T\frac{2T(23T)}{\gamma }\mathrm{and}\mathrm{\Delta }r2\sqrt{2}T\sqrt{\frac{1T}{\gamma }},$$

that is valid when $`r_0\mathrm{\Delta }r`$ or $`\gamma 8T^2/(1T)`$. Expression (15) becomes asymptotically exact as $`\gamma `$ grows, approaching the “semi-circle” distribution with the center at $`r_0`$ and the radius $`\mathrm{\Delta }r`$, see Fig. 2b.

Thermal emission.– As an application of our results we consider thermal emission from random media. In his seminal paper Beenakker Beenakker (1998) has shown that the quantum optical problem of the photon statistics can be reduced to a computation of the $`S`$-matrix of the classical wave equation. In particular, chaotic radiation may be characterized by the effective number $`\nu _{\text{eff}}`$ degrees of freedom as follows: $`\nu _{\text{eff}}/\nu =(1r)^2/(1r)^21`$ Beenakker (1998), with $`\nu _{\text{eff}}=\nu `$ for blackbody radiation Mandel and Wolf (1995). We find using (9) that at any $`\gamma `$ and $`T`$

$$\frac{\nu _{\text{eff}}}{\nu }=\frac{\tau ^2}{\tau ^2}=\frac{(\gamma +T)^2}{\gamma ^2+2(\gamma +T)},$$
(16)
and a mean photocount $`\overline{n}=\nu f\gamma T/(\gamma +T)`$, with $`f`$ being a Bose-Einstein function. The earlier result Beenakker (1998) is reproduced at $`T=1`$. Upon the substitution of $`\gamma `$ with $`\gamma `$ one can use Beenakker (1998) $`\overline{n}`$ and $`\nu _{\text{eff}}/\nu `$ (16) even for amplified spontaneous emission below the laser threshold, $`\gamma <T`$. In this case our mean photocount agrees with the findings of Ref. Hackenbroich et al. (2001), where the general theory of photocount statistics in random amplifying media was developed. In the limit of vanishing absorption or amplification, the ratio $`\nu _{\text{eff}}/\nu `$ is $`T/2`$. The large $`\gamma `$ expansion $`\nu _{\text{eff}}/\nu =12(1T)/\gamma +O(\gamma ^2)`$ shows that the saturation to the blackbody limit gets slower when transmission $`T<1`$.

In conclusion, for many-channel chaotic systems we have derived the general distribution of reflection eigenvalues at arbitrary values of absorption and transmission. We note that due to a duality relation Beenakker (1998); Paasschens et al. (1996), an amplifying system in the linear regime ($`\mathrm{\Gamma }_a<\mathrm{\Gamma }_W`$) is directly linked to the dual absorbing one through the change of the sign of $`\mathrm{\Gamma }_a`$ in (1) and correspondingly thereafter. As a result, the reflection matrices (and their eigenvalues) of dual systems are each other’s reciprocal. Therefore, the analysis presented can straightforwardly be extended to the case of linear amplification that might be also relevant for the rapidly developing field of random lasers Beenakker et al. (1996); Beenakker (1998); Hackenbroich et al. (2001); Paasschens et al. (1996); Cao (2003).

We are grateful to G. Hackenbroich and C. Viviescas for useful discussions. The financial support by the SFB/TR 12 der DFG is acknowledged with thanks.


text:
# Quantum dynamics with non-Markovian fluctuating parameters

## I Introduction

The dynamics of classical and quantum systems in the presence of randomly fluctuating micro-fields of environment presents one of the fundamental problems in physics. Spin relaxation in solids anderson53 ; kubo , exciton transport in molecular systems reineker , single-molecular spectroscopy jbs02 , classical transport processes with fluctuating barriers fbarrier  present a few relevant examples. A popular approach consists in modeling the ambient noise influence by means of a classical stochastic field acting on the dynamical system. In the case of quantum systems, such a phenomenological approach is known under the label of Stochastic Liouville Equation (SLE) approach kubo63 ; fox ; lindenberg . It is suitable in the limit of sufficiently high (formally infinite) temperatures reineker ; lindenberg . In the field of chemical kinetics a similar methodology is known under the label of rate processes with dynamical disorder where the rates of chemical reactions fluctuate zwanzig . Moreover, the addition of a non-equilibrium classical noise into dissipative quantum dynamics can serve to describe the influence of the non-equilibrium environmental degrees of freedom on the transport properties gpm95 .

The tractability of the SLE approach, which allows to arrive at the exact model solutions in several known cases, makes it popular over the years. For example, the case of two-state quantum dynamics subjected to a white Gaussian noise can be treated exactly and the corresponding exact master equations for the averaged parameters of quantum systems can be derived reineker ; lindenberg . However, already in the case of a colored Gaussian noise (the so-called colored noise problem) a perturbation theory must be used which leads generally to an approximate description, e.g., within a generalized master equation approach mukamel .

What to do, however, when the ambient noise has a non-Gaussian statistics and/or it does exhibit long-range, e.g., power law temporal correlations with a very large, practically infinite range? A relevant example is given by a $`1/f^\alpha `$ noise which is ubiquitous in the amorphous solids and other glass-like materials like proteins weismann ; lowen . Any perturbation theory in such situations will certainly fail and we are confronting with a rather difficult problem. Nevertheless, for an arbitrary quantum dynamics the problem of finding the corresponding noise-averaged propagator can be solved exactly, at least formally, for a rather general class of non-Markovian jump processes modeled as a continuous-time random walk (CTRW) montroll ; lax ; hughes  within arbitrary, but finite number of states $`N`$. A similar problem has been already considered in several previous works, notably in Ref. kampen79 ; chvosta . These works did not solve, however, the problem at hand for the case of stationary noise-averaging for a multi-state non-Markovian noise, when the evolution of the considered stochastic process starts in the infinite past, and not simultaneously with the evolution of the considered dynamical system. In the case of a non-Markovian noise this problem is not trivial. Physically this is but the most important and relevant case to study.

In this work we utilize the most general description of the discrete state non-Markovian processes of the continuous time random walk (CTRW) type with uncorrelated jumps. Generally, such processes are defined by the set of probability densities $`\psi _{ij}(\tau )`$ for making transitions among the discrete noise states hughes . The noise-averaged quantum propagator is obtained below for this general case. The approaches of Ref. kampen79  and Ref. chvosta  are unified within this general description. Next, the problem of stationary noise averaging is considered. It is shown in a constructive way how to make the stationary noise averaging in the case of factorized probability densities like in Ref. chvosta , but with the time-independent matrix of transition probabilities $`p_{ij}`$. In this respect, the present work is most close to Ref. burstein86 , where a similar approach has been proposed, however the explicit solution for the stationary noise averaged quantum propagator has not been obtained. The corresponding formally exact expression for the Laplace-transformed averaged propagator is found in this work in the explicit form for the first time. This presents the first main result of this work which possesses an ample range of applications. In particular, the noise-averaged propagator of the Kubo-oscillator with a stochastically modulated frequency, as well as the corresponding line shape form are also obtained. This presents the second main result of this work which can be important, e.g., for the single-molecular spectroscopy.

## II Model and Theory

Let us consider an arbitrary quantum system with the Hamilton operator $`\widehat{H}[\xi (t)]`$ which depends parametrically on the stochastic process $`\xi (t)`$ which in turn can acquire randomly in time $`N`$ discrete values $`\xi _i`$. Accordingly, the Hamiltonian $`\widehat{H}[\xi (t)]`$ can take on $`N`$ different random operator values $`\widehat{H}[\xi _i]`$. The discrete stochastic process is assumed to be a non-Markovian renewal process which is fully characterized by the set of probability densities $`\psi _{ij}(\tau )`$ which describe random transitions among the states $`\xi _i`$. Namely, $`\psi _{ij}(\tau )`$ is the probability density for making transition from the state $`j`$ to the state $`i`$. These probability densities must obviously be positive and obey the normalization conditions

$$\underset{i=1}{\overset{N}{}}_0^{\mathrm{}}\psi _{ij}(\tau )𝑑\tau =1,$$
(1)
for all $`j=\overline{1,N}`$. All random transitions are assumed to be mutually independent remark0 . The residence time distribution (RTD) $`\psi _j(\tau )`$ in the state $`j`$ reads obviously

$$\psi _j(\tau )=\underset{i}{}\psi _{ij}(\tau ).$$
(2)
The survival probability $`\mathrm{\Phi }_j(\tau )`$ of the state $`j`$ follows then as

$$\mathrm{\Phi }_j(\tau )=_\tau ^{\mathrm{}}\psi _j(\tau )𝑑\tau .$$
(3)
This is the most general description used in the CTRW theory hughes .

The problem is to average the quantum dynamics in the Liouville space which is characterized by the Liouville-von-Neumann equation

$$\frac{d}{dt}\rho (t)=i[\xi (t)]\rho (t)$$
(4)
for the density operator $`\rho (t)`$ over the realizations of noise $`\xi (t)`$. $`[\xi (t)]`$ in Eq. (4) stands for the quantum Liouville superoperator, $`[\xi (t)]()=\frac{1}{\mathrm{}}[\widehat{H}[\xi (t)],()]`$. In other words, one has to find the noise-averaged propagator

$$U(t)=𝒯\mathrm{exp}[i_{t_0}^{t_0+t}[\xi (\tau )]𝑑\tau ],$$
(5)
where $`𝒯`$ denotes the time-ordering operator. It is the major advance of this work that we obtain the Laplace-transform of propagator (5) in the exact form for arbitrary non-Markovian processes $`\xi (t)`$ of the discussed form. The results of previous research work done within the framework of SLE approach for the discrete Markovian processes kubo  follows as a particular Markovian limit of the developed non-Markovian theory. Moreover, some previous non-Markovian theories, notably that by van Kampen kampen79  and that by Chvosta and Reineker chvosta  are also included as different interpretations of this general formulation.

In particular, the approach of van Kampen is reproduced by introduction of the time-dependent age-specific rates $`k_{ij}(t)`$ like in the renewal theory cox . The probability densities $`\psi _{ij}(\tau )`$ then read kampen79 ; remark1

$$\psi _{ij}(\tau ):=k_{ij}(\tau )\mathrm{exp}[\underset{i}{}_0^\tau k_{ij}(t)𝑑t].$$
(6)
The Markovian case corresponds to $`k_{ij}(\tau )=const`$. Any deviation of $`\psi _{ij}(\tau )`$ from the strictly exponential form which yields a time-dependence of the transition rates $`k_{ij}(\tau )`$ amounts to a non-Markovian behavior remark\_new . Furthermore, the survival probability $`\mathrm{\Phi }_j(\tau )`$ in the state $`j`$ within the time-dependent rate description is given by

$$\mathrm{\Phi }_j(\tau )=\mathrm{exp}[\underset{i=1}{\overset{N}{}}_0^\tau k_{ij}(t)𝑑t]$$
(7)
and Eq. (6) then can be recast as

$$\psi _{ij}(\tau ):=k_{ij}(\tau )\mathrm{\Phi }_j(\tau ).$$
(8)
The introduction of time-dependent rates is one possible way to describe the non-Markovian effects remark1 . It is not unique. For example, Chvosta and Reineker adopted a quite different and more general standpoint chvosta . Namely, they defined

$$\psi _{ij}(\tau ):=p_{ij}(\tau )\psi _j(\tau )$$
(9)
with $`_ip_{ij}(\tau )=1`$. The interpretation is as follows. The process stays in the state $`j`$ during a random time interval characterized by the probability density $`\psi _j(\tau )`$. At the end of this time interval the process makes jump into the state $`i`$ with a generally time-dependent conditional probability $`p_{ij}(\tau )`$. Indeed, any stochastic process of the considered kind can be interpreted in this way. For some particular applications in chvosta  the probability densities $`\psi _j(\tau )`$ were taken strictly exponential and all the non-Markovian effects were assumed to come from the time-dependent transition probabilities $`p_{ij}(\tau )`$. By equating Eq. (8) and Eq. (9) and taking into account that $`\psi _j(\tau ):=d\mathrm{\Phi }_j(\tau )/d\tau `$ it is easy to see that van Kampen approach can be reduced to that of Chvosta and Reineker with the time-dependent transition probabilities

$$p_{ij}(\tau )=\frac{k_{ij}(\tau )}{_ik_{ij}(\tau )}$$
(10)
and with the non-exponential probability densities $`\psi _j(\tau )`$ which follow as $`\psi _j(\tau )=\gamma _j(\tau )\mathrm{exp}[_0^\tau \gamma _j(t)𝑑t]`$ with $`\gamma _j(\tau ):=_ik_{ij}(\tau )`$.

The description of non-Markovian effects with the time-dependent transition probabilities $`p_{ij}(\tau )`$ seems, however, be merely a theoretical device. It appears to be rather difficult (if possible) to obtain $`p_{ij}(\tau )`$ from the sample trajectories of an experimentally observed random process $`\xi (t)`$. In view of Eq. (10) the same is valid for the concept of time-dependent rates. These rates cannot be measured directly from the sample trajectories. On the contrary, the RTD $`\psi _j(\tau )`$ and the time-independent $`p_{ij}`$ (with $`p_{ii}:=0`$ what is assumed in the following) can be routinely deduced from the sample trajectories measured in a single-molecular experiment. This latter description is definitely more advantageous from the practical point of view. From the experimentally well-defined quantities, $`\mathrm{\Phi }_j(\tau )`$ and $`p_{ij}`$, the corresponding time-dependent rates description can be readily found as

$$k_{ij}(\tau )=p_{ij}\frac{d\mathrm{ln}[\mathrm{\Phi }_j(\tau )]}{d\tau }.$$
(11)
Moreover, as it will be shown below namely this description with constant $`p_{ij}`$ in Eq. (9) does provide the consistent way in order to construct the stationary realizations of the non-Markovian process $`\xi (t)`$ and, therefore, in order to find $`U(t)`$ averaged, correspondingly, over the stationary realizations of $`\xi (t)`$. These circumstances give definite advantages of the approach with factorized $`\psi _{ij}(\tau )`$ and time-independent $`p_{ij}`$ as compare with the formulations in kampen79  and chvosta , even though the technical details are quite similar. The correct quantum-mechanical propagator averaged over the stationary realizations of a non-Markovian discrete-state noise with an arbitrary number of states is obtained in this work for the first time. This corresponds to the situation where the quantum system was prepared at the time $`t_0`$ in a non-equilibrium state described by the density matrix $`\rho (t_0)`$ , but the noise was not specially prepared but it has been already relaxed to its stationary state. The ambient noise evolution was started in the infinite past and one assumes that the initial preparation of the quantum system in a non-equilibrium state has no influence on the noise source. This is the most important physical situation to confront with.

The task of performing the noise-averaging of the quantum dynamics in Eq. (5) can be solved exactly due the piecewise constant character of the noise $`\xi (t)`$ burshtein ; frisch . Namely, let us consider the time-interval $`[t_0,t]`$ and to take a frozen realization of $`\xi (t)`$ assuming $`k`$ switching events within this time-interval at the time-instants $`t_i`$,

$$t_0<t_1<t_2<\mathrm{}<t_k<t.$$
(12)
Correspondingly, the noise takes on the values $`\xi _{j_0},\xi _{j_1},\mathrm{},\xi _{j_k}`$ in the time sequel. Then, the propagator $`U(t,t_0)`$ reads obviously

$$U(t,t_0)=e^{i[\xi _{j_k}](tt_k)}e^{i[\xi _{j_{k1}}](t_kt_{k1})}\mathrm{}e^{i[\xi _{j_0}](t_1t_0)}$$
(13)
Let us assume first that it is known with certainty that at the time instant $`t_0`$ the process $`\xi (t)`$ has just started its sojourn in the state $`j_0`$. In other words, the process $`\xi (t)`$ has been prepared in the state $`j_0`$ at $`t_0`$. Then, the corresponding $`k`$times probability density for such noise realization is

$$P_k(\xi _{j_k},t_k;\xi _{j_{k1}},t_{k1};\mathrm{};\xi _{j_1},t_1|\xi _{j_0},t_0)=\mathrm{\Phi }_{j_k}(tt_k)\psi _{j_kj_{k1}}(t_kt_{k1})\mathrm{}\psi _{j_1j_0}(t_1t_0).$$
(14)
In order to obtain the noise-averaged $`U(t,t_0)_{j_0}`$ conditioned on such nonstationary initial noise preparation one has to average (13) with the probability measure in (14) (for $`k=\overline{0,\mathrm{}}`$). Literally (operationally) this means the following. First, one has to construct the time-ordered product of (13) and (14), i.e.

$$\mathrm{\Phi }_{j_k}(tt_k)e^{i[\xi _{j_k}](tt_k)}\psi _{j_kj_{k1}}(t_kt_{k1})e^{i[\xi _{j_{k1}}](t_kt_{k1})}\mathrm{}\psi _{j_1j_0}(t_1t_0)e^{i[\xi _{j_0}](t_1t_0)}.$$
(15)
Second, one has to perform the $`k`$-dimensional time integration of (15) over the variables $`\{t_k\}`$ within the time-ordered domain (12) and to sum the results over all possible $`\{j_k\}`$. Furthermore, this procedure has to be repeated for every $`k=\overline{0,N}`$ and the results summed at the end. The case $`k=0`$ is special with

$$P_0(\xi _{j_0},t_0)=\mathrm{\Phi }_{j_0}(tt_0).$$
(16)
The just outlined task can be easily done formally by use of the Laplace-transform, denoted in the following as $`\stackrel{~}{A}(s):=_0^{\mathrm{}}\mathrm{exp}(s\tau )A(\tau )𝑑\tau `$ for any time-dependent quantity $`A(\tau )`$. For this goal, let us introduce two auxiliary matrix operators $`\stackrel{~}{A}(s)`$ and $`\stackrel{~}{B}(s)`$ with the matrix elements

$$\stackrel{~}{A}_{kl}(s):=\delta _{kl}_0^{\mathrm{}}\mathrm{\Phi }_l(\tau )e^{(s+i[\xi _l])\tau }𝑑\tau $$
(17)
and

$$\stackrel{~}{B}_{kl}(s):=_0^{\mathrm{}}\psi _{kl}(\tau )e^{(s+i[\xi _l])\tau }𝑑\tau ,$$
(18)
correspondingly. Then all what remains to do is to sum the geometric matrix operator series. The result reads,

$$\stackrel{~}{U}(s)_{j_0}=\underset{i}{}\left(\stackrel{~}{A}(s)[I\stackrel{~}{B}(s)]^1\right)_{ij_0},$$
(19)
where $`I`$ is the unity matrix. It is quite obvious that instead of the quantum Liouville operator in Eq. (4) there could be any linear operator, e.g. $``$ be a matrix and $`\rho (t)`$ be then a vector function. Then, the developed theory can be immediately applied to the averaging of arbitrary linear stochastic differential equations frisch ; kampenBook . Even nonlinear stochastic differential equations can be attempted to deal with by introducing a Liouville equation for the corresponding classical probability density kampenBook . The earlier results in kampen79  and chvosta  are reproduced immediately from Eqs. (17)-(19) by an appropriate modification of the considered problem and specifying the transition probability densities $`\psi _{ij}(\tau )`$ from a most general (non-factorized form) to a particular representation in accord with the above discussion.

The derived result in Eqs.(17)-(19) corresponds to the initial preparation of $`\xi (t)`$ in a particular state $`j_0`$. Experimentally, this presents a quite unusual and strongly non-equilibrium situation. For a stationary environment one has to perform yet an additional averaging of $`\stackrel{~}{U}(s)_{j_0}`$ over the initial distribution $`p_{j_0}(t_0)`$ taken as the stationary distribution, i.e., $`p_{j_0}(t_0)=p_{j_0}^{st}`$, where $`p_j^{st}=lim_t\mathrm{}p_j(t)`$. Indeed, this presents a valid prescription for stationary noise-averaging in the Markovian case. However, in a non-Markovian case this prescription is not sufficient.

Quite generally, the stationarity of noise realizations in the strict sense requires kampenBook  that not only the single-time distribution $`p_j(t)`$, but also all the multi-time joint probability distributions of the given process must be invariant of a simultaneous time shift of all time arguments. However, for many physical applications, the stationarity in a weak sense, i.e., on the level of the two-time, $`P(j,t;j_0,t_0)`$, joint distribution is sufficient. Then, e.g., the stationary power spectrum of the corresponding process can be defined. This two-time joint distribution can be expressed as $`P(j,t;j_0,t_0)=\mathrm{\Pi }_{jj_0}(t|t_0)p_{j_0}(t_0)`$ via the conditional probabilities $`\mathrm{\Pi }_{jj_0}(t|t_0)`$ (propagator of the process). For a stationary process the consistency condition, $`p_j^{st}=_{j_0}\mathrm{\Pi }_{jj_0}(t|t_0)p_{j_0}^{st}`$, must be satisfied for all times, i.e. $`p_j^{st}`$ has to be the fixed point of the corresponding propagator (see Appendix A). The propagator of non-Markovian process having this property can be called quasi-stationary. In the present case, in order to construct such propagator and the corresponding stationary realizations of the noise trajectories the probability density of the first time intervals must differ from the all subsequent ones. Indeed, if a noise state $`j`$ was occupied at $`t=t_0`$ with the stationary probability $`p_j^{st}`$, it is not known for how long this state was already occupied before $`t_0`$. The proper conditioning on and averaging over this unknown time must be made and the corresponding survival probabilities for the first residence time interval $`\tau _0=t_1t_0`$ be introduced remark3 . These survival probabilities read cox ; remark3 ; remark4

$$\mathrm{\Phi }_j^{(0)}(\tau )=\frac{_\tau ^{\mathrm{}}\mathrm{\Phi }_j(\tau ^{})𝑑\tau ^{}}{\tau _j},$$
(20)
where $`\tau _j=_0^{\mathrm{}}\mathrm{\Phi }_j(\tau )𝑑\tau `$ is the mean residence time (MRT) of the noise in the state $`j`$. Only for strictly exponential survival probabilities, i.e., in the Markovian case, $`\mathrm{\Phi }_j^{(0)}(\tau )=\mathrm{\Phi }_j(\tau )`$. Otherwise, this is not the case. The corresponding residence time distributions follow immediately as the negative time derivative of $`\mathrm{\Phi }_j^{(0)}(\tau )`$ in Eq. (20) and are known to be cox ; lax

$$\psi _j^{(0)}(\tau )=\frac{\mathrm{\Phi }_j(\tau )}{\tau _j}.$$
(21)
The first time transition densities $`\psi _{ij}^{(0)}(\tau )`$ then follow as

$$\psi _{ij}^{(0)}(\tau )=p_{ij}\frac{\mathrm{\Phi }_j(\tau )}{\tau _j}.$$
(22)
For the logical consistency of this definition with the consideration pursued in remark remark3  $`p_{ij}`$ must be time-independent constants. Otherwise, a logical problem emerges: How to make the proper conditioning of $`p_{ij}(\tau )`$ on the unknown times before $`t_0`$? This is the reason why within the approaches of time-dependent $`p_{ij}(\tau )`$, or time-dependent rates $`k_{ij}(\tau )`$ it is rather obscure how to solve the problem of stationary noise-averaging. Therefore, these approaches do not seem suit well for this stated purpose.

In accord with the above discussion, the transition density $`\psi _{j_1j_0}(t_1t_0)`$ in Eq. (14) must be replaced by $`\psi _{j_1j_0}^{(0)}(t_1t_0)`$ from Eq. (22). Moreover, $`\mathrm{\Phi }_{j_0}(tt_0)`$ in Eq. (16) must be replaced by $`\mathrm{\Phi }_{j_0}^{(0)}(tt_0)`$ from Eq. (20). To account for these modifications, two auxiliary quantities $`\stackrel{~}{A}_{kl}^{(0)}(s)`$ and $`\stackrel{~}{B}_{kl}^{(0)}(s)`$ are introduced which are given by the expressions similar to Eq. (17) and Eq. (18), but with $`\mathrm{\Phi }_j^{(0)}(\tau )`$ instead of $`\mathrm{\Phi }_j(\tau )`$ and $`\psi _{ij}^{(0)}(\tau )`$ instead of $`\psi _{ij}(\tau )`$, i.e.,

$$\stackrel{~}{A}_{kl}^{(0)}(s):=\frac{\delta _{kl}}{\tau _l}_0^{\mathrm{}}e^{(s+i[\xi _l])\tau }_\tau ^{\mathrm{}}\mathrm{\Phi }_l(\tau ^{})𝑑\tau ^{}𝑑\tau $$
(23)
and

$$\stackrel{~}{B}_{kl}^{(0)}(s):=\frac{p_{kl}}{\tau _l}_0^{\mathrm{}}\mathrm{\Phi }_l(\tau )e^{(s+i[\xi _l])\tau }𝑑\tau =\frac{p_{kl}}{\tau _l}\stackrel{~}{A}_{ll}(s).$$
(24)
The resulting geometric operator series can again be easily summed exactly. For the (stationary) noise-averaged Laplace-transformed propagator $`\stackrel{~}{U}(s)`$ we obtain:

$$\stackrel{~}{U}(s)=\underset{ij}{}\left(\stackrel{~}{A}^{(0)}(s)+\stackrel{~}{A}(s)[I\stackrel{~}{B}(s)]^1\stackrel{~}{B}^{(0)}(s)\right)_{ij}p_j^{st},$$
(25)
where

$$p_j^{st}=\frac{\tau _j}{_k\tau _k}$$
(26)
are the stationary occupation probabilities of $`\xi _j`$ (see Eq. (44) in Appendix A and the corresponding discussion). This result can be brought into a physically more insightful form by using the identity $`_\tau ^{\mathrm{}}\mathrm{\Phi }_j(\tau ^{})𝑑\tau ^{}=\tau _j_0^\tau \mathrm{\Phi }_j(\tau ^{})𝑑\tau ^{}`$ and upon introducing two new auxiliary quantities

$$\stackrel{~}{C}_{kl}(s):=\delta _{kl}_0^{\mathrm{}}e^{(s+i[\xi _l])\tau }_0^\tau \mathrm{\Phi }_l(\tau ^{})𝑑\tau ^{}𝑑\tau $$
(27)
and

$$\stackrel{~}{D}_{kl}(s):=\delta _{kl}_0^{\mathrm{}}\psi _l(\tau )e^{(s+i[\xi _l])\tau }𝑑\tau .$$
(28)
Finally we obtain

$$\stackrel{~}{U}(s)=\stackrel{~}{U}(s)_{static}\frac{1}{T}\underset{ij}{}\left(\stackrel{~}{C}(s)\stackrel{~}{A}(s)[IP\stackrel{~}{D}(s)]^1P\stackrel{~}{A}(s)\right)_{ij},$$
(29)
where $`\stackrel{~}{U}(s)_{static}`$ is the Laplace-transform of the statically averaged propagator

$$U(\tau )_{static}:=\underset{k}{}e^{i[\xi _k]\tau }p_k^{st},$$
(30)
$`T=_k\tau _k`$ is the sum of mean residence times, and $`P`$ is the matrix of transition probabilities $`p_{ij}`$ (“scattering matrix” of the random process $`\xi (t)`$). The result in Eqs. (29), (30) together with Eqs. (17), (27), (28), presents the cornerstone result of this work which can be used in numerous applications.

## III Application: Kubo-oscillator

As a simplest practical example we consider the averaging of the so-called Kubo-oscillator

$$\dot{x}(t)=i\omega [\xi (t)]x(t).$$
(31)
This particular problem appears in the theory of optical line shapes, in the nuclear magnetic resonance kubo ; anderson53 , and in the single molecular spectroscopy jbs02 . In Eq. (31), $`\omega [\xi (t)]`$ presents a stochastically modulated frequency of quantum transitions between the levels of a “two-state atom”, or between the eigenstates of a spin 1/2 which are caused by the action of a resonant laser, or magnetic field, respectively. The spectral line shape is determined through the corresponding stochastically averaged propagator of Kubo-oscillator as kubo

$$I(\omega )=\frac{1}{\pi }\underset{ϵ+0}{lim}\mathrm{Re}[\stackrel{~}{U}(i\omega +ϵ)].$$
(32)
Note that the limit $`ϵ+0`$ in Eq. (32) is necessary for the regularization of the corresponding integral in the quasi-static limit $`T\mathrm{}`$. By identifying $`[\xi _k]`$ with $`\omega _k`$ in Eq. (29) we obtain after some algebra

$`\stackrel{~}{U}(s)`$ $`=`$ $`{\displaystyle \underset{k}{}}{\displaystyle \frac{p_k^{st}}{si\omega _k}}{\displaystyle \frac{1}{_k\tau _k}}{\displaystyle \underset{k}{}}{\displaystyle \frac{1\stackrel{~}{\psi }_k(si\omega _k)}{(si\omega _k)^2}}`$
$`+`$ $`{\displaystyle \frac{1}{_k\tau _k}}{\displaystyle \underset{n,l,m}{}}{\displaystyle \frac{1\stackrel{~}{\psi }_l(si\omega _l)}{si\omega _l}}\left({\displaystyle \frac{1}{IP\stackrel{~}{D}(s)}}\right)_{lm}p_{mn}{\displaystyle \frac{1\stackrel{~}{\psi }_n(si\omega _n)}{si\omega _n}}`$
with $`\stackrel{~}{D}_{nm}(s)=\delta _{nm}\stackrel{~}{\psi }_m(si\omega _m)`$. The corresponding line shape follow immediately from Eq. (III) by virtue of Eq. (32). This result presents a non-Markovian generalization of the earlier result by Kubo kubo  for arbitrary $`N`$-state discrete Markovian processes. The generalization consists in allowing for arbitrary non-exponential RTDs $`\psi _k(\tau )`$, or, equivalently, in accordance with Eq. (11) for time-dependent transition rates $`k_{ij}(\tau )`$. This generalization is obtained here for the first time and presents one of our main results. Let us further simplify the result in Eq. (III) for the case of two-state non-Markovian noise with $`p_{12}=p_{21}=1`$. Then, Eq. (III) yields after some simplifications:

$`\stackrel{~}{U}(s)`$ $`=`$ $`{\displaystyle \underset{k=1,2}{}}{\displaystyle \frac{1}{si\omega _k}}{\displaystyle \frac{\tau _k}{\tau _1+\tau _2}}`$ (34)
$`+`$ $`{\displaystyle \frac{(\omega _1\omega _2)^2}{(\tau _1+\tau _2)(si\omega _1)^2(si\omega _2)^2}}{\displaystyle \frac{[1\stackrel{~}{\psi }_1(si\omega _1)][1\stackrel{~}{\psi }_2(si\omega _2)]}{1\stackrel{~}{\psi }_1(si\omega _1)\stackrel{~}{\psi }_2(si\omega _2)}}.`$
With (34) in (32) one obtains the result for the corresponding spectral line shape which is equivalent to one obtained recently in Ref. jung  using a different method. It is reproduced here as a simplest application of our more general approach.

## IV Summary

In this work the problem of the stochastic averaging of a quantum dynamics with non-Markovian fluctuating parameters has been investigated within the trajectory description of continuous time random walk theory. The formally exact expression for the stochastically averaged quantum-mechanical propagator is obtained for the most general CTRW with uncorrelated jumps and for a non-equilibrium noise preparation. The problem of stationary noise averaging has been solved for the practically relevant formulation with the time-independent matrix of transition probabilities $`p_{ij}`$. Especially, the formally exact expression for the stationary averaged quantum propagator has been found in an explicit form. This general expression has been used on order to find the stationary propagator of the Kubo oscillator describing the spectral line shape of optical transitions in a two-state atom. Further applications, such as decoherence of a two-state quantum dynamics driven by two-state non-Markovian noises, including $`1/f^\alpha `$ noise case, are in progress.

###### Acknowledgements.

This work has been supported by the Deutsche Forschungsgemeinschaft via the Sonderforschungsbereich SFB-486, Manipulation of matter on the nanoscale, project No. A10.

## Appendix A Propagators and generalized master equations

In this Appendix, both the nonstationary propagator and the propagator for quasi-stationary initial preparations for the considered non-Markovian processes are obtained, along with the corresponding generalized master equations (GMEs). The propagator $`\mathrm{\Pi }_{ij}`$ of the process $`\xi (t)`$, or the matrix of conditional probabilities connects the initial probability vector $`\stackrel{}{p}(t_0)`$ with the final one, $`\stackrel{}{p}(t_0+\tau )`$, i.e.,

$`p_i(t_0+\tau )={\displaystyle \underset{j}{}}\mathrm{\Pi }_{ij}(t_0+\tau |t_0)p_j(t_0).`$ (35)
The expression for the Laplace-transform $`\stackrel{~}{\mathrm{\Pi }}_{ij}(s)`$ can be obtain in a way similar to the averaging of quantum propagator in Eqs. (19) and (25). Basically, one has to put there $`0`$. For the nonstationary propagator of $`\xi (t)`$, i.e., when the process $`\xi (t)`$ starts its evolution at $`t=t_0`$ in a particular state the result reads for the general case of non-factorized probability densities $`\psi _{ij}(\tau )`$ as follows:

$`\stackrel{~}{\mathrm{\Pi }}_{ij}(s)=\stackrel{~}{\mathrm{\Phi }}_i(s)\left([I\stackrel{~}{\mathrm{\Psi }}(s)]^1\right)_{ij},`$ (36)
where $`\mathrm{\Psi }`$ is the matrix of $`\stackrel{~}{\psi }_{ij}(s)`$ and $`\stackrel{~}{\mathrm{\Phi }}_j(s)=[1\stackrel{~}{\psi }_j(s)]/s`$, $`\stackrel{~}{\psi }_j(s)=_i\stackrel{~}{\psi }_{ij}(s)`$. For quasi-stationary initial preparations the propagator of non-Markovian process is generally different hanggi ; goychuk03 . In the factorized case $`\psi _{ij}(\tau )=p_{ij}\psi _j(\tau )`$ it reads for the considered process (cf. Eq. (25)):

$`\stackrel{~}{\mathrm{\Pi }}_{ij}^{st}(s)=\stackrel{~}{\mathrm{\Phi }}_i^{(0)}(s)\delta _{ij}+\stackrel{~}{\mathrm{\Phi }}_i(s){\displaystyle \underset{k}{}}\left([I\stackrel{~}{\mathrm{\Psi }}(s)]^1\right)_{ik}p_{kj}{\displaystyle \frac{\stackrel{~}{\mathrm{\Phi }}_j(s)}{\tau _j}},`$ (37)
where $`\stackrel{~}{\mathrm{\Phi }}_j^{(0)}(s)=1/s[1\stackrel{~}{\psi }_j(s)]/(s^2\tau _j)`$. The stationary populations follow as $`p_i^{st}=lim_{s0}(s\stackrel{~}{\mathrm{\Pi }}_{ij}(s))`$. The quasi-stationary propagator $`\mathrm{\Pi }_{ij}^{st}(\tau |0)`$ must satisfy the consistency condition $`p_i^{st}=_j\mathrm{\Pi }_{ij}^{st}(\tau |0)p_j^{st}`$ for all times $`\tau `$, i.e., $`\stackrel{}{p}^{st}`$ is the fixed point of $`\mathrm{\Pi }_{ij}^{st}(\tau |0)`$. Let us prove this fact and to find $`p_i^{st}`$. It is more convenient to do both tasks by finding first the corresponding GMEs for $`p_i(t)`$. These generalized master equations are of substantial interest per se.

In order to find the corresponding GMEs the procedure of burstein86  can be applied to a more general present case of non-factorized $`\psi _{ij}(\tau )`$. Indeed, let us consider the conditional probability $`P_k(j|j_0)(t)`$ for making $`k`$-jumps within the time interval $`[0,t]`$ starting at $`t=0`$ in the state $`j_0`$ with the probability $`p_{j_0}(0)`$ and finishing in the state $`j`$ with the probability $`p_j(t)`$. This probability is given by a corresponding $`k`$-dimensional integral of (14)) (see the discussion below Eq. (14)) with the summation made over $`j_{k1},j_{k2},\mathrm{},j_1`$. The corresponding Laplace-transform $`\stackrel{~}{P}_{jj_0}^{(k)}(s)`$ reads:

$`\stackrel{~}{P}_{jj_0}^{(k)}(s)={\displaystyle \underset{j_{k1}}{}}\mathrm{}{\displaystyle \underset{j_1}{}}\stackrel{~}{\mathrm{\Phi }}_j(s)\stackrel{~}{\psi }_{jj_{k1}}(s)\mathrm{}\stackrel{~}{\psi }_{j_2j_1}(s)\stackrel{~}{\psi }_{j_1j_0}^{(0)}(s)`$ (38)
for $`k>0`$. Here, Eq. (14) was used in a slightly modified form with $`\psi _{ij}^{(0)}(\tau )`$ for the probability densities of the first time intervals. Furthermore, for $`k=0`$, $`\stackrel{~}{P}_{jj_0}^{(0)}(s)=\stackrel{~}{\mathrm{\Phi }}_j^{(0)}(s)\delta _{jj_0}`$, where $`\stackrel{~}{\mathrm{\Phi }}_j^{(0)}(s)`$ is the Laplace-transform of the corresponding survival probability $`\mathrm{\Phi }_j^{(0)}(\tau )=_i_\tau ^{\mathrm{}}\psi _{ij}^{(0)}(\tau ^{})𝑑\tau ^{}`$ of the first time-interval. For $`k2`$, The quantities $`\stackrel{~}{P}_{jj_0}^{(k)}(s)`$ satisfy obviously the following recurrence relation:

$`\stackrel{~}{P}_{jj_0}^{(k)}(s)=\stackrel{~}{\mathrm{\Phi }}_j(s){\displaystyle \underset{n}{}}\stackrel{~}{\psi }_{jn}(s){\displaystyle \frac{\stackrel{~}{P}_{nj_0}^{(k1)}(s)}{\stackrel{~}{\mathrm{\Phi }}_n(s)}}.`$ (39)
Furthermore, the Laplace-transform of propagator $`\stackrel{~}{\mathrm{\Pi }}_{jj_0}(s)`$ is expressed in terms of $`\stackrel{~}{P}_{jj_0}^{(k)}(s)`$ as $`\stackrel{~}{\mathrm{\Pi }}_{jj_0}(s)=_{k=0}^{\mathrm{}}\stackrel{~}{P}_{jj_0}^{(k)}(s)`$. Then, by virtue of Laplace-transformed Eq.(35) with $`t_0=0`$,

$`{\displaystyle \underset{k=2}{\overset{\mathrm{}}{}}}{\displaystyle \underset{j_0}{}}\stackrel{~}{P}_{jj_0}^{(k)}(s)p_{j_0}(0)`$ $`=`$ $`\stackrel{~}{p}_j(s){\displaystyle \underset{j_0}{}}\stackrel{~}{P}_{jj_0}^{(1)}(s)p_{j_0}(0)\stackrel{~}{P}_{jj}^{(0)}(s)p_j(0)`$ (40)
$`=`$ $`\stackrel{~}{\mathrm{\Phi }}_j(s){\displaystyle \underset{n}{}}\stackrel{~}{\psi }_{jn}(s){\displaystyle \frac{1}{\stackrel{~}{\mathrm{\Phi }}_n(s)}}{\displaystyle \underset{k=2}{\overset{\mathrm{}}{}}}{\displaystyle \underset{j_0}{}}\stackrel{~}{P}_{nj_0}^{(k1)}(s)p_{j_0}(0)`$
where the recurrence relation (39) has been used. The use of $`_{k=2}^{\mathrm{}}_{j_0}\stackrel{~}{P}_{nj_0}^{(k1)}(s)p_{j_0}(0)=\stackrel{~}{p}_n(s)\stackrel{~}{\mathrm{\Phi }}_n^{(0)}(s)p_n(0)`$ in Eq. (40) finally yields

$`\stackrel{~}{p}_j(s)=\stackrel{~}{\mathrm{\Phi }}_j(s){\displaystyle \underset{n}{}}\stackrel{~}{\psi }_{jn}(s){\displaystyle \frac{\stackrel{~}{p}_n(s)}{\stackrel{~}{\mathrm{\Phi }}_n(s)}}+\stackrel{~}{\mathrm{\Phi }}_j^{(0)}(s)p_j(0)`$
$`+\stackrel{~}{\mathrm{\Phi }}_j(s){\displaystyle \underset{n}{}}\left(\stackrel{~}{\psi }_{jn}^{(0)}(s)\stackrel{~}{\psi }_{jn}(s){\displaystyle \frac{\stackrel{~}{\mathrm{\Phi }}_n^{(0)}(s)}{\stackrel{~}{\mathrm{\Phi }}_n(s)}}\right)p_n(0).`$ (41)
Let us consider now the case when the noise has been prepared at $`t_0=0`$ in a particular state with the probability one. Then, $`\psi _{ij}^{(0)}(\tau )\psi _{ij}(\tau )`$ and the last term in (A) vanishes. For this class of non-equilibrium initial preparations, the inversion of (A) yields:

$`\dot{p}_j(t)={\displaystyle \underset{n}{}}{\displaystyle _0^t}\mathrm{\Gamma }_{nj}(tt^{})p_j(t^{})𝑑t^{}+{\displaystyle \underset{n}{}}{\displaystyle _0^t}\mathrm{\Gamma }_{jn}(tt^{})p_n(t^{})𝑑t^{}`$ (42)
where the Laplace-transformed kernels reads

$`\stackrel{~}{\mathrm{\Gamma }}_{jn}(s)={\displaystyle \frac{s\stackrel{~}{\psi }_{jn}(s)}{1\stackrel{~}{\psi }_n(s)}}`$ (43)
with $`\stackrel{~}{\psi }_n(s)=_j\stackrel{~}{\psi }_{jn}(s)`$. The just derived GME (42),(43) is the most general GME for the continuous time random walk processes with uncorrelated jumps for the given class of initial preparations. In the case of factorized (but still non-separable) CTRW with $`\psi _{ij}(\tau )=p_{ij}\psi _j(\tau )`$ it reduces to the GME of Burstein, Zharikov and Temkin burstein86 . Moreover, in assumption that all $`\psi _j(\tau )`$ are equal (separable CTRW of Montroll and Weiss), it reduces further to the GME of Kenkre, Montroll and Shlesinger kenkre . The stationary populations can be obtained from (A) as $`p_j^{st}=lim_{s0}(s\stackrel{~}{p}_j(s))`$. Assuming that the mean residence times $`\tau _j`$ exist, i.e. $`\stackrel{~}{\psi }_{ij}(s)=\alpha _{ij}st_{ij}+o(s)`$ with $`_i\alpha _{ij}=1`$ and $`_it_{ij}=\tau _j`$, Eq. (A) yields the system of linear algebraic equations for the stationary populations

$`{\displaystyle \frac{p_j^{st}}{\tau _j}}={\displaystyle \underset{n}{}}\alpha _{jn}{\displaystyle \frac{p_n^{st}}{\tau _n}}`$ (44)
For an ergodic process the stationary probability to find the process in a particular state should be proportional to the time which the process spends in this particular state on average, i.e., be given by Eq. (26). It is easy to verify that Eq. (26) provides the solution of (44) for $`_j\alpha _{ij}=1`$ (for the factorized case $`\alpha _{ij}`$ coincide obviously with $`p_{ij}`$). This latter condition can be also expressed as

$`{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle _0^{\mathrm{}}}\psi _{ij}(\tau )𝑑\tau =1.`$ (45)
For factorized case $`\psi _{jn}(\tau )=p_{jn}\psi _n(\tau )`$ with $`\psi _j^{(0)}(\tau )=\mathrm{\Phi }_j(\tau )/\tau _j`$ remark3 , the inversion of (A) yields the GME for the discussed non-Markovian process with quasi-stationary initial preparations:

$`\dot{p}_j(t)={\displaystyle _0^t}\mathrm{\Gamma }_j(tt^{})[p_j(t^{})p_j(0)]𝑑t^{}+{\displaystyle \underset{n}{}}p_{jn}{\displaystyle _0^t}\mathrm{\Gamma }_n(tt^{})[p_n(t^{})p_n(0)]𝑑t^{}`$
$`{\displaystyle \frac{p_j(0)}{\tau _j}}+{\displaystyle \underset{n}{}}p_{jn}{\displaystyle \frac{p_n(0)}{\tau _n}}`$ (46)
with the kernels given by $`\stackrel{~}{\mathrm{\Gamma }}_j(s)=s\stackrel{~}{\psi }_j(s)/(1\stackrel{~}{\psi }_j(s))`$. By choosing $`p_n(0)=p_n^{st}`$ which satisfy Eq. (44) it becomes obvious that $`p_n(t)=p_n^{st}`$ is the solution of (A) for all times $`t>0`$. This means that the stationary $`\stackrel{}{p}^{st}`$ provides indeed the fixed point of the corresponding quasi-stationary propagator.