text:
# Local scale invariance as dynamical space-time symmetry in phase-ordering kinetics

## Abstract

The scaling of the spatio-temporal response of coarsening systems is studied through simulations of the $`2D`$ and $`3D`$ Ising model with Glauber dynamics. The scaling functions agree with the prediction of local scale invariance, extending dynamical scaling to a space-time dynamical symmetry.

Ageing phenomena are observed in a broad variety of systems with slow relaxation dynamics. Useful insight may often be gained via the consideration of simple ferromagnetic models (rather than genuinely glassy systems), quenched to a temperature $`T<T_c`$ from a disordered initial state into the low-temperature phase. For their study two-time quantities, rather than the more usual one-time quantities, are particularly useful, see  for recent reviews. Those display dynamical self-similarity in the ageing regime where the order-parameter autocorrelator $`C(t,s)=\varphi (t)\varphi (s)`$ decays from its plateau value $`q_{\mathrm{EA}}=M_{\mathrm{eq}}^2`$ to zero ($`M_{\mathrm{eq}}`$ is the spontaneous magnetization). In this regime both $`s`$ and $`\tau =ts`$ are much larger than the microscopic time scale (set to unity). Similarly, the autoresponse function, defined as $`R(t,s)=\delta \varphi (t)/\delta h(s)`$ with $`t>s`$ shows a scaling behaviour such that the scaling laws

$$C(t,s)M_{\mathrm{eq}}^2f_C(t/s),R(t,s)s^{1a}f_R(t/s)$$
(1)
with $`f_C(1)=1`$ are found to hold for a broad range of models . For $`x1`$, the scaling functions usually fall off as $`f_C(x)x^{\lambda _C/z}`$ and $`f_R(x)x^{\lambda _R/z}`$, where $`z`$ is the dynamic critical exponent, $`\lambda _C`$ is the autocorrelation exponent  and $`\lambda _R`$ is the autoresponse exponent . We shall focus here on phase-ordering in the Glauber-Ising model in $`d>1`$ dimensions where $`a=1/z=1/2`$ is expected . Recent arguments  leading to $`a=1/4`$ in the $`2D`$ Glauber-Ising model have been rejected through a detailed study of the scaling of $`R(t,s)`$, the results of which reconfirm $`a=1/2`$ .

Eq. (1) states that the two-time quantities evaluated at the same spatial location transform covariantly under a global rescaling of time $`tbt`$ with $`b`$ constant. Recently, it has been proposed that the response functions should transform covariantly under local scale transformations with $`b=b(t)`$, but with time-translations excluded . By analogy with conformal invariance, in particular covariance of $`R`$ under the so-called ‘special’ transformations which transform time as $`tt/(1+\gamma t)`$, is assumed. If that is the case, the exact autoresponse scaling function becomes

$$f_R(x)=r_0x^{1+a\lambda _R/z}(x1)^{1a}$$
(2)
and depends on the universality classes only through the values of the exponents $`a`$ and $`\lambda _R/z`$, while $`r_0`$ merely is a normalization constant. This prediction of local scale invariance (LSI) has been confirmed in a variety of physically very different models, see  and references therein. Still, better evidence in favour of LSI than mere phenomenology would be desirable.

The origin of LSI can be understood in the special case $`z=2`$ (which is realized in all cases of phase-ordering kinetics we are going to consider ). A Ward identity can be written down such that if the system is known to be (i) scale-invariant with $`z=2`$ and (ii) Galilei-invariant, then its invariance under the full group of local scale transformations follows . Therefore, scale invariance and Galilei invariance appear as the building blocks for LSI. So far, no direct test of Galilei invariance in phase-ordering systems was ever performed.

Direct tests of Galilei invariance require the study of the time- and space-dependence of response functions. For phase-ordering, it is natural to expect a scaling behaviour of the linear response function

$$R(t,s;𝐫)=\frac{\delta \varphi _𝐫(t)}{\delta h_\mathrm{𝟎}(s)}|_{h=0}=s^{1a}F_R(\frac{t}{s},\frac{𝐫}{(ts)^{1/z}})$$

with the scaling function $`F_R(x,𝐮)`$ and $`F_R(x,\mathrm{𝟎})=f_R(x)`$. For $`z=2`$, Galilei invariance (combined with the Ward identity of LSI) predicts

$$R(t,s;𝐫)=R(t,s)\mathrm{exp}\left(\frac{}{2}\frac{𝐫^2}{ts}\right)$$
(3)
where $`R(t,s)`$ is the autoresponse function given by eqs. (1,2) and $``$ is a direction-dependent non-universal constant. Eq. (3) gives the full spatio-temporal scaling of the linear response. We shall sucessfully test it in the $`2D`$ and $`3D`$ kinetic Glauber-Ising models and provide thereby the first direct evidence of Galilei invariance in a phase-ordering system. Since the $`2D`$/$`3D`$ Glauber-Ising model cannot be reduced to a free-field theory, it is highly non-trivial that its exact response function takes the simple free-field form eq. (3).

As precise simulational data for the autoresponse $`R(t,s)`$ are difficult to obtain, it is convenient to study instead the thermoremanent magnetization

$$TM_{\mathrm{TRM}}(t,s)/h_{(0)}=\rho (t,s)=T_0^sduR(t,u)$$
(4)
We consider the Ising model on a hypercubic lattice, with periodic boundary conditions and the Hamiltonian $`=_{(𝐢,𝐣)}\sigma _𝐢\sigma _𝐣`$. We use heat-bath dynamics defined through the stochastic rule

$$\sigma _𝐢(t+1)=\pm 1\text{with probability }\frac{1}{2}\left[1\pm \mathrm{tanh}(h_𝐢(t)/T)\right]$$
(5)
with the local field $`h_𝐢(t)=_{𝐚(𝐢)}\sigma _𝐚(t)`$ and where $`𝐚(𝐢)`$ runs over the nearest neighbours of the site $`𝐢`$. $`M_{\mathrm{TRM}}(t,s)`$ is measured by applying a quenched spatially random magnetic field $`\pm h_{(0)}`$ for times between the quench at $`t=0`$ and the waiting time $`s`$ . The presence of this external field then changes the local field in eq. (5) to $`h_𝐢(t)=_{𝐚(𝐢)}\sigma _𝐚(t)\pm h_{(0)}`$.

It has been understood recently that the scaling behaviour of $`M_{\mathrm{TRM}}`$ for $`s1`$ is not a simple power law but rather shows a cross-over behaviour . It is of the form

$$\rho (t,s)=r_0s^af_M(t/s)+r_1s^{\lambda _R/z}g_M(t/s)$$
(6)
provided the system was initially prepared at infinite temperature. Here the correction term can be estimated as $`g_M(x)x^{\lambda _R/z}`$, while the scaling function $`f_M(x)`$ can be found from local scale invariance using eqs. (1,2,4) and is given explicitly by eq. (5.47) in . Finally, $`r_{0,1}`$ are normalization constants. In practice, it turns out that $`\lambda _R`$ is quite close to its lower bound $`d/2`$ and in particular in $`2D`$, the correction to the leading scaling behaviour is sizeable. Before any meaningful study of the spatio-temporal response can be carried out, the correction term must be subtracted off.

As explained in , this can be carried out by fixing a value of $`x=t/s`$. Then $`r_{0,1}`$ are obtained by fitting eq. (6) to the computed thermoremanent magnetization $`M_{\mathrm{TRM}}`$. In order to illustrate the quality of this procedure, we show in figure 1 data for the scaling function $`f_M(x)`$ obtained after subtraction of the correction term $`r_1s^{\lambda _R/z}g_M(t/s)`$ from the integrated autoresponse $`\rho `$ of the $`2D`$ and $`3D`$ Glauber-Ising model. Statistical error bars in $`2D`$ are smaller than the symbol size and in $`3D`$ are of the order of the scatter in the data. In table I, we list the values of $`\lambda _R`$  and the constants $`r_{0,1}`$ which were determined at $`x=7`$. We find a nice data collapse and a clear agreement with the LSI prediction.

Having thus checked the expected scaling of the autoresponse function, we can now turn towards the spatio-temporal response. Again, we consider the integrated response rather than $`R(t,s;𝐫)`$ because it is considerably less affected by noise. In order to fix the non-universal parameter $``$ in (3), we form

$$\rho _0(t,s;𝐫)=T_0^sd\tau R(t,s\tau ;𝐫)$$
(7)
which is measured on the lattice by keeping a small random magnetic field $`\pm h_{(0)}`$ until the waiting time $`s`$ at the site $`\mathrm{𝟎}`$ and then observing the thermoremanent magnetization $`M_{\mathrm{TRM}}(t,s;𝐫)=h_{(0)}\rho _0(t,s;𝐫)/T=\overline{h_\mathrm{𝟎}\sigma _𝐫(t)}/T`$ at a different site $`𝐫`$ (where $`\mathrm{}`$ indicates an average over the thermal noise whereas the bar means an average over the random field distribution). We computed the spatially and temporally integrated response for the Ising model with Glauber dynamics in two and three dimensions. The two-dimensional systems usually contained $`300^2`$ spins on a square lattice, whereas in three dimensions cubic lattices with typically $`60^3`$ sites were considered. Some Monte Carlo simulations were also done for other system sizes in order to check against finite-size effects. Since $`\rho _0(t,s;𝐫)`$ is very noisy, one must average over a large number of runs with different realizations of the thermal noise and of the spatially random magnetic field. For every waiting time we avaraged over at least 50000 different runs. The whole study took about $`210^5`$ CPU hours on a SGI Origin 3800 parallel computer.

Following the same lines as in , we arrive at the scaling form, where $`𝐫`$ varies along a fixed direction on the lattice

$`\rho _0(t,s;𝐫)`$ $``$ $`\rho _{\mathrm{}}(ts;𝐫)+r_0s^aF_0({\displaystyle \frac{t}{s}},{\displaystyle \frac{r^2}{s}})`$ (9)
$`+r_1s^{\lambda _R/z}G_0({\displaystyle \frac{t}{s}},{\displaystyle \frac{r^2}{s}})`$
where the universal scaling functions read

$`F_0(x,y)`$ $`=`$ $`{\displaystyle _0^1}dv\mathrm{exp}[y/2(x1+v)]h_R(x,v)`$ (10)
$`G_0(x,y)`$ $``$ $`x^{\lambda _R/z}e^{y/2x}`$ (11)
with $`h_R(x,v):=f_R\left(\frac{x}{1v}\right)(1v)^{1a}`$ where $``$ is a direction-dependent new parameter. In this work we choose $`T`$ such as to keep anisotropy effects small and also to avoid the cross-over to critical dynamics at $`T=T_c`$. This is illustrated in figure 2, for the $`2D`$ case at $`T=1.5`$ where we plot $`F_0`$ as function of $`r^2/s`$. For the displayed data points $`r`$ varies between 1 and $`\sqrt{s}`$ when going along the (10) direction and between $`\sqrt{2}`$ and $`\sqrt{2s}`$ when going along the (11) direction. At the considered temperature any anisotropies are very small and the values of $``$, determined separately in the two directions, coincide, $`_{(10)}=_{(11)}=`$. Similar results also hold in $`3D`$ at $`T=3`$. We fix $``$ at $`x=7`$ (see table I and inset) using (9) and then obtain a parameter-free prediction for other values of $`x`$. This is shown in the inset of figure 2, where the LSI-prediction is compared to numerical data for two additional values of $`x`$ and two different waiting times. Since by now all non-universal parameters are fixed, this constitutes by itself a quantitative test of local scale invariance.

A fuller test of the spatio-temporal response is obtained by considering the spatially and temporally integrated response function

$$\frac{\mathrm{d}\rho (t,s;\mu )}{\mathrm{d}\mathrm{\Omega }}=T_0^sdu_0^{\sqrt{\mu s}}drr^{d1}R(t,u;𝐫)$$
(12)
where the space integral is along a straight line of length $`\mathrm{\Lambda }=\sqrt{\mu s}`$ but we do allow for the possibility of an anisotropy as a function of the solid angle $`\mathrm{\Omega }`$. Such anisotropies are known to occur if $`T<T_c`$ . As before, we derive from (3) the scaling form

$`{\displaystyle \frac{\mathrm{d}\rho (t,s;\mu )}{\mathrm{d}\mathrm{\Omega }}}=\rho ^{(1)}(ts,𝐫^2)`$ (13)
$`+`$ $`r_0s^{d/2a}\rho ^{(2)}(t/s,\mu )+r_1s^{d/2\lambda _R/z}\rho ^{(3)}(t/s,\mu )`$ (14)
with the explicit scaling functions

$`\rho ^{(2)}(x,\mu )`$ $`=`$ $`{\displaystyle \frac{\mu ^{d/2}}{d}}{\displaystyle _0^1}dvh_R(x,v)_d\left({\displaystyle \frac{\mu }{x1+v}}\right)`$ (15)
$`\rho ^{(3)}(x,\mu )`$ $``$ $`{\displaystyle \frac{\mu ^{d/2}}{d}}x^{\lambda _R/z}_d\left({\displaystyle \frac{\mu }{x}}\right)`$ (16)
$`_d(y)`$ $`=`$ $`e^{y/2}{}_{1}{}^{}\text{F}_{1}^{}(1,1+{\displaystyle \frac{d}{2}};{\displaystyle \frac{y}{2}})`$ (17)
This is the general expression for the scaling of the spatio-temporally integrated response function. If we fix $`\mu `$ and let $`x=t/s`$ vary, the form of the scaling function of $`s^{ad/2}\mathrm{d}\rho /\mathrm{d}\mathrm{\Omega }`$ merely depends on $`\mu `$. We stress that the exponent $`\lambda _R`$ and the free parameters $`r_0,r_1,`$ are now all fixed such that there remains no free fitting parameter at all when comparing (13) with our numerical data.

It is of interest to compare the maximal distance $`r_{\mathrm{max}}`$ accessed by our simulation with the physical length scale of the problem, for example the domain size $`L(s)`$ at the instant the magnetic field is turned off. From the correlation function, in $`2D`$ we estimate $`L(s)3,6,9.5`$ for $`s=25,100,225`$, respectively, with $`h_{(0)}=0.05`$. In $`3D`$, we find $`L(s)2.1,3.3,4`$ for $`s=25,64,100`$. As expected $`L(s)s^{1/z}`$, with $`z2`$. The value $`\mu =4`$ ($`\mu =2`$) corresponds to $`r_{\mathrm{max}}/L(s)3.2`$ (3.4) in $`2D`$ ($`3D`$), respectively and we see that the data probe the large-distance region, well beyond the domain size $`L(s)`$.

The integrated spatio-temporal response in $`2D`$ at $`T=1.5`$ is shown in figure 3 for two values of $`\mu `$ and where the correction $`\rho ^{(3)}`$ was already subtracted off. We find a nice scaling behaviour over the whole range of waiting times we could consider. We stress that the agreement between local scale invariance (3) and our data for several values of $`\mu `$ is a real test, since no free parameter remains. In particular, both the height and the position of the maximum of the scaling function for $`\mu =2`$ and $`4`$ is completely fixed. We can conclude that LSI is fully vindicated. This is the first time that the exact functional form of the scaling function of the spatio-temporal response of a generic non-equilibrium spin system is found.

Similarly, the integrated spatio-temporal response in $`3D`$ at $`T=3`$ is shown in figure 4. Again, a nice scaling behaviour in perfect agreement with LSI is found.

The prediction (3) is also verified in a few exactly solvable models with a non-conserved order parameter and undergoing phase-ordering kinetics with $`z=2`$. In particular, (3) holds in the $`d`$-dimensional kinetic spherical model with $`d>2`$ quenched to either $`T<T_c`$ or $`T=T_c`$ and independently of initial correlations .

Summarizing, we presented the first quantitative study of the scaling of the spatio-temporal response in coarsening systems. For the first time, we find direct evidence for Galilei invariance in the ageing regime of a phase-ordering system. This is a new dynamic symmetry whose presence in a non-equilibrium phase transition was not anticipated. Our result provides a strong indication that LSI is indeed a true space-time dynamical symmetry of statistical systems undergoing phase-ordering kinetics.

We thank B.P. Vollmayr-Lee and C. Godrèche for discussions. This work was supported by the Bayerisch-Französisches Hochschulzentrum (BFHZ) and by CINES Montpellier (projet pmn2095).


format: text

text:
# Parrondo’s games as a discrete ratchet

## Abstract

We write the master equation describing the Parrondo’s games as a consistent discretization of the Fokker–Planck equation for an overdamped Brownian particle describing a ratchet. Our expressions, besides giving further insight on the relation between ratchets and Parrondo’s games, allow us to precisely relate the games probabilities and the ratchet potential such that periodic potentials correspond to fair games and winning games produce a tilted potential.

Parrondo paradox, ratchets, master equations

The Parrondo’s paradoxha99 ; ha02  shows that the alternation of two losing games can lead to a winning game. This surprising result is nothing but the translation into the framework of very simple gambling games of the ratchet effectr02 . In particular, the flashing ratchetab94 ; pcpa94  can sustain a particle flux by alternating two relaxational potential dynamics, none of which produces any net flux. Despite that this qualitative relation between the Parrondo paradox and the flashing ratchet has been recognized from the very beginning (and, in fact, it constituted the source of inspiration for deriving the paradoxical games), only very recently there has been some interest in deriving exact relations between bothaa03 ; hkk02 .

In this paper, we rewrite the master equation describing the evolution of the probabilities of the different outcomes of the games in a way that shows clearly its relation with the Fokker–Planck equation for the flashing ratchet. In this way, we are able to give an expression for the dynamical potentials in terms of the probabilities defining the games, as well as an expression for the current. Similarly, given a ratchet potential we are able to construct the games that correspond to that potential.

The Parrondo’s paradox considers a player that tosses different coins such that a unit of “capital” is won (lost) if heads (tails) show up. Although several possibilities have been proposedpha00 ; t01 ; t03 ; na02 ; blp02 ; m00 ; mb02  , in this paper we consider the original and easiest version in which the probability of winning, $`p_i`$, depends on the actual value of the capital, $`i`$, modulus a given number $`L`$. A game is then completely specified by giving the set or probabilities $`\{p_0,p_1,\mathrm{},p_{L1}\}`$ from which any other value $`p_k`$ can be derived as $`p_k=p_{k\mathrm{mod}L}`$. A fair game, one in which gains and losses average out, is obtained if $`_{i=0}^{L1}p_i=_{i=0}^{L1}(1p_i)`$. The paradox shows that the alternation (either random or periodic) of two fair games can yield a winning game. For instance, the alternation of game A defined by $`p_ip=1/2,i`$, and game B defined by $`L=3`$ and $`p_0=1/10,p_1=p_2=3/4`$ produces a winning game although both $`A`$ and $`B`$ are fair games.

A discrete time $`\tau `$ can be introduced by considering that every coin toss increases $`\tau `$ by one. If we denote by $`P_i(\tau )`$ the probability that at time $`\tau `$ the capital is equal to $`i`$, we can write the general master equation

$$P_i(\tau +1)=a_1^iP_{i1}(\tau )+a_0^iP_i(\tau )+a_1^iP_{i+1}(\tau )$$
(1)
where $`a_1^i`$ is the probability of winning when the capital is $`i1`$, $`a_1^i`$ is the probability of losing when the capital is $`i+1`$, and, for completeness, we have introduced $`a_0^i`$ as the probability that the capital $`i`$ remains unchanged (a possibility not considered in the original Parrondo games). Note that, in accordance with the rules described before, we have taken that the probabilities $`\{a_1^i,a_0^i,a_1^i\}`$ do not depend on time. It is clear that they satisfy:

$$a_1^{i+1}+a_0^i+a_1^{i1}=1$$
(2)
which ensures the conservation of probability: $`_iP_i(\tau +1)=_iP_i(\tau )`$.

It is a matter of straightforward algebra to write the master equation in the form of a continuity equation:

$$P_i(\tau +1)P_i(\tau )=\left[J_{i+1}(t)J_i(t)\right]$$
(3)
where the current $`J_i(\tau )`$ is given by:

$$J_i(\tau )=\frac{1}{2}\left[F_iP_i(\tau )+F_{i1}P_{i1}(\tau )\right]\left[D_iP_i(\tau )D_{i1}P_{i1}(\tau )\right]$$
(4)
and

$$F_i=a_1^{i+1}a_1^{i1},D_i=\frac{1}{2}(a_1^{i+1}+a_1^{i1})$$
(5)
This form is a consistent discretization of the Fokker–Plank equationhl84  for a probability $`P(x,t)`$

$$\frac{P(x,t)}{t}=\frac{J(x,t)}{x}$$
(6)
with a current

$$J(x,t)=F(x)P(x,t)\frac{[D(x)P(x,t)]}{x}$$
(7)
with general drift, $`F(x)`$, and diffusion, $`D(x)`$. If $`\mathrm{\Delta }t`$ and $`\mathrm{\Delta }x`$ are, respectively, the time and space discretization steps, such that $`x=i\mathrm{\Delta }x`$ and $`t=\tau \mathrm{\Delta }t`$, it is clear the identification

$$F_i\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}F(i\mathrm{\Delta }x),D_i\frac{\mathrm{\Delta }t}{(\mathrm{\Delta }x)^2}D(i\mathrm{\Delta }x)$$
(8)
The discrete and continuum probabilities are related by $`P_i(\tau )P(i\mathrm{\Delta }x,\tau \mathrm{\Delta }t)\mathrm{\Delta }x`$ and the continuum limit can be taken by considering that $`M=\underset{\mathrm{\Delta }t0,\mathrm{\Delta }x0}{lim}{\displaystyle \frac{(\mathrm{\Delta }x)^2}{\mathrm{\Delta }t}}`$ is a finite number. In this case $`F_iM^1\mathrm{\Delta }xF(i\mathrm{\Delta }x)`$ and $`D_iM^1D(i\mathrm{\Delta }x)`$.

From now on, we consider the case $`a_0^i=0`$. Since $`p_i=a_1^{i+1}`$ we have

$$D_iD=1/2F_i=1+2p_i$$
(9)
and the current $`J_i(\tau )=(1p_i)P_i(\tau )+p_{i1}P_{i1}(\tau )`$ is nothing but the probability flux from $`i1`$ to $`i`$.

The stationary solutions $`P_i^{st}`$ can be found solving the recurrence relation derived from (4) for a constant current $`J_i=J`$ with the boundary condition $`P_i^{st}=P_{i+L}^{st}`$:

$$P_i^{st}=N\mathrm{e}^{V_i/D}\left[1\frac{2J}{N}\underset{j=1}{\overset{i}{}}\frac{\mathrm{e}^{V_j/D}}{1F_j}\right]$$
(10)
with a current

$$J=N\frac{\mathrm{e}^{V_L/D}1}{2_{j=1}^L\frac{\mathrm{e}^{V_j/D}}{1F_j}}$$
(11)
$`N`$ is the normalization constant obtained from $`_{i=0}^{L1}P_i^{st}=1`$. In these expressions we have introduced the potential $`V_i`$ in terms of the probabilities of the games111In this, as well as in other similar expressions, the notation is such that $`_{j=1}^0=0`$. Therefore the potential is arbitrarily rescaled such that $`V_0=0`$.

$$V_i=D\underset{j=1}{\overset{i}{}}\mathrm{ln}\left[\frac{1+F_{j1}}{1F_j}\right]=D\underset{j=1}{\overset{i}{}}\mathrm{ln}\left[\frac{p_{j1}}{1p_j}\right]$$
(12)
The case of zero current $`J=0`$, implies a periodic potential $`V_L=V_0=0`$. This reproduces again the condition $`_{i=0}^{L1}p_i=_{i=0}^{L1}(1p_i)`$ for a fair game. In this case, the stationary solution can be written as the exponential of the potential $`P_i^{st}=N\mathrm{e}^{V_i/D}`$. Note that Eq.(12) reduces in the limit $`\mathrm{\Delta }x0`$ to $`V(x)=M^1F(x)𝑑x`$ or $`F(x)=M\frac{V(x)}{x}`$, which is the usual relation between the drift $`F(x)`$ and the potential $`V(x)`$ with a mobility coefficient $`M`$.

The inverse problem of obtaining the game probabilities in terms of the potential requires solving Eq. (12) with the boundary condition $`F_0=F_L`$222The singularity appearing for a fair game $`V_L=V_0`$ in the case of an even number $`L`$ might be related to the lack of ergodicity explicitely shown in hkk02  for $`L=4`$:

$$F_i=(1)^i\mathrm{e}^{V_i/D}\left[\frac{_{j=1}^L(1)^j[\mathrm{e}^{V_j/D}\mathrm{e}^{V_{j1}/D}]}{(1)^L\mathrm{e}^{(V_0V_L)/D}1}+\underset{j=1}{\overset{i}{}}(1)^j[\mathrm{e}^{V_j/D}\mathrm{e}^{V_{j1}/D}]\right]$$
(13)
These results allow us to obtain the stochastic potential $`V_i`$ (and hence the current $`J`$) for a given set of probabilities $`\{p_0,\mathrm{},p_{L1}\}`$, using (12); as well as the inverse: obtain the probabilities of the games given a stochastic potential, using (13). Note that the game resulting from the alternation, with probability $`\gamma `$, of a game A with $`p_i=1/2,i`$ and a game $`B`$ defined by the set $`\{p_0,\mathrm{},p_{L1}\}`$ has a set of probabilities $`\{p_0^{},\mathrm{},p_{L1}^{}\}`$ with $`p_i^{}=(1\gamma )\frac{1}{2}+\gamma p_i`$. For the $`F_i`$’s variables, this relation yields:

$$F_i^{}=\gamma F_i,$$
(14)
and the related potential $`V^{}`$ follows from (12).

We give now two examples of the application of the above formalism. In the first one we compute the stochastic potentials of the fair game $`B`$ and the winning game $`B^{}`$, the random combination with probability $`\gamma =1/2`$ of game $`B`$ and a game $`A`$ with constant probabilities, in the original version of the paradoxha99 . The resulting potentials are shown in figure 1. Notice how the potential of the combined game clearly displays the asymmetry under space translation that gives rise to the winning game.

The second application considers as input the potential

$$V(x)=A\left[\mathrm{sin}\left(\frac{2\pi x}{L}\right)+\frac{1}{4}\mathrm{sin}\left(\frac{4\pi x}{L}\right)\right]$$
(15)
which has been widely used as a prototype for ratchetsr02 . Using (13) we obtain a set of probabilities $`\{p_0,\mathrm{},p_{L1}\}`$ by discretizing this potential with $`\mathrm{\Delta }x=1`$, i.e. setting $`V_i=V(i)`$. Since the potential $`V(x)`$ is periodic, the resulting game $`B`$ defined by these probabilities is a fair one and the current $`J`$ is zero. Game $`A`$, as always is defined by $`p_i=p=1/2,i`$. We plot in figure 2 the potentials for game $`B`$ and for the game $`B^{}`$, the random combination with probability $`\gamma =1/2`$ of games $`A`$ and $`B`$. Note again that the potential $`V_i^{}`$ is tilted as corresponding to a winning game $`B^{}`$. As shown in figure 3, the current $`J`$ depends on the probability $`\gamma `$ for the alternation of games $`A`$ and $`B`$.

In summary, we have written the master equation describing the Parrondo’s games as a consistent discretization of the Fokker–Planck equation for an overdamped Brownian particle. In this way we can relate the probabilities of the games $`\{p_0,\mathrm{},p_{L1}\}`$ to the dynamical potential $`V(x)`$. Our approach yields a periodic potential for a fair game and a tilted potential for a winning game. The resulting expressions, in the limit $`\mathrm{\Delta }x0`$ could be used to obtain the effective potential for a flashing ratchet as well as its current.

The work is supported by MCyT (Spain) and FEDER, projects BFM2001-0341-C02-01, BMF2000-1108. P.A. is supported by a grant from the Govern Balear.


text:
# Photoelectric heating and [CII] cooling in translucent clouds: results for cloud models based on simulations of compressible MHD turbulence

## 1 Introduction

In diffuse interstellar clouds ($`A_\mathrm{V}1`$) the gas is heated mainly by photoelectrons expelled from dust grains by far–ultraviolet (FUV) photons (de Jong, 1977; Draine, 1978; Bakes & Tielens, 1994; Weingartner & Draine, 2001b) in the energy range 6 eV$`<h\nu <`$13.6 eV. The high energy limit corresponds to the cut–off in the FUV radiation field caused by the hydrogen absorption ($`h\nu `$=13.6 eV), while the low energy limit corresponds to the energy needed to free electrons from the grains ($`h\nu `$6 eV). In the cold neutral medium ($`T_{\mathrm{kin}}200`$ K) photoelectric heating accounts for most of the heating, the X-ray and cosmic ray heating rates being more than an order of magnitude smaller (Wolfire et al., 1995). In a relatively dense neutral medium ($`n_\mathrm{H}`$100 cm-3), where a significant fraction of carbon is in the neutral form, carbon ionization becomes an important heating source, but it is still not comparable to the photoelectric effect.

The efficiency of the photoelectric heating, $`ϵ`$, is defined as the ratio between the kinetic energy of photoelectrons available to heat the gas and the FUV energy absorbed by the grains. Bakes & Tielens (1994) estimated the photoelectric yield to be close to 10%, meaning that of all the photons absorbed by grains only one in ten ionizes them. Furthermore, since the average absorbed photon has an energy of $``$8 eV and the work function for a neutral grain is $``$5.5 eV, only one third of the absorbed energy goes into the kinetic energy of the electron. The overall efficiency is therefore expected to be close to 3%, at least in the cold neutral medium. The rate decreases as the grains become positively charged. In the intense FUV radiation field of photodissociation regions, for example, the efficiency can be more than a factor of ten lower than in the cold neutral medium.

Ingalls et al. (2002) have recently studied a sample of translucent, high latitude clouds (HLCs). The \[C+\] 158 $`\mu `$m line $`{}_{}{}^{2}P_{3/2}^{}^2P_{1/2}`$ was observed with ISO and this was correlated with FIR intensities. In the cold neutral medium the \[CII\] line dominates the cooling (Wolfire et al., 1995), neutral carbon being the next most important coolant. Ingalls et al. (2002) estimated that for these clouds the CI cooling rate is at most 30% of the \[CII\] cooling rate and the CO rate is more than one order of magnitude smaller. Therefore, for the observed HLCs the \[C+\] 158$`\mu `$m line should account for at least 60% of the cooling and in most lines of sight the percentage should be considerably higher.

To be a tracer of photoelectric heating the \[CII\] emission must originate from the same regions where the photoelectric effect is the main heating source and, conversely, all regions where the photoelectric heating is dominant must be cooled predominantly by the \[CII\] emission. As discussed by Ingalls et al. (2002), both conditions are approximately satisfied. The ionization potential of carbon, 11.3 eV, is above the work function of grains, and therefore C+ exists only in regions where photoelectric heating is possible. Conversely, in regions where C+ is not present and C and CO are significant coolants, the photoelectric heating rate must be much lower than in regions where C+ is present. As in Ingalls et al. (2002) we therefore assume that the above conditions are satisfied and the photoelectric heating is exactly balanced by the \[CII\] line cooling.

Since all the FUV radiation in the energy range 6 eV$`<h\nu <`$13.6 eV absorbed by grains contributes to the photoelectric heating, the efficiency $`ϵ`$ of this heating mechanism can be computed as the ratio of the observed \[CII\] line emission and the absorbed FUV radiation. If dust properties are assumed to be known, the absorbed FUV radiation can be estimated through the observed FIR emission. The problem of estimating $`ϵ`$ is therefore reduced to the problem of computing the correlation between the absorbed FUV radiation and the observable FIR emission. Such correlation depends on both the intensity of the radiation field and the column density and can be determined with model calculations. The degeneracy between the effect of the column density and of the intensity of the the external radiation field could in principle be broken using the information contained in the shape of the spectral energy distribution. This requires, however, quite precise modeling of the dust emission.

Ingalls et al. (2002) modeled their observations using a plane parallel geometry. The FIR emission was calculated with a dust model consisting only of grains at equilibrium temperature. A comparison with \[CII\] and FIR observations resulted in an estimate of $`ϵ=0.043`$ for the efficiency of the photoelectric heating. The average radiation field, $`\chi `$, was estimated to be above the local interstellar radiation field (ISRF), $`\chi 1.6\chi _0`$, where $`\chi _0`$ refers to the sum of the 2.7 K cosmic background and the spectrum given by Mezger et al. (1982) and Mathis et al. (1983).

In the present work we follow the same approach as in Ingalls et al. (2002), in the sense of deriving theoretical relations between the FUV absorption and the FIR emission and then using the empirical relation between \[CII\] and FIR intensities to estimate the average efficiency of the photoelectric heating. However, our models of the FUV absorption and FIR emission is significantly more detailed and realistic than in Ingalls et al. (2002). We compute the radiative transfer on cloud models based on three dimensional numerical simulations of compressible magneto–hydrodynamic turbulence. The density distribution of these turbulent flows provides a realistic model for the density inhomogeneity of interstellar clouds (Padoan et al., 1997a, b, 1998, 1999; Padoan & Nordlund, 1999). The penetration of FUV radiation into the cloud depends on the cloud structure. Inside an inhomogeneous cloud the intensity of short wavelength radiation is higher than inside a homogenous cloud. Furthermore, density variations generate a significant scatter in the relation between the local FUV absorption and FIR emission. Our model can therefore be used to estimate what fraction of the observed scatter can be attributed to the inhomogeneous nature of the density field. Other factors that could also contribute to the observed scatter, such as anisotropy in the radiation field or abundance variations, are not considered. Other improvements of our work, compared with Ingalls et al. (2002), are related to the employed dust model. We use the three component model of Li & Draine (2001). The FIR intensities are calculated using the method of Juvela et al. (2003) and the emission from transiently heated dust grains is included.

## 2 The cloud models

The density distributions of the models are the result of three dimensional simulations of super–Alfvénic, compressible, magneto–hydrodynamic (MHD) turbulence. Three models are used, from three simulations with different values of the rms sonic Mach number of the flow, $`M_\mathrm{S}=0.6`$, 2.5 and 10.0 in model $`A`$, $`B`$, and $`C`$, respectively.

The simulations are carried out on a staggered grid of 2503 computational cells, with periodic boundary conditions. Turbulence is set up as an initial large scale random and solenoidal velocity field (generated in Fourier space with power only in the range of wavenumbers $`1K2`$) and maintained with an external large scale random and solenoidal force, correlated at the largest scale turn–over time. The initial density and magnetic field are uniform and the gas is assumed to be isothermal. Details about the numerical method are given in Padoan & Nordlund (1999).

Experiments are run for approximately 10 dynamical times in order to achieve a statistically relaxed state. The cloud models used in this work correspond to the final snapshot of each simulation. The value of $`M_\mathrm{S}`$ is varied in different experiments by varying the thermal energy. The initial rms Alfvénic Mach number remains unchanged, $`M_\mathrm{A}=10.0`$, from run to run. The volume–averaged magnetic field strength is constant in time because of the imposed flux conservation. The magnetic energy is instead amplified. The initial value of the ratio of average magnetic and dynamic pressures is $`P_\mathrm{m}_{\mathrm{in}}/P_\mathrm{d}_{\mathrm{in}}=0.005`$ for all runs, so all the runs are initially super–Alfvénic. The value of the same ratio at later times is larger, due to the magnetic energy amplification, but still significantly lower than unity (0.21, 0.14 and 0.12 for models $`A`$, $`B`$ and $`C`$ respectively). The turbulence is therefore super–Alfvénic at all times.

Supersonic and super–Alfvénic turbulence of an isothermal gas generates a highly inhomogeneous gas density distribution, with a density contrast of several orders of magnitude. It has been shown to provide a good description of the dynamics of molecular clouds and of their highly fragmented nature (Padoan et al., 1997a, b, 1998, 1999; Padoan & Nordlund, 1999). Translucent high latitude clouds (HLCs), such as the ones observed by Ingalls et al. (2002), have been studied thanks to their CO emission (Magnani et al., 1996). Despite their lower gas density and presumably larger fraction of gas in the atomic form, the clouds can still be approximately modeled with an isothermal equation of state. At densities above 1 cm-3 the gas temperature is expected to remain mostly in the range 50–100 K (e.g. Wolfire et al., 1995). The observed clouds are characterized by supersonic turbulent motions like denser molecular clouds (see e.g. Miesch et al., 1999). It seems therefore appropriate to use supersonic turbulence as a model not only of dense molecular clouds, but also of more diffuse clouds, such as the ones studied in this work. A subsonic model (model $`A`$) is also used in this work for the purpose of studying the effect of varying $`M_\mathrm{S}`$ over a large range of values.

To compute the radiative transfer the data cubes of the density distributions are scaled to physical units by fixing the length of the computational grid, $`L=6.0`$ pc, and the average hydrogen density, $`n_\mathrm{H}=100`$ cm-3. However, the radiative transfer calculations apply to all values of $`L`$ and $`n_\mathrm{H}`$ satisfying the condition that their product is constant, $`Ln_\mathrm{H}=600`$ pc cm-3= $`1.8\times 10^{21}`$ cm-2. The column density in the models varies between different lines of sight, because of the inhomogeneous density field. Variations are particularly large in model $`C`$, where the range of column density values is comparable to that found in the Ingalls et al. (2002) sample. This can be seen for example by comparing the FIR intensities in the models and in the observations (see Sect. 4). When results are compared with observations, the angular size of the models (or their distance) must be fixed as well, in order to take into account the spatial resolution of the observations. One cell, 6.0/250 pc, is assumed to correspond to 0.25$`\mathrm{}`$, setting the model clouds at the distance of 330 pc.

## 3 Calculation of FUV absorption and FIR emission

Radiative transfer calculations are needed to obtain estimates for the \[CII\] emission and the total FIR emission along different lines of sight through the model clouds. Calculations are carried out on the original grid of 2503 cells. The external radiation field is assumed to be the local interstellar radiation field near the Sun (Mathis et al., 1983). The dust model is taken from Li & Draine (2001). This includes silicate grains (sizes $`a>3.5`$ Å), graphite grains ($`a50`$ Å) and PAHs (from $`a=3.5`$ Å to $`50`$ Å) with ionization state corresponding to physical conditions in the cold neutral medium (see also Juvela et al., 2003).

Following Ingalls et al. (2002) we assume explicitly a balance between the cooling by \[CII\] line emission and the photoelectric heating. The latter is taken to be equal to the dust absorption in the energy range 6–13.6 eV multiplied by the efficiency of the photoelectric emission, $`ϵ`$. The penetration of the external radiation field is simulated with Monte Carlo methods, and the absorption in the given energy range is stored in each cell. The emerging \[CII\] intensity is calculated by summing the corresponding emission on selected lines–of–sight, assuming the emission is optically thin. The efficiency of the photoelectric heating is unknown but since it is assumed to be constant troughout the cloud it can be taken into account later in the analysis.

The dust emission was calculated using the ‘library method’ discussed in Juvela et al. (2003). First a mapping between the strength of the local radiation field at a few reference frequencies and the resulting dust emission is established. The radiation field at the reference frequencies is then simulated and the emission is obtained using the existing mapping. The emission, including the contribution from transiently heated particles, can be estimated with an accuracy of a few per cent even in the case of large three–dimensional models (see Juvela et al., 2003). The infrared emission from the clouds is estimated using this method. Ingalls et al. (2002) calculated FIR emission based on 60$`\mu `$m and 100$`\mu `$m values using the formula

$$I_{\mathrm{FIR}}=12.6\frac{I(100\mu m)2.58I(60\mu m)}{\mathrm{MJy}\mathrm{sr}^1}\times 10^6\frac{\mathrm{erg}}{\mathrm{cm}^2\mathrm{s}\mathrm{sr}}.$$
(1)
The formula gives the FIR emission between 42.5 $`\mu `$m and 122.5 $`\mu `$m i.e. the wavelength range covered by the IRAS 60$`\mu `$m and 100$`\mu `$m filters (Helou et al., 1985). In our models the formula was found to be accurate to $`23`$ per cent. However, the quoted FIR intensities are here calculated by direct integration over the computed spectra.

Both the FUV absorption (equivalent to the \[CII\] emission divided by the efficiency, $`I_{\mathrm{CII}}/ϵ`$) and the FIR intensity, $`I_{\mathrm{FIR}}`$, are calculated towards three directions perpendicular to the faces of the cubic cloud. For each direction, maps of 250$`\times `$250 positions are generated. Results are compared with observations presented by Ingalls et al. (2002). The ISO LWS observations of the \[C+\] 158$`\mu `$m line have a resolution of 71$`\mathrm{}`$. The FIR data are based on IRAS ISSA maps with a resolution of $`45\mathrm{}`$. We assume the cells of our models have an angular size of 0.25$`\mathrm{}`$, and the computed \[CII\] and FIR maps are convolved to the resolution of the observations.

## 4 Results

### 4.1 FIR and \[CII\] emission in model clouds

For each model the FUV absorption and the FIR intensity, $`I_{\mathrm{FIR}}`$ are calculated towards three directions perpendicular to the faces of the density field data cubes. Since the FUV absorption is assumed to be equivalent to the \[CII\] line intensity, $`I_{\mathrm{CII}}`$, divided by the efficiency, $`ϵ`$, in the following we will refer to the quantity $`I_{\mathrm{CII}}/ϵ`$ (with $`ϵ`$ still unknown), instead of the computed FUV absorption.

Fig. 1 shows the distribution of the ratio $`I_{\mathrm{CII}}/(ϵI_{\mathrm{FIR}})`$ in the three models. The figure is obtained at the original resolution (250$`\times `$250 pixels), and neither $`I_{\mathrm{FIR}}`$ nor $`I_{\mathrm{CII}}/ϵ`$ are spatially convolved. Model $`A`$ is based on a flow with subsonic turbulence, with only small density fluctuations. This is reflected in the top panel of Fig. 1 showing a smooth distribution of $`I_{\mathrm{CII}}/(ϵI_{\mathrm{FIR}})`$ in model $`A`$. The density contrast increases with increasing $`M_\mathrm{S}`$. Model $`B`$, with $`M_\mathrm{S}=2.5`$, shows a more ‘clumpy’ distribution of $`I_{\mathrm{CII}}/(ϵI_{\mathrm{FIR}})`$ (middle panel of Fig. 1). Model $`C`$, with $`M_\mathrm{S}=10`$, has clearly the most inhomogeneous and filamentary distribution (bottom panel of Fig. 1).

In the case of the nearly homogeneous cloud $`A`$, both $`I_{\mathrm{FIR}}`$ and $`I_{\mathrm{CII}}/ϵ`$ are determined mainly by the distance to the cloud surface. Since $`I_{\mathrm{CII}}/ϵ`$ is assumed to be equivalent to the FUV absorption, $`I_{\mathrm{CII}}/ϵ`$ depends only on photons in the energy range 6–13.6 eV, while $`I_{\mathrm{FIR}}`$ depends on a much wider wavelength range responsible for heating the dust grains. The dust extinction increases towards shorter wavelengths, and this leads to a reduction in the ratio $`I_{\mathrm{CII}}/(ϵI_{\mathrm{FIR}})`$ inside the clouds. At the cloud surface the ratio $`I_{\mathrm{CII}}/(ϵI_{\mathrm{FIR}})`$ is above one, and it decreases to $``$0.8 inside the cloud, in all three models. The total range of values is wider in the more inhomogeneous clouds, with the lowest values reached in the highest density peaks. The average ratios are, however, very similar in all three models, $`I_{\mathrm{CII}}/(ϵI_{\mathrm{FIR}})`$0.9. The average value depends, of course, on the total column density and it would become lower in more opaque clouds. The models have similar column densities as the clouds in the sample of Ingalls et al. (2002), and if the clouds are subjected to similar radiation field as the models the exact correspondence can be checked by comparing the FIR intensities.

In Fig. 2 we plot the ratio $`I_{\mathrm{CII}}/(ϵI_{\mathrm{FIR}})`$ against the visual extinction of the corresponding line of sight. The average optical depth of the model clouds is below one, but the total range of $`A_\mathrm{V}`$ values depends on the degree of inhomogeneity. While in model $`A`$ the visual extinction stays below $``$1m, the other models include a number of lines of sight with $`A_\mathrm{V}`$ in excess of 2m (model $`B`$) and 3m (model $`C`$). The figure shows the results for both unconvolved data and for the case where data have been convolved to the resolution of the observations (71$`\mathrm{}`$ for \[CII\] and $``$4$`\mathrm{}`$ for FIR data). The $`A_\mathrm{V}`$ data has not been convolved. This would correspond to a situation where the extinction is derived from photometry of an individual background star, probing only a single line of sight inside the larger beam employed in the other observations.

The negative correlation between extinction and the ratio $`I_{\mathrm{CII}}/(ϵI_{\mathrm{FIR}})`$ is clearly visible in model $`A`$ (Fig. 2a). In model $`B`$ the correlation can be seen in the unconvolved data, while in the convolved data the large difference between the resolution of the \[CII\] and FIR observations increases the scatter and the correlation becomes weak. At high extinction values the convolution increases systematically the $`I_{\mathrm{CII}}/(ϵI_{\mathrm{FIR}})`$ values. The \[CII\] intensity depends on FUV absorption, and because of the high optical depth the \[CII\] emission is likely to saturate near the column density maxima. Since the intensity distribution is already flat the convolution has little effect on the peak \[CII\] intensities. The FIR emission is less saturated and, more importantly, the FIR intensities are convolved with a beam that is several times larger. As the the FIR peaks are smoothed away the ratio $`I_{\mathrm{CII}}/(ϵI_{\mathrm{FIR}})`$ increases. The effect is even higher in model $`C`$ because of the stronger density contrast. However, in model $`C`$ the $`A_\mathrm{V}`$–dependence is quite weak even in the unconvolved data. This can be attributed to other effects related to the increased density contrast. For example, dense filaments and sheets are common, with the result that the $`A_\mathrm{V}`$ along one line of sight can be much higher than the average extinction by which the material is shielded towards other directions. In the figure this results in a larger scatter along the $`A_\mathrm{V}`$ axis.

In a homogeneous model the extinction of the external radiation field proceeds in a predictable fashion, and the ratio between optical and FUV intensities only depends on the distance from the cloud surface. In an inhomogeneous cloud there is no longer a one–to–one relation between the shape of the spectral energy distribution and the intensity of the radiation. The FUV radiation penetrates deeper into the cloud than in the homogeneous model, and the variation in the ratio between optical and FUV intensities becomes smaller. If the cloud consisted of optically very thick clumps, the incoming radiation would have exactly the same spectral energy distribution at the centre of the cloud as on the surface – only the intensity would be smaller. MHD models are somewhere between the two extremes, that is between a homogeneous cloud and a cloud with optically thick clumps. In different parts of the cloud the spectrum of the radiation field may be different even for the same ‘effective extinction’ and this increases the scatter along the $`I_{\mathrm{CII}}/(ϵI_{\mathrm{FIR}})`$–axis.

In Fig. 3 the dependence on the visual extinction along the line of sight is shown separately for $`I_{\mathrm{CII}}/ϵ`$ and $`I_{\mathrm{FIR}}`$. This shows explicitly the saturation of the \[CII\] emission, which at high column densities is caused by the total absorption of incoming FUV intensity. The \[CII\] line itself is assumed to be optically thin. The effect is particularly clear in model $`B`$. It is present also in model $`C`$, but it is partially masked by the scatter resulting from the more inhomogeneous cloud structure and the effects of the convolution.

Inhomogeneity is expected to lead to a decrease in total energy absorption (see Juvela et al., 2003) and, therefore, to lower the emitted intensities. The FIR emission in models $`B`$ and $`C`$ is respectively 2% and 5% lower than in model $`A`$. The FUV absorption, or the equivalent quantity $`I_{\mathrm{CII}}/ϵ`$, is 3.6% and 9.0% lower in models $`B`$ and $`C`$ respectively than in model $`A`$. The larger variation of $`I_{\mathrm{CII}}/ϵ`$ relative to $`I_{\mathrm{FIR}}`$ is the result of a difference in the optical depths relevant to the two components. For optically thin radiation the absorption depends only on the total amount of absorbing material while for optically thick radiation the absorbed energy is proportional to the area filling factor. Since the optical depth is higher for the FUV radiation than for the longer wavelengths heating the grains, the quantity $`I_{\mathrm{CII}}/ϵ`$ is more sensitive to the density contrast (filling factor) of the cloud. The equilibrium dust temperature for large grains, and consequently the FIR emission, depends on longer wavelengths that are partly optically thin. Therefore, in more clumpy clouds the reduction in $`I_{\mathrm{FIR}}`$ is more modest than in $`I_{\mathrm{CII}}/ϵ`$. However, the value $`I_{\mathrm{CII}}/(ϵI_{\mathrm{FIR}})`$ decreases from model $`A`$ to model $`C`$ by less than 2%. The interpretation of observations is affected by several other uncertainties, such as the imprecise knowledge of the dust properties or the strength of the external radiation field (see Sect. 5). A variation of 2% is insignificant compared with these other uncertainties.

### 4.2 Efficiency of the photoelectric heating

In Fig. 4 the quantity $`I_{\mathrm{CII}}/ϵ`$ is plotted versus the FIR intensity, $`I_{\mathrm{FIR}}`$, for the models $`A`$, $`B`$ and $`C`$ (left to right) observed along the $`y`$–direction. We show the relations for both the original data and the data convolved to the resolution of the observations discussed in Ingalls et al. (2002), that is $`1.2\mathrm{}`$ for \[CII\] and $`4\mathrm{}`$ for the FIR observations. The relations are close to linear, and Fig. 4 shows the least squares lines fitted to the convolved data. The range of intensities and the scatter in the relations both increase from model $`A`$ to model $`C`$. The slope of the relation decreases, albeit very slowly with increasing Mach number of the model. The fitted slopes for each of the three orthogonal directions of the lines of sight are given in Table 1.

In Sect 4.1 the \[CII\] intensities were found to saturate at high column densities when all incoming FUV radiation was absorbed. As a result, at high intensities the points for unconvolved data fall below the least squares line. However, at the location of intensity maxima the convolution decreases the FIR intensities more than the \[CII\] intensities, and in the figure the convolved points are shifted to the left more than down. This is seen especially in model $`C`$ that has the largest density contrast. In model $`C`$ the relation $`I_{\mathrm{CII}}/ϵ`$ versus $`I_{\mathrm{FIR}}`$ for the convolved data turns slightly upwards at higher intensities, relative to the relation for the unconvolved data. In model $`B`$ the convolution raises the points at high-$`A_\mathrm{V}`$ only slightly closer to the fitted least squares line.

A comparison between the ratios $`I_{\mathrm{CII}}/(ϵI_{\mathrm{FIR}})`$ obtained from the models and the empirical relation $`I_{\mathrm{CII}}/I_{\mathrm{FIR}}`$ allows to estimate the efficiency of the photoelectric heating, $`ϵ`$. In Ingalls et al. (2002) the average ratio of CII and FIR intensities is found to be $`I_{\mathrm{CII}}/I_{\mathrm{FIR}}`$=(2.5$`\pm `$0.9)$`\times 10^2`$, for a sample of high latitude clouds observed with ISO. The slope they derive for large scale correlation using DIRBE and FIRAS observations for $`|b|>5\mathrm{deg}`$ is also very similar, $`I_{\mathrm{CII}}`$=(2.54$`\pm `$0.03)$`\times 10^2I_{\mathrm{FIR}}`$ \- (1.1365$`\pm `$0.0004)$`\times 10^6`$ erg cm-2 s-1 sr-1.

We have calculated the efficiencies of the photoelectric emission, $`ϵ`$, assuming $`I_{\mathrm{CII}}/I_{\mathrm{FIR}}=2.5\times 10^2`$. The results are given in the fourth column of Table 1. We obtain an average value of $`ϵ=2.88\times 10^2`$. There is very little variation between the models or the direction from which the model clouds are observed. However, as seen in Fig. 1, the variation within individual maps can be quite significant. In the case of unconvolved data in Fig. 1 the variation ranges from $``$40% in model $`A`$ to almost 90% in model $`C`$.

## 5 Discussion

The observations presented by Ingalls et al. (2002) contain scatter in the plot of $`I_{\mathrm{CII}}`$ versus $`I_{\mathrm{FIR}}`$ exceeding the uncertainty of individual observations. Differences between clouds can be attributed largely to differences in the dust composition. The photoelectric heating and therefore \[CII\] emission depends on the abundance of very small dust grains and PAH–molecules while the far–infrared emission is caused mainly by large dust grains. The ratio between the abundances of very small grains and large grains is known to show significant spatial variations. The variations are more noticeable in dense clouds (e.g. Laureijs et al., 1991; Juvela et al., 2002; Stepnik et al., 2003), but changes have been observed even in cirrus type clouds (Bernard et al., 1999; Cambrésy et al., 2001). Variations in the radiation field are probably less important.

In the first approximation an increase in the strength of the radiation field affects FIR and \[CII\] intensities in the same way. However, as the increased dust temperature moves dust emission to shorter wavelengths the FIR intensity is reduced relative to the \[CII\] emission (FUV absorption). Assuming that the FIR intensity in the models is fixed and constrained by the observations, if the interstellar radiation field is increased the dust column density must be decreased to maintain the FIR intensity unchanged. As a result of the decreased column density, the \[CII\] emission is less saturated, and the ratio $`I_{\mathrm{CII}}/I_{\mathrm{FIR}}`$ is increased. Variations within individual clouds could similarly be caused by variations in the local radiation field or in the properties of the gas or dust.

In the present models dust properties are assumed to be constant throughout the clouds. The fractional abundance of C+ is also kept constant and the \[CII\] emission is assumed to be responsible for all of the gas cooling. Even with these simplifying assumptions the computed ($`I_{\mathrm{CII}}`$,$`I_{\mathrm{FIR}}`$) values contain a significant scatter. The scatter is smallest in the most homogeneous model $`A`$. Only in a completely homogeneous model the points would fall on one line as both \[CII\] and FIR emission would depend solely on the distance from the cloud surface. The scatter is largest in model $`C`$, with the largest value of $`M_\mathrm{S}`$ and density contrast. In this model, both the range of $`I_{\mathrm{CII}}`$ and $`I_{\mathrm{FIR}}`$ values and the scatter in their relation are similar to the observations presented by Ingalls et al. Such a comparison is relevant only if the observed positions were selected randomly within the clouds. This is not always the case, but the scatter of the observed points should still be close to the true variation within the observed clouds.

In Fig. 5 we compare our model $`C`$ with the observation presented by Ingalls et al. (2002). We have plotted in the same figure our model $`C`$ (direction $`y`$) for which the \[CII\] intensities have been calculated using the value $`ϵ=2.9\times 10^2`$ derived in Sect. 4.2. The model covers a large fraction of the observed range of $`I_{\mathrm{CII}}`$ and $`I_{\mathrm{FIR}}`$ values. After subtracting a linear fit, the remaining scatter in all observations is about twice the scatter in model $`C`$, even when the observational errors are taken into account. On the other hand, in individual clouds the variation is very similar to that found in the model. Our results show that most of the observed scatter may originate from the inhomogeneous density structure of the clouds, and the observed scatter is consistent with that predicted with supersonic and super–Alfvénic MHD turbulence with sonic rms Mach number $`M_\mathrm{S}2`$. In subsonic models both the scatter and the total range of \[CII\] and FIR intensities is instead significantly smaller than in the observed clouds.

Compared with model $`C`$ the scatter is significantly larger only in the cloud MBM-12, which also has the highest FIR intensities. MBM-12 has an HI column density of $``$1.1$`\times 10^{21}`$ cm-2, i.e. some 40% below the average column density of the models. Timmermann et al. (1998) suggested the cloud is subjected to stronger than average radiation field, $`\chi `$1–2 $`\chi _0`$, and this could also explain the high FIR intensities. On the other hand, since the cloud has significant molecular line emission (Magnani et al., 1996) the total column density may not be very different from the model, and therefore the observed FIR intensities do not imply a strong radiation field by themselves. Furthermore, if one would accept the hypothesis of a strong radiation field, the low FIR colour temperature (see Ingalls et al., 2002, Fig. 4) would require a very low abundance of small dust grains.

The scatter in the observed \[CII\] and FIR intensities may depend on two effects not included in the models: i) CI cooling and ii) optical depth of the \[CII\] line. Although the models assume all the cooling is from C+, CI emission may be important on some lines of sight, increasing the observed scatter. According to Ingalls et al. (2002) the maximum contribution to the cooling from CI is approximately 30%. The scatter due to the the presence of CI cooling on some lines of sight should therefore be well below 15%. Furthermore, the CI contribution can be significant only for high column density, and could not explain the scatter at low FIR intensities. The models also assume that the \[CII\] line is optically thin. If the line became optically thick at some lines of sight, the main effect would be to flatten the $`I_{\mathrm{CII}}`$ vs. $`I_{\mathrm{FIR}}`$ relation at high column densities, while the scatter would not be significantly affected. In the models, the saturation of $`I_{\mathrm{CII}}`$ is caused by the total absorption of the FUV radiation that is converted into \[CII\] emission, not by self–absorption of the \[CII\] line itself. Since the observations do not show strong \[CII\] saturation with increasing FIR intensity, the optical depth of the \[CII\] line cannot be very high.

The derived value of the photoelectric heating efficiency, $`ϵ`$, is close to theoretically derived values (de Jong, 1977; Bakes & Tielens, 1994). Bakes & Tielens (1994) performed theoretical calculations using a MRN type grain size distributions ($`na^{3.5}`$, Mathis et al., 1977). They found that about half of the heating rate is due to grains with sizes below 15 Å and the contribution from grains larger than 100 Å is insignificant. The net efficiency was calculated as a function of the parameter $`G_0(T/\mathrm{K})^{1/2}(n_\mathrm{e}/\mathrm{cm}^3)^1`$, where $`G_0`$ is the intensity of the far-UV radiation field relative to Habing’s field (Habing, 1968), $`T`$ is gas temperature and $`n_\mathrm{e}`$ the electron density. The estimated efficiency is constant $`ϵ3\times 10^2`$ for $`G_0(T/\mathrm{K})^{1/2}(n_\mathrm{e}/\mathrm{cm}^3)^110^3`$. This prediction should apply to the cloud sample in Ingalls et al. (2002). Bakes & Tielens found, however, that the predicted heating rate could be doubled if the exponent of the grain size distribution was decreased from -3.5 to -4.0. Since $`A_\mathrm{V}`$ and FIR intensity depend on large grains and the heating mainly on PAHs the derived efficiency $`ϵ`$ is model dependent.

In this paper we have used the dust model of Li & Draine (2001). Weingartner & Draine (2001b) based their studies on the same dust model and considered further different size distributions discussed in Weingartner & Draine (2001a). For cold neutral medium ($`T`$100 K) and size distributions consistent with the diffuse cloud extinction law, $`R_\mathrm{V}`$3.1, their predictions for the photoelectric heating rate are close to the values given by Bakes & Tielens (1994). The results depend on the carbon abundance which is limited by extinction measurements to values below $`6\times 10^5`$ carbon atoms per H nucleus (Weingartner & Draine, 2001a, b; Li & Draine, 2001). The models of Li & Draine (2001) and Weingartner & Draine (2001a) favoured carbon abundances close to this upper limit and the resulting rate of the photoelectric heating is some 25% above the value given by Bakes & Tielens (Weingartner & Draine, 2001b).

Habart et al. (2001) derived values of $`ϵ`$=2–3% for the efficiency of the photoelectric heating across the cloud L1721. Compared with the Ingalls et al. sample, the cloud has higher visual extinction, $`A_\mathrm{V}3^m`$, and, due to the proximity of a B2 star, it is subjected to a stronger radiation field, $`\chi 5\chi _0`$. By adopting the dust model of Désert (1990), Habart et al. were able to derive efficiencies of the photoelectric emission separately for PAH, very small grain and large grain components. The derived efficiencies were $`3`$% for PAHs, $``$1% for very small grains and $``$0.1% for large grains. Part of the observed variation was attributed to changes in the abundance of the dust components. The clouds in the Ingalls sample have lower visual extinction, and abundance variations are expected to be correspondingly smaller.

Another source of uncertainty is the balance between the column density and the strength of the radiation field. As discussed above, if a stronger radiation field or a smaller column density were adopted in the models, the calculated ratio $`I_{\mathrm{CII}}/(ϵI_{\mathrm{FIR}})`$ would increase and the derived value of $`ϵ`$ decrease. According to Fig. 1, the ratio $`I_{\mathrm{CII}}/(ϵI_{\mathrm{FIR}})`$ increases up to 30% above its average value when moving from the cloud center to the cloud edges, where the radiation field is not attenuated. Therefore, if the clouds had significantly lower extinction, the efficiency of the photoelectric heating could be lower by the same amount. We modified the model $`C`$ by reducing its column density by a factor of three and by multiplying the intensity of the external radiation field by a factor of two. The estimated average efficiency of the photoelectric heating decreased from 2.88$`\times 10^2`$ to 2.25$`\times 10^2`$, i.e. only by 22%. Based on HI and molecular line data the column densities are in most of the observed clouds within a factor of two from the average column density of the models. In this paper we have used for the intensity of the local ISRF values given by Mezger et al. (1982) and Mathis et al. (1983). According to Draine (1978) and Parravano et al. (2003) the actual ISRF could be up to $``$70% stronger. However, the difference is very unlikely to be a factor of two. The actual errors caused by the uncertainty of the column density and intensity values is therefore less than 22%.

Compared with Ingalls et al. (2002) our estimate of the photoelectric heating efficiency is smaller by one third. Most of the difference is due to a difference in the predicted FIR intensity. For example, for a model with optical depth $`A_\mathrm{V}`$=1.0 mag through the cloud Ingalls et al. obtained a 100 $`\mu `$m surface brightness of slightly more than 15 MJy sr-1. For a homogeneous, spherically symmetric cloud with equal optical depth we obtain a surface brightness of $``$9.5 MJy sr-1. This is close to what Bernard et al. (1992) obtained for a similar cloud model using the dust model of Désert (1990). The difference between our result and Ingalls et al. cannot be due to model geometry. One would expect more emission (higher $`I_{\mathrm{FIR}}`$) from a three dimensional cloud heated from all directions than from a plane parallel cloud heated on two surfaces only. The most likely explanation for the high FIR intensities obtained by Ingalls et al. (2002) is their assumption of thermal equilibrium. This means that all grains are constantly at temperatures close to 20 K and almost all absorbed energy is re-radiated in the FIR. In our case a larger fraction of the absorbed energy is radiated at shorter wavelengths because of the small particles that are temporarily heated to much higher temperatures.

## 6 Conclusions

We have calculated FIR emission and \[CII\] line emission for three–dimensional density distributions of compressible magneto–hydrodynamic turbulent flows with rms sonic Mach numbers $`M_\mathrm{S}=`$ 0.6, 2.5, and 10.0. The dust emission, $`I_{\mathrm{FIR}}`$, is computed with full radiative transfer calculations. The \[CII\] emission is estimated assuming the photoelectric heating caused by FUV photons between 0.0912$`\mu `$m and 0.2066$`\mu `$m is balanced by \[CII\] line emission. The FUV absorption is determined by the radiative transfer simulations, and is assumed to be equal to the \[CII\] line intensity divided by the unknown efficiency of the photoelectric heating, $`I_{\mathrm{CII}}/ϵ`$.

The average ratio $`I_{\mathrm{CII}}/(ϵI_{\mathrm{FIR}})`$ is in all models between 0.85 and 0.88, showing that its dependence on the density contrast (Mach number) is weak. However, dense filaments are visible in the maps as regions with lower value of $`I_{\mathrm{CII}}/(ϵI_{\mathrm{FIR}})`$. The degree of correlation between $`I_{\mathrm{CII}}/(ϵI_{\mathrm{FIR}})`$ and visual extinction decreases in more inhomogeneous clouds. In the case of simulated observations convolved to different resolutions (71$`\mathrm{}`$ for \[CII\] and $``$4$`\mathrm{}`$ for FIR) most of the correlation is lost.

The scatter in the observational ($`I_{\mathrm{FIR}}`$,$`I_{\mathrm{CII}}/ϵ`$) plot can be reproduced by models with rms Mach number $`M_\mathrm{S}`$2 (supersonic turbulence), showing that most of the scatter may be due to the inhomogeneous nature of the clouds (likely due to the turbulence). In subsonic models ($`M_\mathrm{S}`$1) both the scatter and the total range of FIR and \[CII\] intensities become smaller than in the observed clouds. Using the empirical value $`I_{\mathrm{CII}}/I_{\mathrm{FIR}}`$=2.5$`\times 10^2`$ found for high latitude clouds (Ingalls et al., 2002) the efficiency of the photoelectric heating is found to be $`ϵ2.9\times 10^2`$. The value is very close to the theoretical predictions for the cold neutral medium.

MJ acknowledges the support from the Academy of Finland Grants no. 174854 and 175068. The work of PP was partially performed while PP held a National Research Council Associateship Award at the Jet Propulsion Laboratory, California Institute of Technology. The work of RJ is supported in part by NSF grant AST-0206031.


text:
# Experimental multi-photon-resolving detector using a single avalanche photodiode

## 1 Introduction

There is an increasing need for identifying the number of photons that is contained within a weak light pulse or time interval of a weak cw light field. A device capable of photon-number resolution would contribute both to fundamental research in the area of quantum optics and to more-or-less practical quantum communication systems, such as quantum key distribution schemes . In the former case it would significantly help in the analysis and preparation of quantum states with prescribed statistics , in the latter it would prolong communication distances and increase the rate of secure communication . A detector capable to resolve the vacuum state, one-photon Fock state and states containing more than one photon was recently constructed , but it requires operation under extreme conditions at present and it hardly becomes a common laboratory tool at least in near future. A similar problem occurs with a photon-number resolving detector  based on superconducting transition-edge sensor microcalorimeter technology that needs mK temperatures for its operation.

Common detectors of weak light fields (avalanche photodiodes and photomultipliers) cannot provide photon-number resolution but they have high quantum efficiencies in the visible  and near infrared  regions. The use of a cascade of such detectors behind a $`1\times N`$ multiport  then provides a device capable to resolve photon-numbers to some extent. Assuming a lossless device, ideal detectors, and provided that the mean number of photons in the signal $`\mu N`$, a part of any multi-photon signal containing $`k`$ photons gets split with a high probability to the arms of the multiport in such a way that $`k`$ detections at different detectors occur. (We note that in the limit $`N\mathrm{}`$ this probability is one.) However, it was shown in  that the performance of such a device is severely inflicted by losses in the device and imperfections of the detectors. To achieve a reasonable performance, large array of detectors with high quantum efficiency and low noise is required , so that such photon-number resolving device becomes unacceptably complex. The use of this photon-number resolving device in a source of one-photon Fock states based on entangled photon pairs and postselection was analyzed in . It was shown that photon fields with the Fano factor around 0.7 can be generated by this source under real conditions. A quantum-key distribution system using this source can have a secure communication distance up to 120 km and can provide higher values of gain . These results motivate the endeavour to simplify the above discussed photon-number resolving device. We note that also single-photon sources utilizing NV centers  and quantum dots  are perspective for quantum-key distribution systems. However, NV centers have low efficiencies of single-photon generation ($`10^3`$) and quantum-dot sources have to be cooled to 5 K at present.

In this letter we propose and test a variant of such cascading device in which we replace the $`1\times N`$ multiport and $`N`$ detectors with a fiber-loop delay line and a single avalanche-photodiode detector. This decreases the complexity as well as the cost of the device to a reasonable measure. A similar device has been recently suggested also by other authors . The number of detectors $`N`$ is given by the number of time windows (channels) we detect using the time-of-flight spectrometer. We use a variable-ratio coupler at the entrance to the fiber-loop delay line so that we can control the distribution of probabilities of detection in the individual channels. In order to show the performance of this fiber-loop detection device, we theoretically determine and experimentally test the best setting of the device with respect to distinguishing multi-photon states (containing two or more photons) from the one-photon state and the vacuum state.

## 2 Description of the device

The scheme of our device is plotted in Fig. 1. Quantum state to be analyzed is injected to port 1 of the single-mode variable ratio coupler (SIFAM SVR-82, referred as SVR hereafter). Port 3 outputs to a detector module based on a silicon avalanche photodiode with active quenching (Perkin-Elmer SPCM-AQ-141-FC). This detector yields TTL output upon incident photon with quantum efficiency $`\eta `$ of 60% (at 830 nm, including coupling-optics losses) and is not capable of resolving the number of incident photons. The dark count rate is about 40 counts per second. Ports 2 and 4 are interconnected with 10 m long single-mode fiber thus forming a delay line longer than the dead time of the detector which is 50 ns. The outputs of the detector are registered as stop pulses in the time-of-flight spectrometer (Fast ComTec TOF7885) that is connected to a PC directly or through a multichannel buffer. Start pulses for the time-of-flight spectrometer are generated by trigger pulses of the source of quantum states to be examined.

Denoting the intensity transmission coefficients from port $`i`$ to port $`j`$ of the SVR by $`t_{ij}`$, the device input coupling transmission by $`t_0`$ (due to loss at the input connector), the SVR transmission by $`\theta `$ (due to SVR excess loss), the total transmission of the fiber loop by $`t_l`$ and the transmission from port 3 of the SVR to the detector including detector quantum efficiency by $`\eta `$, we arrive at the transmission coefficients $`h_1,h_2,\mathrm{}`$ for the detection channels:

$`h_1`$ $`=`$ $`t_0\theta t_{13}\eta ,`$ (1)
$`h_k`$ $`=`$ $`t_0t_{14}\theta ^kt_l^{k1}t_{23}t_{24}^{k2}\eta ,k2.`$ (2)
If a single photon enters the device, transmission coefficients give also the probabilities of photon detection.

A question arises, what is the optimum setting of the SVR for the detection of multi-photon states. The highest probability to distinguish multi-photon states from the one-photon state and the vacuum state is achieved provided that the probabilities $`h_1,\mathrm{},h_N`$ have the same value. The reason is that the energy in a detected photon field is equally distributed over all detectors and then the intensities of the photon fields in all detectors reach the lowest possible value. This case can be also distinguished by maximizing Shannon entropy. We cannot reach this best case with our fiber-loop detection device. Best conditions from the point of view of photon-number resolution obtainable by our fiber-loop detection device can be found only numerically in principle. However, numerical calculations as well as experimental results indicated that the principle of maximization of Shannon entropy is valid for this task and can provide these conditions (see Fig. 4).

We first simplify the real SVR to an idealized unitary device setting $`\theta =t_0=1`$, $`t_{13}+t_{14}=1`$, $`t_{23}+t_{24}=1`$, and $`t_{13}+t_{23}=1`$. This idealized coupler can then be represented with a single variable-division-ratio parameter $`r`$. Upon replacing $`t_{13}=t_{24}r`$ and $`t_{14}=t_{23}(1r)`$ we get:

$`h_1`$ $`=`$ $`\eta r,`$ (3)
$`h_k`$ $`=`$ $`\eta (1r)^2t_l^{k1}r^{k2},k2.`$ (4)
Assuming further $`t_l=\eta =1`$ in the ideal case, we can evaluate Shannon entropy

$$E=\underset{i}{}h_i\mathrm{ln}(h_i)$$
(5)
of the ideal channel detector as follows:

$$E=2r\mathrm{ln}(r)2(1r)\mathrm{ln}(1r).$$
(6)
This entropy is maximized for $`r=1/2`$, i.e., a balanced SVR is optimal. In the more general case $`0<t_l,\eta ,\theta <1`$ the condition for maximum entropy needs to be evaluated numerically and it can be found that for realistic values of $`t_0,t_l,\eta ,`$ and $`\theta `$ the maximum value of entropy is decreased and the maximum occurs for $`r<1/2`$. Thus, in a real device we can expect best performance of the detector when SVR is unbalanced in favor of port 4. For the parameters of our device we get the maximum value of entropy at $`r=0.446`$.

## 3 Experimental results

For an experimental test of the detector we have fed it with a source of faint laser pulses with Poissonian statistics and variable mean photon number. They have been obtained from a laser diode (SHARP LT015) yielding 4 ns pulses that have been subsequently attenuated by a digital variable attenuator (OZ Optics DA-100) to a single-photon level.

A typical time spectrum measured by the time-of-flight spectrometer at the mean photon number $`\mu 2.13`$ (see Eq. (10) below for the definition of $`\mu `$) is shown in Fig. 2. Multiple peaks can be resolved at time distances of about 60 ns. The background between the peaks is caused by afterpulses, i.e. false detections that occur after dead time of the detector. This is why they are not present between the first two peaks. In our detector, the total afterpulse probability is of the order of $`p_{ap}^{peak}8\times 10^3`$. The envelope of the afterpulse probabilities coming from individual peaks then forms the background pattern visible in Fig. 2. At low intensities the probability of detection due to an afterpulse is considerably lower so that the main source of noise are the dark counts of the detector whose probability is $`2\times 10^7`$ per 5 ns bin.

It is important to know the losses in the device. The total transmission of the device may be evaluated using Eqs. (1) and (2) as

$$T=\underset{k=1}{\overset{\mathrm{}}{}}h_k=\eta t_0\left[\frac{\theta (t_{13}t_{24}t_{14}t_{23})}{t_{24}}\frac{t_{14}t_{23}\theta }{t_{24}(t_lt_{24}\theta 1)}\right].$$
(7)
The quantity $`T/\eta `$ can be measured directly and we have found $`T/\eta =0.78\pm 0.01`$. Since $`\eta `$ is specified by the manufacturer of the detector, the total transmission $`T`$ of the device is known. However, upon inspecting Eq. (7) we can get more insight into the structure of the losses.

We have observed experimentally that the total transmission $`T`$ depends only weakly on the SVR setting. We therefore again replace the four $`t_{ij}`$ coefficients by a single parameter $`r`$ and arrive at a simplified expression

$$T\frac{\eta t_0(2t_l\theta 1)}{t_l}\frac{\eta t_0(t_l\theta 1)^2}{t_l(rt_l\theta 1)}.$$
(8)
Since the dependence of $`T`$ on $`r`$ is weak, we may deduct that $`T`$ is given dominantly by the first term in Eq. (8). The value of $`\theta `$ was found in an independent measurement to be $`\theta =0.955`$. The transmission of the fiber loop $`t_l`$ can be found from the measured values of the normalized channel probabilities $`H_k=h_k/T`$ ($`_{k=1}^{\mathrm{}}H_k=1`$). In particular, the ratio $`H_{k+1}/(H_kH_1)`$ evaluates to

$$\frac{H_{k+1}}{H_kH_1}2\theta t_l1,$$
(9)
where only the first term of Eq. (8) has been used. From experimental results we find the value $`H_{k+1}/(H_kH_1)0.80`$ to be almost independent of $`r`$ and $`k`$ (thus justifying the approximations used). From Eq. (9) we then get $`t_l0.94`$ (0.27 dB). Finally, the value of input coupling transmission $`t_0`$ is evaluated from Eq. (8) as $`t_00.92`$ (0.36 dB). These values are rather low and are caused mainly by the wear-out of the connectors used in our laboratory setup. The application of new connectors or fused fibers would further improve the performance of the device.

The division ratio of the SVR is set by a micrometer screw. The distribution of detection probabilities to the first six channels based on the SVR position is shown in Fig. 3. The small white area in the top part of the plot corresponds to the sum of the higher channels ($`k>6`$). We can see that various settings are possible. We can, e.g., set the device to the regime of one dominant channel with other channels weak but a moderate number of them (5 or 6) nonnegligible (their weight is above 1 %; see positions 5-10 in Fig. 3), or to a different regime of fewer but relatively balanced channels (see, e.g., position 19 in Fig. 3 where the first three channels are in the ratio 39%:42%:13%).

Let us now have a look at the capability to resolve multi-photon states using our device. We characterize the measured Poissonian signal by the mean photon number $`\mu `$. The mean photon number $`\mu `$ is obtained from the measured probability of detection $`p`$ assuming detection of the whole Poissonian signal by the detector with efficiency $`T`$, for which $`p=_{i=1}^{\mathrm{}}p_n=1p_0=1\mathrm{exp}(T\mu )`$, i.e.

$$\mu =\mathrm{ln}(1p)/T.$$
(10)
Our measured signal had $`\mu `$ equal to 4.26 photons per pulse. The probability of vacuum state detection $`p_0`$, single-channel detection $`p_1`$, and multichannel detection $`p_M=_{k=2}^{\mathrm{}}p_k`$ were extracted from the measured data. At each setting, the probabilities were obtained by detecting large number of laser pulses ($`10^5`$). We characterize the fraction of multi-photon states detected by our detector with another quantity, namely the multi-photon content $`c_M`$, $`c_M=p_M/(p_1+p_M)`$ ($`c_M`$ gives the probability that non-vacuum pulses contain more than 1 photon) . Figure 4 shows the detected multi-photon content $`c_M`$ at the mean photon number $`\mu =4.26`$ photons per pulse. For comparison, a real (theoretical) multi-photon content $`c_M`$ in a Poissonian laser pulse in front of the device is shown as well ($`c_M`$ is not constant for different SVR settings because the laser intensity slightly changed from one setting to another). The highest values of $`c_M`$ have been observed for SVR positions 16-20 with maximum at position 17. The value of $`c_M`$ grows with the number of channels considered. Due to losses in the device and limited quantum efficiency of the detector, the measured multi-photon content is lower than the real (theoretical) multi-photon content. Nevertheless, the value of the measured multi-photon content was lower by less than 4 % in comparison with the value of the real (theoretical) multi-photon content at best performance of the device (SVR position 17, 15 channels) for this pulse-energy level. It is also worth noting that the use of a large number of channels is useful mainly in the regime with the first dominant channel while it brings only a little advantage in case of fewer but relatively balanced channels. Values of entropy for each setting of the SVR are also plotted in Fig. 4 and they clearly show that entropy of the multichannel detector is the right indication of the multi-photon-resolution capability. The optimum performance of the device is achieved close to the value $`r=0.453`$ as estimated theoretically in the preceding section.

The performance of the detector in identification of the multi-photon content over a wide range of the mean photon number $`\mu `$ of the input Poissonian state is shown in Fig. 5. The measured values of the multi-photon content $`c_M`$ (black rectangles) follow closely the (dashed) curve of the real (theoretical) multi-photon content of the input state. The ratio of the measured values of the multi-photon content $`c_M`$ to the real (theoretical) ones shown in Fig. 5 determines the efficiency of identification of multi-photon states. It can be seen in Fig. 5 that the ratio exceeds 60% for weak coherent states and reaches almost unity for strong coherent states.

Reduction of the multi-photon content $`c_{M,\mathrm{in}}`$ of a beam can be reached in the following scheme. A source provides perfectly correlated photon pairs (e.g., such pairs are generated by spontaneous parametric downconversion). One of the correlated beams is postselected by the measurement on the other beam using our fiber-loop detection device. Postselection process provides a beam with lower values of the multi-photon content $`c_{M,\mathrm{out}}`$. The ratio $`w_M`$ ($`w_M=c_{M,\mathrm{out}}/c_{M,\mathrm{in}}`$) then determines efficiency of this reduction. The dependence of the ratio $`w_M`$ on the mean photon number $`\mu `$ in Fig. 5 shows that the value of the ratio $`w_M`$ can be reduced down to 40 % for input Poissonian states with $`\mu 5`$. The most interesting input states for applications utilizing multi-photon reduction are those containing a more-or-less balanced mixture of vacuum, single-photon and multi-photon contributions. The reason is that weak Poissonian states contain only a small fraction of multi-photon states (e.g., $`c_M0.005`$ for $`\mu =0.01`$) whereas the fraction of a one-photon state is practically negligible in strong coherent states. The reduction down to 40 % can be ideally reached for such states. We note that a real source would be characterized by higher values of the ratio $`w_M`$ due to imperfect coupling of photons in a pair .

It should be mentioned that in fact our measured value of the multi-photon content is slightly influenced by the noises in the device. The noise counts may cause multichannel detection even when only a single photon from the laser pulse has been detected. It is difficult to determine exactly the contribution of these false multi-detections. However, it is possible to estimate the upper limit of false detections due to afterpulses as $`p_M^{false}<p_1p_{ap}^{peak}q`$, where $`q`$ is a duty factor of the detector channels at the time-of-flight spectrometer (i.e., the ratio of the channel time window to the time distance between channels; $`q=0.17`$ in our case). Then the contribution to the multi-photon content due to afterpulses $`c_M^{false}<p_{ap}^{peak}q1.4\times 10^3`$. For very weak pulse energies false multiple detections might rather stem from dark counts of the detector. Nevertheless, the latter would become notable only at mean-photon-number levels below $`10^5`$ photon per pulse. Therefore, the influence of noises is negligible in our device.

The analyzed device might help in the preparation of photon-number-squeezed states by postselection from entangled photon pairs obtained by spontaneous parametric downconversion . Vacuum states can easily be filtered out similarly to the case of postselection with a single detector  and the fraction of multi-photon states can be significantly reduced using multichannel detection. The obtained states are potentially useful, e.g., for quantum-key distribution systems. Our device is a cheap, though not better as for the performance, alternative to the cascading detector  for this purpose. Other devices like those based on NV centers, quantum dots or using detectors with a superconducting microcalorimeter are quite demanding even for a common laboratory practice, at least at present.

In this work we use only the information whether single or multiple detections occurred regardless of the information in which channels the photons were detected. Since the channels exhibit different detection probabilities, there is an additional information which might be used. In general, it is possible to obtain the photon-number statistics of a detected quantum state to some extent. This is a topic of a forthcoming publication.

In conclusion, we have built a configurable multichannel detector using a single avalanche photodiode. We have characterized the best regimes for its operation as multi-photon resolving detector. Our analysis shows, that this device can efficiently distinguish multi-photon states from the one-photon state and the vacuum state over a wide range of energies of the input state.

## Acknowledgments

The authors acknowledge the support by the projects Research Center for Optics (LN00A015), CEZJ-14/98 and RN19982003013 supported by the Ministry of Education of the Czech Republic. They also thank Pavel Trojek for his contribution to the experimental setup.