text:
# Comparison between Theory and Some Recent Experiments on Quantum Dephasing

## I Maximum dephasing times

For all the samples reported in Refs.  the mean free path $`l`$ deduced from the values $`D`$ was 2 to 5 times shorter than both $`t`$ and $`w`$ for quasi-$`1d`$ samples and $`t`$ for four quasi-$`2d`$ samples of Ref. . Therefore, even though the condition $`lw,t`$ was not very well satisfied in these experiments, it would be more appropriate to use eq. (4) for our comparison. It is also easy to observe from eqs. (5) that, since $`l`$ is not drastically smaller than $`w`$ and $`t`$, our $`1d`$ and $`2d`$ expressions will yield the same order of magnitude estimates for $`\tau _\phi ^{max}`$ for all the samples in question (also see below).

The maximum dephasing times were reported in Ref.  for ten (six quasi-$`1d`$ and four quasi-$`2d`$) samples fabricated from the same material ($`AuPd`$) with practically the same resistivity. For all ten samples the maximum dephasing times were found to be nearly universal $`\tau _\phi ^{max}(0.8÷2)\times 10^{11}`$ sek independently of the sample geometry. At $`T80`$ mK these values yield $`T\tau _\phi 0.09÷0.22`$. Since both the material and its properties were the same an important test for our theory is to fit the data  for $`\tau _\phi ^{max}`$ to eq. (4) which only depends on the material parameters and not on geometry.

In experiments  the maximum dephasing times were found to be 2 to 3 orders of magnitude bigger as compared to those in Ref. . Furthermore, $`\tau _\phi ^{max}`$ for the silver sample  was about an order of magnitude higher than that for the copper sample  with similar (although not identical) parameters. Observing this difference the authors  suggested that the low temperature saturation of $`\tau _\phi `$ in disordered metal wires cannot be universal and can be material dependent.

For our comparison we will use the data for 8 samples  with nominally identical resistivity $`\rho 24`$ $`\mu \mathrm{\Omega }`$ cm and $`D1.5\times 10^3`$ m2/sek (samples C to F and H to K), 2 samples  ($`Ag`$ and $`Cu`$) and one sample ($`AuMSU`$) . Since all the maximum dephasing times  were nearly the same, for the sake of brevity we only quote an average value for $`\tau _\phi ^{max}`$ obtained by averaging over samples C to F and H to K. In order to estimate $`\tau _{\phi 0}`$ for $`Ag`$ and $`Cu`$ samples we used the parameters given in the Table 1 of Ref. . Experimentally observed values of $`\tau _\phi ^{max}`$ and our theoretical predictions for $`\tau _{\phi 0}`$ (eq. (4) with $`a_3=1`$) are summarized in the following Table:

| sample | $`\tau _\phi ^{max}`$ (sek) | $`\tau _{\phi 0}`$ (sek) |
| --- | --- | --- |
| C to F and H to K  (averaged) | 1.3$`\times 10^{11}`$ | 0.3$`\times 10^{11}`$ |
| $`Ag`$  | 10$`\times 10^9`$ | 2$`\times 10^9`$ |
| $`Cu`$ | 1.8$`\times 10^9`$ | 0.3$`\times 10^9`$ |
| $`AuMSU`$  | 8$`\times 10^9`$ | 4$`\times 10^9`$ |

These results demonstrate that for all the above samples our formula (4) with $`a_3=1`$ gives correct order-of-magnitude estimates for $`\tau _\phi ^{max}`$ even though these times differ by 2 to 3 orders of magnitude depending on the sample and experiment. Already this agreement is fairly good since in weak localization measurements $`\tau _\phi `$ is defined up to a numerical prefactor. Furthermore, in accordance with our expectations, in all cases our theoretical estimates involving the cutoff parameter $`a_3=1`$ are smaller than the experimental values by a numerical factor ranging from 2 to 6. This implies that by choosing the cutoff parameter in the range $`a_30.3÷0.6`$ one can exactly reproduce all the above experimental values for $`\tau _\phi ^{max}`$ from the same formula (4) which only depends on the material parameters, like $`\rho `$ and $`D`$.

One can also eliminate the ambiguity related to the cutoff parameter $`a_d`$ in (4) by calculating the ratio between maximum dephasing times for different samples. Making use of the experimental results  and our theoretical estimates based on eq. (4) one finds:

|  | experiment | theory (eq. (4)) |
| --- | --- | --- |
| $`\tau _{\phi 0}^{AuPd}/\tau _{\phi 0}^{Ag}`$ | $`1.3\times 10^3`$ | $`1.5\times 10^3`$ |
| $`\tau _{\phi 0}^{AuPd}/\tau _{\phi 0}^{Cu}`$ | $`0.7\times 10^2`$ | $`1\times 10^2`$ |
| $`\tau _{\phi 0}^{AuPd}/\tau _{\phi 0}^{AuMSU}`$ | $`1.6\times 10^3`$ | $`0.8\times 10^3`$ |

The agreement appears to be even much better than one could possibly expect given all experimental and theoretical uncertainties in determination of $`\tau _\phi `$. We emphasize that our theoretical predictions presented in both above Tables involve no fit parameters. Previously  we already reported on a similar agreement with the results . All these observations strongly favour the conclusion about universality of the low temperature saturation of $`\tau _\phi `$ in disordered conductors.

To complete this issue we mention that also our $`1d`$ and $`2d`$ expressions for the maximum dephasing time  give correct order-of-magnitude estimates for $`\tau _\phi `$ for the above samples. For instance, for $`2d`$ samples  the second eq. (5) yields (for $`a_2=1`$) $`\tau _{\phi 0}^{(2d)}0.6\times 10^{11}`$ sek. This estimate is only $`2÷3`$ times smaller than the experimental values. Again, although such an accuracy is already more than sufficient, the remaining difference can easily be removed simply by choosing $`a_20.3÷0.5`$.

For comparison, the authors  estimated the value of $`\tau _\phi `$ which follows from the theory  for their $`2d`$ samples as $`\tau _\phi 6\times 10^9`$ sek at $`T=1K`$. At the lowest measurement temperature $`T80`$ mK this estimate translates into $`\tau _\phi `$ in the range $`10^7`$ sek which is $`3.5÷4`$ orders of magnitude bigger than the values measured in Ref. .

## II Temperature Dependence

We now analyze the temperature dependence of $`\tau _\phi `$ detected in experiments . For quasi-$`1d`$ samples eq. (1) takes the form

$$\frac{1}{\tau _\phi }=\frac{1}{\tau _{\phi 0}}+\frac{\rho T}{R_qtw}\sqrt{2D\tau _\phi },$$
(6)
where $`R_q=\pi /2e^26453`$ $`\mathrm{\Omega }`$. An essential difference between this result and that of Ref.  is that in our case the temperature-independent contribution to the dephasing rate $`1/\tau _{\phi 0}`$ is not equal to zero. If one neglects this contribution, eq. (6) would immediately yield the standard result  $`1/\tau _{\phi AAK}^{(1d)}T^{2/3}`$. Similarly, in quasi-$`2d`$ samples one gets  $`1/\tau _{\phi AAK}^{(2d)}T\mathrm{ln}T`$. In order to test these predictions it is appropriate to present experimental data for $`\tau _\phi (T)`$ on a log-log plot, exactly as it was done in Refs. . From Ref.  one expects ($`i`$) to observe straight lines with the slope becoming steeper as one crosses over from quasi-$`1d`$ to quasi-$`2d`$ systems and ($`ii`$) for two quasi-$`1d`$ samples with the same $`\rho `$ and $`D`$ but with different cross-sections $`tw`$ the straight lines corresponding to $`T^{2/3}`$-dependence of $`\tau _\phi `$ should be parallel.

None of these two features was observed in the experiments . Just on the contrary, the opposite trend was detected: In all cases the slope of the curves decreased with increasing sample width $`w`$ (see Fig. 3 of Ref. ). The curves for quasi-$`1d`$ samples with different $`w`$ were not parallel (actually the data points for the sample F on Fig. 3 of  clearly saturate and can hardly be fitted to a straight line at all) and, on top of that, the $`T`$-dependence of the decoherence rate for quasi-$`2d`$ samples – instead of being linear in $`T`$ – was very weak already at $`T1K`$.

Now let us see if these features can be accounted for within our theory. Since the low temperature behavior of $`\tau _\phi `$ predicted in Ref.  is different from a power-law one  it appears to be more useful to replot the same data on a linear scale. Then, according to eq. (6) at sufficiently low $`T`$ one expects to observe a linear dependence of the decoherence rate on $`T`$ shifted upwards by the value $`1/\tau _{\phi 0}`$. For the same parameters $`\rho `$ and $`D`$ (and the same $`\tau _{\phi 0}`$) the slope $`d(1/\tau _\phi )/dT`$ for quasi-$`1d`$ samples should depend only on the cross section $`tw`$. In particular, since $`w`$ for the sample $`F`$ was reported in  to be 4 times larger than for the samples A to E, one expects $`d(1/\tau _\phi )/dT`$ to be 4 times smaller for the sample F.

The experimental data for the samples A, B and F  are presented in Fig. 1. We observe that ($`a`$) the dephasing rate $`1/\tau _\phi `$ increases linearly with $`T`$, ($`b`$) at $`T0`$ these linear dependencies extrapolate to a nonzero (and practically the same) value $`1/\tau _\phi `$ for all three samples and ($`c`$) the slope $`d(1/\tau _\phi )/dT`$ is indeed about 4 times smaller for the sample $`F`$ as compared to the samples A and B. All these features are fully consistent with eq. (6). The magnitude of the slope is several ($`4÷6`$) times larger than predicted by eq. (6). The origin of this difference is not clear to us at the moment, but possibly it can be due to non-uniformity of the wires. The data for quasi-$`2d`$ samples I and K  also show the trend qualitatively consistent with our eq. (1). Also in that case the temperature effect seems to be more pronounced than it is predicted within our theory.

We conclude that our theory  accounts for all the main experimental observations of Ref. .

Now we will address the results of Refs. . In Fig 3. of  the authors presented their data for $`Au`$, $`Ag`$ and $`Cu`$ samples with similar parameters. Both the magnitude and the temperature dependence of $`\tau _\phi `$ differ drastically for these three samples. This observation led the authors  to conclude that the behavior of $`\tau _\phi `$ and, in particular, its saturation at low $`T`$ may be material dependent.

Later it was demonstrated by the same group  that unusually low values of $`\tau _\phi `$ in gold  were most likely due to high concentration of magnetic impurities. A pure gold sample ($`AuMSU`$) was fabricated and similar behavior of $`\tau _\phi `$ was observed  as previously for $`Ag`$ sample : $`\tau _\phi `$ did not saturate down to $`T45÷50`$ mK on a log-log plot. At the same time it was confirmed  that also very pure $`Cu`$ samples showed saturation similarly to an earlier $`Cu`$ sample .

The material dependence of $`\tau _\phi `$ cannot be ruled out in general. However, with minimum efforts one can demonstrate that no material dependence needs to be assumed in order to quantitatively explain seemingly different behavior of $`Ag`$ (no saturation down to $`T49`$ mK) and $`Cu`$ (clear saturation already at $`T700`$ mK) samples .

Let us rescale eq. (6) as

$$1/\stackrel{~}{\tau }_\phi =A+B\stackrel{~}{T}\sqrt{\stackrel{~}{\tau }_\phi },$$
(7)
where we introduced dimensionless variables

$$\stackrel{~}{\tau }_\phi =\frac{\tau _\phi }{\tau _{\phi 0}},\stackrel{~}{T}=\frac{T}{T_0},T_0=\frac{R_qtw}{\rho \tau _{\phi 0}^{3/2}\sqrt{2D}}.$$
(8)
Here $`A=B=1`$ both for $`Cu`$ and $`Ag`$ samples provided one performs the scaling to the values $`\tau _{\phi 0}`$ and $`T_0`$ calculated for the respective sample. The temperature $`T_0`$ sets the scale at which thermal and quantum contributions to the dephasing rate become comparable. Making use of the parameter values for $`Cu`$ and $`Ag`$ samples one finds

$$T_0^{Cu}/T_0^{Ag}20.$$
(9)
Note, that this estimate is not sensitive to various ambiguities related, e.g., to the cutoff parameter $`a_d`$. Eq. (9) just follows from eqs. (4), (6) for the values of $`\rho `$, $`D`$, $`t`$ and $`w`$ presented in the Table 1 of Ref. . The estimate (9) remains unchanged if, instead of our eq. (4), one uses experimental values for $`\tau _\phi ^{max}`$ . Hence, the temperature range where quantum effects set in can easily vary by one order of magnitude or more for different samples even if their macroscopic parameters are similar .

What remains is to make both log-log and linear scale plots of eq. (7), see Fig. 2. In these plots we scaled both curves to the parameters of one ($`Cu`$) sample, i.e. $`A_{cu}=B_{Cu}=1`$, $`A_{Ag}=\tau _{\phi 0}^{Cu}/\tau _{\phi 0}^{Ag}0.2`$ and $`B_{Ag}=A_{Ag}^{1/2}T_0^{Cu}/T_0^{Ag}9.5`$. Making use of the values for $`\tau _\phi ^{max}`$  we find $`T_0^{Cu}730`$ mK and $`T_0^{Ag}40`$ mK. These numbers as well as the temperature dependence of $`\tau _\phi `$ are in the agreement with experimental findings . E.g. on the log-log plot we observe that, while $`Cu`$ sample shows clear saturation around $`\stackrel{~}{T}1`$ (i.e. $`T700`$ mK), practically no signs of saturation are visible for $`Ag`$ sample even at temperatures one order of magnitude lower. Similar analysis can be performed, e.g., for $`Au6`$  and $`AuMSU`$  samples (cf. right panel of Fig. 7 in Ref. ). In that case one can estimate $`T_0^{Au6}/T_0^{AuMSU}6`$.

It is obvious from eq. (6) that by a proper choice of the sample parameters one can tune the effective temperature range at which saturation feature occurs on the log-log plot. E.g. thinner samples with identical other parameters will show no saturation down to lower temperatures. The lack of saturation for such samples even at relatively low $`T`$ would, of course, not imply that the electron dephasing rate is zero at $`T=0`$.

Let us replot the ”non-saturating” data  on a linear scale. In Fig. 3 we collected the data points for different samples studied in different experiments ($`Au1`$, $`Au2`$, $`Au6`$ of Ref. , $`Ag`$ of Ref.  and $`AuMSU`$ of Ref. ). All these samples show essentially the same behavior at low temperatures. Both the slope $`d(1/\tau _\phi )/dT`$ and the values $`\tau _{\phi 0}`$ slightly differ from sample to sample, but these differences are unimportant and can easily be accounted for by different values of $`D`$ and $`tw`$ as well as by measurement uncertainties. Much more importantly, at low temperatures $`1/\tau _\phi `$ depends linearly on $`T`$ and clearly extrapolates to a nonzero value for all the samples presented in Fig. 3. This is exactly what one expects from our eq. (6). Also the similarity with the results  (our Fig. 1) is quite obvious.

In summary, we demonstrated that the experimental results  are in a quantitative agreement with our theory  which predicts the low temperature saturation of the electron dephasing time in diffusive conductors. No material dependence of $`\tau _\phi `$ needs to be assumed in order to explain the data  for the samples $`Ag`$, $`Cu`$ and $`AuMSU`$.


format: text

text:
# Asymmetric Beam Combination for Optical Interferometry

## 1 Introduction

Spatial filtering is key for precise calibration of atmospheric turbulence for ground-based optical interferometers. After sidereal delay compensation, the stellar light from each telescope is focused onto the tip of a single-mode fiber. Up to $``$78% of the energy (Shaklan & Roddier 1988) of a perfect telescope beam can be coupled into a single-mode fiber. In practice, this coupling efficiency is significantly less, depending on how well the incoming electric field distribution overlaps with the fiber beam profile. This is the act of “spatial filtering,” whereby only the “clean” part of each wavefront (one spatial mode) is transmitted by the single-mode fiber. In essence, phase irregularities of the incoming wavefronts are converted into amplitude fluctuations, trading a difficult to measure quantity for one easily monitored (e.g., Shaklan 1989).

The light from each fiber can then be efficiently and coherently combined in a balanced 2x2 coupler, a fiber optic component which mixes 50% of the light from each telescope coherently, equivalent to a beamsplitter. In some cases, approximately one-half the light is split-off from each fiber before the coupler to act as a “photometric tap,” a direct and real-time monitor of the coupling efficiency for each telescope. In this scheme the fringe amplitude measured at the outputs of the coupler can be directly calibrated based on the strength of the signals in the photometric taps (e.g., Coude Du Foresto et al. 1997).

There are two distinct ways that spatial filtering improves calibration. The most important effect is that the mean coupling into each fiber is a direct measure of atmospheric turbulence, with poor seeing resulting in low coupling efficiency. When seeing inevitably changes in the time between observing the target and calibrator, the coupling efficiency into the fibers will change, and this will be seen as a change in average signal level. However, since the interference efficiency in the 2x2 coupler is not a function of atmospheric conditions, the overall calibration of the system visibility is maintained. Heterodyne detection is another kind of spatial filtering which utilizes this principle, and is one reason why the Infrared Spatial Interferometer (e.g., Hale et al. 2000) can maintain accurate calibration despite the dearth of bright mid-infrared calibrators.

The second way spatial filtering can improve measurement precision is by correcting for short-term fluctuations in coupling efficiency. The fringe amplitude is proportional to the harmonic mean of the two telescope intensity levels (e.g., see Eq. 3), while the DC level is proportional to the sum of the two telescope intensity levels. This means that the fringe contrast, or visibility, decreases when the signal from one telescope is significantly greater than the other. While this effect on the system visibility can be statistically calibrated by using a calibrator star, the use of photometric taps allows each individual fringe measurement to be independently calibrated. This second use of spatial filtering allows high quality fringe visibility estimates with very few measurements, since one does not have to statistically average over these fluctuations (for specific application and discussion, see Coude Du Foresto et al. 1997). Note that this second aspect of spatial filtering is of limited utility for faint sources, since it can only be applied when the signal levels can be measured to reasonable accuracy within a single atmospheric coherence time.

Initial experiments using spatial filters in infrared interferometry have led to dramatic improvement in the precision of visibility measurements. Using flouride single-mode fibers with in-line 2x2 couplers and photometric taps, Perrin et al. (1999) have reported visibility measurements with 0.3% precision of bright M-giants in the K-band, 10 to 20$`\times `$ better than the precision reported using traditional techniques (e.g., Dyck et al. 1996; Nordgren et al. 1999).

In this short paper, I discuss the use of asymmetric 2x2 couplers which can permit the calibration improvements derived from photometric taps but without diverting valuable photons from the fringe amplitude measurement. I will detail the signal-to-noise ratio advantages in both read-noise and Poisson-noise limiting cases. Lastly, I will briefly comment on how these ideas can be applied to other beam combination schemes, such as when spatial filtering occurs after beam combination, for pair-wise and “all-in-one” combiners for more than two beams, for free-space combiners with no spatial filtering, and for fringe-tracking systems.

## 2 Balanced and Asymmetric Couplers

As mentioned in the Introduction, current fiber optic correlators, e.g. FLUOR instrument on the IOTA interferometer (Coude Du Foresto et al. 1997), VINCI on VLTI (Foresto, private communication), use a photometric tap which removes $``$50% of the light from each telescope before beam combination in a “Balanced” (symmetric 50/50) coupler (hereafter, “B-coupler”). This is schematically indicated in the top panel of Figure 1. However, this method is not efficient for resolved or faint sources since roughly half the photons are being diverted from the fringe measurement for the photometry; it is clearly futile to monitor photometric fluctuations when the fringes can not even be detected.

The bottom panel of Figure 1 shows a novel method for recovering the photometric signals without requiring dedicated photometric taps, thus allowing improved sensitivity while maintaining high calibration precision. At the center of the design is an “Asymmetric” coupler (hereafter, “A-coupler”). An A-coupler is much like a beamsplitter manufactured to reflect a different fraction than it transmits. It will be shown that the photometric signals from each telescope can be reconstructed through appropriate linear combinations of the interferometric signals after the A-coupler, assuming the properties of the coupler are well-known. While the fringe contrast from an A-coupler is slightly decreased by having unequal contributions from each telescopes, the signal strength is increased since no light is wasted on dedicated photometric taps.

For the following discussion and calculations I will refer to Figure 1 which labels the various input and outputs to the beam combiners under consideration. The incident photon fluxes (photons$``$s<sup>-1</sup>) coming from telescope 1 and 2 will be denoted as $`F_1`$ and $`F_2`$ respectively. In order to take into account the finite quantum efficiency $`\eta `$ of the detector, substitute $`\eta F_i`$ for $`F_i`$ in the subsequent formulas. Both combiners under consideration (“B-Coupler” and “A-Coupler”) have so-called “interferometric outputs,” labeled $`I_1`$ and $`I_2`$, while the “B-Coupler” also has two photometric outputs $`P_1`$ and $`P_2`$.

For this calculation we will assume idealized systems. The B-coupler scheme will divert 50% of the incident light into the photometric channels, and we will also assume the coupler is exactly balanced and loss-less. Under these assumptions the output signals can be expressed as:

$`P_1`$ $`=`$ $`{\displaystyle \frac{F_1}{2}}`$ (1)
$`P_2`$ $`=`$ $`{\displaystyle \frac{F_2}{2}}`$ (2)
$`I_1`$ $`=`$ $`{\displaystyle \frac{F_1}{4}}+{\displaystyle \frac{F_2}{4}}+2\sqrt{{\displaystyle \frac{F_1}{4}}{\displaystyle \frac{F_2}{4}}}\gamma (t)={\displaystyle \frac{F_1}{4}}+{\displaystyle \frac{F_2}{4}}+{\displaystyle \frac{1}{2}}\sqrt{F_1F_2}\gamma (t)`$ (3)
$`I_2`$ $`=`$ $`{\displaystyle \frac{F_1}{4}}+{\displaystyle \frac{F_2}{4}}2\sqrt{{\displaystyle \frac{F_1}{4}}{\displaystyle \frac{F_2}{4}}}\gamma (t)={\displaystyle \frac{F_1}{4}}+{\displaystyle \frac{F_2}{4}}{\displaystyle \frac{1}{2}}\sqrt{F_1F_2}\gamma (t)`$ (4)
In the last two equations, $`\gamma (t)`$ is the mutual coherence function and encodes the fringe visibility, the quantity we wish to measure. Usually $`\gamma `$ is temporally modulated by adjusting the relative path lengths in the two arms of the interferometer; this will be discussed further when we calculate the signal-to-noise ratios for the different coupler designs.

The A-coupler will be defined so that the fraction $`\alpha `$ of $`F_1`$ is mixed with $`\beta `$ of $`F_2`$ for interferometric output channel $`I_1`$, while the fraction $`\beta `$ of $`F_1`$ is mixed with fraction $`\alpha `$ of $`F_2`$ for interference in $`I_2`$. If the coupler is loss-less, then $`\beta =1\alpha `$, but I will keep the $`\alpha `$ and $`\beta `$ notation for equation simplicity. To summarize, the outputs of the A-coupler are given as:

$`I_1`$ $`=`$ $`\alpha F_1+\beta F_2+2\sqrt{\alpha F_1\beta F_2}\gamma (t)`$ (5)
$`I_2`$ $`=`$ $`\beta F_1+\alpha F_22\sqrt{\alpha F_1\beta F_2}\gamma (t)`$ (6)
### 2.1 Determining photometric signals with A-coupler

While determination of $`F_1`$ and $`F_2`$ is straightforward for the case of the B-coupler (see Eqs. 1 & 2), let us now discuss how to retrieve this same information using the A-coupler. We can split the $`I_i`$ signal into the sum of an incoherent and a coherent term, $`I=I_{\mathrm{incoherent}}+I_{\mathrm{coherent}}`$. The coherent part of this equation contains the interferometric “fringe” signal, and is usually separated from the incoherent part by temporal modulation of $`\gamma (t)`$ at a frequency $`\nu `$ greater than $`\frac{1}{\tau _0}`$, where $`\tau _0`$ is the coherence time of the atmosphere. In order to estimate the photometric flux from each telescope at the same time as we are detecting the fringe we must be able to estimate the incoherent part of the $`I`$-channels. For the rest of this discussion, we will denote this incoherent component of the $`I`$-channels as $`I_i`$, where the time-averaging is sufficient to remove the fringe component but not atmospheric fluctuations; this is equivalent to setting $`\gamma (t)=0`$ in Eqs. (5) and (6). The following SNR calculations for the A-coupler do assume one can separate the coherent and incoherent parts perfectly.

In order to estimate $`F_i`$ using the A-coupler scheme we must use linear combinations of $`I_1`$ and $`I_2`$. From inspection of Eqs. (5) and (6), we see that:

$`F_1`$ $`=`$ $`{\displaystyle \frac{\alpha I_1\beta I_2}{\alpha ^2\beta ^2}}`$ (7)
$`F_2`$ $`=`$ $`{\displaystyle \frac{\beta I_1\alpha I_2}{\beta ^2\alpha ^2}}`$ (8)
Note that these equations become indeterminate when $`\alpha =\beta `$, corresponding to the balanced coupler case, $`\alpha =\frac{1}{2}`$. This tells us that $`\alpha `$ must be significantly larger than $`\beta `$ in order to avoid “amplifying” the measurement uncertainty when estimating $`F_i`$.

## 3 Signal-to-Noise Ratio Calculation

In this section, I will compare the Signal-to-Noise Ratio (SNR) for the proposed asymmetric coupler system (A-coupler) versus the traditional setup using a 50% photometric tap and balanced couplers for the interference (B-coupler).

### 3.1 Read-Noise Limit (Infrared)

In the read-noise limited situation, we can imagine each output channel having the same read-noise $`\sigma _{\mathrm{read}}`$. Since the B-coupler scheme has twice the number of output channels, we see that more read-noise is introduced and we can expect degraded SNR. In the case of background-limited observations, the considerations are similar to the read-noise dominated case, since read-noise and background noise are both independent of signal strength and average as the root of the number of independent measurements. However, the discussion is complicated since background radiation can originate both before and after the combiner. The background-limited case is discussed briefly in §3.3.

We proceed to calculate the SNR for the read-noise limited case.

#### 3.1.1 Photometric Channels

The amplitudes of the photometric and interferometric signals were determined in §2 (Eq.1-8). For the B-coupler, the noise in each photometric channel is obviously $`\sigma _{\mathrm{read}}`$ per sample. For a total integration time of $`N\mathrm{\Delta }t`$ (N samples each with integration time of $`\mathrm{\Delta }t`$), it is trivial to calculate the estimated signal-to-noise ratio (SNR) and this result has been placed in Table 1.

For the A-coupler we form the linear combination described in §2.1, Eqs. (7) and (8). In estimating $`F_i`$ we must combine the errors from $`I_1`$ and $`I_2`$ in quadrature, weighted by the appropriate factors. Let us also assume here we are averaging $`N`$ separate samples, each with integration time of $`\mathrm{\Delta }t`$, to take into account averaging over the fringe modulation (see §2.1). In this case, the SNR for $`F_1`$ can be written as:

$`\mathrm{SNR}F_1{\displaystyle \frac{\mathrm{Signal}}{\mathrm{Noise}}}`$ $`=`$ $`{\displaystyle \frac{F_1N\mathrm{\Delta }t}{\frac{\sqrt{N}\sigma _{\mathrm{read}}\sqrt{\alpha ^2+\beta ^2}}{|\alpha ^2\beta ^2|}}}`$ (9)
$`=`$ $`{\displaystyle \frac{|\alpha ^2\beta ^2|F_1\sqrt{N}\mathrm{\Delta }t}{\sigma _{\mathrm{read}}\sqrt{\alpha ^2+\beta ^2}}}`$ (10)
A similar derivation for $`F_2`$ completes the first half of Table 1.

#### 3.1.2 Interferometric Channels

It is more straightforward to calculate the SNR for the fringe amplitude. The signal strength is simply equal to the coherent term of $`I_i`$, while the noise per sample is $`\sigma _{\mathrm{read}}`$. Eq. (5) can be used to calculate the signal-to-noise ratio for Fringe 1 as follows:

$`\mathrm{SNR}\mathrm{Fringe}1{\displaystyle \frac{\mathrm{Coherent}\mathrm{part}\mathrm{of}I_1}{\mathrm{Read}\mathrm{Noise}}}`$ $`=`$ $`{\displaystyle \frac{2N\mathrm{\Delta }t\sqrt{\alpha F_1\beta F_2}|\gamma |}{\sigma _{\mathrm{read}}\sqrt{N}}}`$ (11)
$`=`$ $`{\displaystyle \frac{2\sqrt{N}\mathrm{\Delta }t}{\sigma _{\mathrm{read}}}}\sqrt{\alpha F_1\beta F_2}|\gamma |`$ (12)
Substituting $`\alpha =\beta =\frac{1}{4}`$, this result can be applied for the case of B-coupler. These results now appear in Table 1.

The SNR advantage of the A-coupler approach, $`\mathrm{SNRA}=\frac{\mathrm{A}\mathrm{Coupler}\mathrm{SNR}}{\mathrm{B}\mathrm{Coupler}\mathrm{SNR}}`$, has been included in the last column of Table 1 and will discussed in §3.3.

### 3.2 Photon-Noise Limit (Visible)

In the photon-limited regime, there is no explicit punishment for using extra detectors or for multiple reads, and therefore we do not expect the A-coupler advantage to be as significant here. We will keep all the other notation the same. The formulation for SNR for A- and B-couplers developed in the last section still applies, except we must use $`\sigma _{\mathrm{photon}}`$ instead of $`\sigma _{\mathrm{read}}`$. For pure Poisson noise, the rms fluctuation in the measurement of $`F_1N\mathrm{\Delta }t`$ photons is simply $`\sqrt{F_1N\mathrm{\Delta }t}`$ photons. We now proceed to create Table 2 containing the SNR for the B-coupler and A-coupler under photon-noise dominated conditions.

#### 3.2.1 Photometric Channels

For the B-coupler, the noise in each photometric channel is equal to the square root of the number of detected photons. For a total integration time of $`N\mathrm{\Delta }t`$, the estimated SNR is straightforward to to derive, and appears in Table 2.

As before, use of the A-coupler requires the formation of linear combinations (Eqs. 7 and 8) of $`I_1`$ and $`I_2`$ to estimate $`F_i`$. To correctly propagate the measurement uncertainties of $`I_1`$ and $`I_2`$ into our estimate of $`F_i`$, we combine the Poisson errors from $`I_1`$ and $`I_2`$ in quadrature, weighted by the appropriate factors. As before, we average $`N`$ separate samples, each with integration time of $`\mathrm{\Delta }t`$, to take into account averaging over the fringe modulation. In this case, the SNR for $`F_1`$ can be written as:

$`\mathrm{SNR}F_1{\displaystyle \frac{\mathrm{Signal}}{\mathrm{Noise}}}`$ $`=`$ $`{\displaystyle \frac{F_1N\mathrm{\Delta }t}{\frac{\left[\left(\alpha \sqrt{I_1N\mathrm{\Delta }t}\right)^2+\left(\beta \sqrt{I_2N\mathrm{\Delta }t}\right)^2\right]^{\frac{1}{2}}}{|\alpha ^2\beta ^2|}}}`$ (13)
$`=`$ $`{\displaystyle \frac{F_1\sqrt{N\mathrm{\Delta }t}}{\frac{\sqrt{\alpha ^2(\alpha F_1+\beta F_2)+\beta ^2(\beta F_1+\alpha F_2)}}{|\alpha ^2\beta ^2|}}}`$ (14)
$`=`$ $`{\displaystyle \frac{|\alpha ^2\beta ^2|F_1\sqrt{N\mathrm{\Delta }t}}{\sqrt{\alpha ^2(\alpha F_1+\beta F_2)+\beta ^2(\beta F_1+\alpha F_2)}}}`$ (15)
I have not collected terms further in order to make it more apparent that for the nominal case of a loss-less coupler ($`\alpha +\beta =1`$) and for equal fluxes from each telescope ($`F_1=F_2`$), then $`\mathrm{SNR}F_1=\frac{|\alpha ^2\beta ^2|\sqrt{F_1N\mathrm{\Delta }t}}{\sqrt{\alpha ^2+\beta ^2}}`$, which has the same $`(\alpha ,\beta )`$ dependences as for the read-noise limited case (see Eq. 10). The calculation for $`\mathrm{SNR}F_2`$ follows directly, and we place these results in Table 2.

#### 3.2.2 Interferometric Channels

As for the read-noise limited case, the SNR of the fringe amplitude is more straightforward to derive than for the photometric signals. The signal strength is simply equal to the coherent term of $`I_i`$, while the Poisson noise should be calculated based on the incoherent term $`I_i`$, corresponding to the average flux. Eq. (5) can be used to calculate the signal-to-noise ratio for Fringe 1 as follows:

$`\mathrm{SNR}\mathrm{Fringe}1{\displaystyle \frac{\mathrm{Coherent}\mathrm{part}\mathrm{of}I_1}{\mathrm{Poisson}\mathrm{Noise}\mathrm{of}I_1}}`$ $`=`$ $`{\displaystyle \frac{2N\mathrm{\Delta }t\sqrt{\alpha F_1\beta F_2}|\gamma |}{\sqrt{(\alpha F_1+\beta F_2)N\mathrm{\Delta }t}}}`$ (16)
$`=`$ $`{\displaystyle \frac{2\sqrt{N\mathrm{\Delta }t}\sqrt{\alpha F_1\beta F_2}|\gamma |}{\sqrt{(\alpha F_1+\beta F_2)}}}`$ (17)
This result also applies for the B-combiner, substituting $`\alpha =\beta =\frac{1}{4}`$.

As for the photometric SNR using the A-coupler, this result (Eq. 17) can be greatly simplified in the case of a loss-less combiner and equal telescope fluxes, i.e. $`\mathrm{SNR}\mathrm{Fringe}1=2\sqrt{F_1N\mathrm{\Delta }t}\sqrt{\alpha \beta }|\gamma |`$, which has the same $`(\alpha ,\beta )`$ dependence as the read-noise limited case.

### 3.3 Comparison

In order to facilitate comparison between the performance of the A-coupler and the B-coupler, we will make a few simplifying assumptions. We will assume that $`F_1=F_2`$ and that the A-coupler is loss-less ($`\alpha +\beta =1`$). In addition, we will introduce a new quantity, $`A=\frac{\alpha }{\beta }`$, the “Asymmetry Factor,” defined greater than or equal to unity. Note that for a balanced system, the Asymmetry Factor is equal to Unity.

For this simplified configuration, the SNRs for the two photometric determinations are identical, as are the SNRs for Fringes 1 and Fringes 2. Therefore, we will stop differentiating between the two and speak only of the SNR of the photometric signal (i.e., the flux amplitude) and of the interferometric signal (i.e., the fringe amplitude). Table 3 contains the SNR advantage of the A-coupler over the B-coupler as a function for Asymmetry Factor ($`A`$) for the read-noise and photon-noise limit cases based on results contained in the Tables 1 and 2. Note that the SNR advantage between the read-noise-limited and photon-limited regimes differs by exactly $`\sqrt{2}`$ for both the flux and the fringe amplitude estimators.

It is interesting to note that the background-limited case falls somewhere in-between these two limits. The photometric channels and the non-unity coupling efficiency will introduce background radiation at the temperature of the fiber optic components. For a free-space system, this is equivalent to the fact that you can not “split-off” a fraction of a beam without introducing new lines-of-sight to the detector not coupled to the sky; thus additional background degrades the measurement. However if the combiner is cryogenically cooled, the background-limited case will be identical to the photon-noise limited case. Alternatively, if emission from the combiner dominates sky+telescope emission then the background-limited case will share the same SNR as the read-noise limited case. In practice, the result will intermediate.

If we assume the SNRs of the flux amplitude and fringe amplitude are equally important for estimating the fringe visibility (or contrast), which is probably only true for relatively bright sources, then we can estimate the combined SNR Advantage (SNRA) of the A-coupler over the B-coupler by applying equal weights and adding the errors in quadrature. Mathematically, this is done as follows:

$$SNRA_{\mathrm{Visibility}}\frac{\sqrt{2}}{\sqrt{\frac{1}{SNRA_{\mathrm{Flux}}^{}{}_{}{}^{2}}+\frac{1}{SNRA_{\mathrm{Fringe}}^{}{}_{}{}^{2}}}}$$
(18)
These results are shown graphically in Figures 2 and 3. In the read-noise limited case, we see that the A-coupler design is superior for Asymmetry factors greater than 1.78 and less than 25.63, with maximum improvement of 53% for $`A=4.61`$. In the photon-noise limit, the A-coupler is only slightly superior between Asymmetry factors of 2.82 and 8.71, with a maximum improvement of 8.2%, also at $`A=4.61`$.

Strictly speaking, the photometric signal can be averaged for longer than the fringes since the fringes must be modulated faster than the atmospheric time scale. This means that the SNR of the B-coupler setup can be improved by diverting a smaller percentage of the flux to the photometric taps; the optimal ratio would depend on various instrumental and observing factors, such as the inherent fringe visibility of your source and the temporal coherence time of the atmosphere. However by the same argument, the SNR of the A-coupler setup can be also improved by using a smaller asymmetry factor. These details have been neglected in the above calculation, but are not expected to qualitatively change the results.

## 4 Discussion

We see that the A-coupler offers a significant improvement for the read-noise limited case, and a slight improvement in the photon-noise limit. This method will work with a FLUOR-type system where beam combination occurs in a fiber optic coupler or integrated optic (Berger et al. 1999). Fringe-tracking systems can benefit from this scheme as well, since the mean signal in each interferometric channel can be separately determined by binning all data in each fringe-modulation period, but this does require both interferometric channels to be measured simultaneously. Asymmetric beam combination is also useful when there is no spatial filtering. Other effects can cause modulation of the beam intensities from telescopes, including scintillation and guiding errors, and the methods outlined here can be used to correct for these effects.

For interferometric arrays with three or more telescopes, other possibilities exist for determining the individual telescope photometric strengths. The A-coupler scheme can be generalized for “all-in-one” combiners (e.g., Baldwin et al. 1994), for instance using four asymmetric couplers allows one to differentiate the four input signals using the four output signals. Interestingly for the three-telescope (or more) pair-wise combination scheme, the photometric signals can be determined even with standard, symmetric couplers using linear combinations similar to the ones derived here.

Another advantage to the A-coupler design is decreased complexity of the combiner, although the data analysis will be more complicated. This is especially important in the visible, because of the expensive cost of high quantum efficiency, photon-counting detectors, such as avalanche photo-diodes (APDs). One apparent disadvantage to an A-coupler design is that one can no longer subtract the two interferometric signals to remove “common-mode” fluctuations such as scintillation. However, when the fringes are modulated faster than scintillation time scales, these variations can be completely corrected by forming the appropriate linear combinations as already discussed herein.

One caveat should be mentioned here. When spatial filtering occurs after (free-space) beam combination (e.g., Colavita et al. 1999), these methods may not work well. The reason is that optical misalignments leading to the single-mode fibers can cause decorrelations between the coupling efficiencies on either side of the free-space beamsplitter. For instance, a wavefront tilt introduced by the atmosphere could cause an increase in coupling efficiency for one fiber, but a decrease for the other fiber if the optical axes of both fibers are not aligned identically with the telescope axis. Hence, the photometric estimators will not be a fixed linear combination of the interferometric outputs, degrading the utility of the technique.

An experiment is underway at the Harvard-Smithsonian Center for Astrophysics to test these ideas, using a 780 nm asymmetric coupler, asymmetry factor 2.3. If lab tests prove satisfactory, the equipment will be incorporated into IOTA (Carleton et al. 1994), utilizing much of the existing visible light detection system for measurements of stellar diameters. The ultimate goal of this experiment is to yield high precision diameters measurements ($`<<`$1%), allowing the pulsation amplitude to be measured for calibration of the Cepheid distance scale.

## 5 Conclusions

Asymmetric beam combination in optical interferometry can improve sensitivity while reducing combiner complexity and cost, especially in the read-noise limited regime encountered for most infrared observing.

I appreciate discussions with Rafael Millan-Gabet, Jean-Phillipe Berger, Wes Traub, and other IOTA group members, and thank Peter Tuthill and David Mozurkewich for comments on the manuscript. Nat Carleton pointed out how photometric correction for a three-telescope, pair-wise combination scheme can be achieved with symmetric couplers, using appropriate linear combinations of the output channels. After the manuscript was submitted, the author was informed that a similar idea was proposed in 1995 by Jean Gay of Observatoire de la Cote d’Azur in Nice. JDM acknowledges support from a Center for Astrophysics Fellowship at the Harvard-Smithsonian Center for Astrophysics.


text:
# References

Photon production in relativistic heavy-ion collisions using rates with two-loop calculations from quark matter

Dinesh Kumar Srivastava

Variable Energy Cyclotron Centre,
1/AF Bidhan Nagar, Calcutta 700 064, India

Erratum: Eur. Phy. J. C 10 (1999) 487.

It has been shown recently  shown that the values of $`J_T`$ and $`J_L`$ which appear in Eq.(2) and (5) of the above paper and which are taken from the work of Aurenche et al are too large by a factor of 4. Correcting for this changes the Fig. 1–4 (see below). The basic result of the paper, that the emissions from the quark matter can outshine those from the hadronic matter when the photon rates upto two-loop level are used, remains valid though the range of $`p_T`$ over which it happens is reduced to $`p_T<`$ 0.5, 1, and 2 GeV/$`c`$ respecitively at SPS, RHIC, and LHC energies.


text:
# Efficient Monte Carlo Simulation Methods in Statistical Physics

## 1 Statistical Physics

Statistical mechanics  is a theory developed at the end of nineteenth century to deal with physical systems from an atomistic point of view. In principle the properties of bulk matter, which may contain $`10^{22}`$ atoms, can be worked out from the motion of atoms following the basic equations of Newtonian mechanics or quantum mechanics. However, such detailed information is not available or not really necessary. A probabilistic point of view, in the form of statistical ensemble is taken. The theory is extremely economical and successful in dealing with equilibrium problems.

There are a number of equivalent formulations of the theory. But statistical mechanics in a nutshell is the following concise formula that connects statistical mechanics with thermodynamics. First the partition function is defined as

$$Z=\underset{X}{}\mathrm{exp}\left(\frac{E(X)}{kT}\right).$$
(1)
where $`X`$ is a “state” of the system; the summation is carried over all possible states. $`T`$ is absolute temperature and $`k`$ is the Boltzmann constant ($`1.38\times 10^{23}`$ Joule/Kelvin). The free energy of the system at temperature $`T`$ is given by

$$F=kT\mathrm{ln}Z.$$
(2)
Free energy is a useful thermodynamic quantity in dealing with phase transitions. Other macroscopic observables are averages of the corresponding microscopic quantities over the Boltzmann-Gibbs weight $`\mathrm{exp}\left(E(X)/kT\right)`$,

$$g(X)=\frac{1}{Z}\underset{X}{}g(X)\mathrm{exp}\left(E(X)/(kT)\right).$$
(3)
The following three quantities are perhaps the most important ones, internal energy $`U=E`$, heat capacity $`C=dU/dT=\left(E^2E^2\right)/(kT^2)`$, and entropy $`S=(UF)/T=k\mathrm{ln}p(X)`$, where the average is over the Gibbs probability density, $`p(X)=\mathrm{exp}\left(E(X)/kT\right)/Z`$.

## 2 Monte Carlo Method

The essential task in statistical mechanics is to do the multi-dimensional integrations (for continuous systems) or summations (for discrete systems). To fix the notation and language, we briefly introduce the basics of Monte Carlo method . Consider the computation of a statistical average,

$$G=_\mathrm{\Omega }g(X)p(X)𝑑X=g,$$
(4)
where the probability density obeys $`p(X)0`$, and $`_\mathrm{\Omega }p(X)𝑑X=1`$. Suppose that we can generate samples $`X_i`$ according to the probability $`p(X)`$, then the integral can be estimated from an arithmetic mean over the samples,

$$G_M=\frac{1}{M}\underset{i=1}{\overset{M}{}}g(X_i).$$
(5)
The random variable $`G_M`$ has mean $`G`$ and standard deviation $`\sqrt{\tau (gG)^2/M}`$, $`\tau `$ being decorrelation time. Thus in the limit of large number of samples, the estimate converges to the exact value.

The most general method to generate $`X`$ according to $`p(X)`$ is given by a Markov chain  with a transition probability $`P(X^{}|X)=W(XX^{})`$, satisfying the conditions stated below. This is the probability of generating a new state $`X^{}`$ given that the current state is $`X`$. Such process converges if

$$W^n(X^{}X)>0,\text{for all }X^{}\text{ and }X\text{, and }n>n_0.$$
(6)
The equilibrium distribution of the Markov chain satisfies

$$p(X)=\underset{X^{}}{}p(X^{})W(X^{}X).$$
(7)
In constructing a Monte Carlo algorithm, it is convenient to consider a much stronger condition, the detailed balance

$$p(X^{})W(X^{}X)=p(X)W(XX^{}).$$
(8)
One of the most famous and widely used Monte Carlo algorithms is the Metropolis importance sampling algorithm . It takes a simple choice of the transition matrix:

$$W(XX^{})=S(XX^{})\mathrm{min}(1,\frac{p(X^{})}{p(X)})$$
(9)
where $`XX^{}`$, and $`S`$ is a conditional probability of choosing $`X^{}`$ given that the current value is $`X`$, and it is symmetric, $`S(XX^{})=S(X^{}X)`$. Usually $`S(XX^{})=0`$ unless $`X^{}`$ is in some “neighborhood” of $`X`$. The diagonal term of $`W`$ is fixed by the normalization condition $`_X^{}W(XX^{})=1`$.

## 3 Cluster Algorithms

In order to facilitate the discussion we first introduce the Ising model. Ising model is an interacting many-particle model for magnets. A state consists of a collection of variables $`\sigma _i`$ taken on only two possible values $`+1`$ and $`1`$, signifying the spin up and spin down states. The spins are on a lattice. The energy of the state $`\sigma =\{\sigma _i\}`$ is given by

$$E(\sigma )=J\underset{i,j}{}\sigma _i\sigma _j,$$
(10)
where $`J`$ is a constant which fixes the energy scale, the summation is over the nearest neighbor pairs. When temperature $`T`$ is specified, the states are distributed according to

$$p(\sigma )=Z^1\mathrm{exp}\left(\frac{E(\sigma )}{kT}\right).$$
(11)
In a local Monte Carlo dynamics (Metropolis algorithm), one picks a site $`i`$ at random, i.e., choosing a site with equal probability (this specifies or realizes $`S`$). Then, the energy increment if the spin is flipped, $`\sigma _i\sigma _i`$, is computed, which has the result $`\mathrm{\Delta }E=2J\sigma _i_j\sigma _j`$. The flip is accepted with probability $`\mathrm{min}(1,\mathrm{exp}(\mathrm{\Delta }E/kT))`$. If the flip is rejected, the move is also counted and the state remains unchanged. One Monte Carlo step is defined as $`N`$ moves (trials) for a system of $`N`$ spins, such that each spin is attempted to flip once on average.

The local algorithm of Metropolis type has some salient features: (1) it is extremely general. Less assumption is made of the specific form of probability distribution. (2) Each move involves $`O(1)`$ operations and $`O(1)`$ degrees of freedom. (3) The dynamics suffers from critical slowing down. The correlation time $`\tau `$ diverges as a critical temperature is approached. We shall elaborate more on this in the following.

The statistical error, using estimator Eq. (5), is given by

$$ϵ\sigma _G\sqrt{\frac{\tau }{M}},$$
(12)
where $`\mathrm{Var}(G_1)=\sigma _G^2=g^2g^2`$ is the variance of the observable, $`M`$ is the number of Monte Carlo steps, and $`\tau `$ is the decorrelation time. We can take the point of view that the above equation defines $`\tau `$, i.e.,

$$\tau =\underset{M\mathrm{}}{lim}\frac{M\mathrm{Var}(G_M)}{\mathrm{Var}(G_1)}.$$
(13)
Perhaps it is appropriate to call $`\tau `$ decorrelation time, since $`(\tau 1)/2`$ is sometimes called correlation time in the literature. The decorrelation time $`\tau `$ is the minimum number of Monte Carlo steps needed to generate effectively independent, identically distributed samples in the Markov chain. The smallest possible value for $`\tau `$ is 1, which represents independent sample by every step. The usual integrated autocorrelation time  differs from our definition by a factor of 2, $`\tau _{int}=\tau /2`$.

The critical slowing down manifests itself by the fact that $`\tau L^z`$ at the critical temperature $`T_c`$ where a second order phase transition occurs. Here $`L`$ is the linear dimension of the system ($`N=L^d`$ in $`d`$ dimensions). For the local algorithms for many models and in any dimensions, $`z2`$. This suggests bad convergence, specially for large systems. At a first-order phase transition, where some thermodynamic variables change discontinuously, the situation is even worse — $`\tau `$ diverges exponentially with system sizes.

For the two-dimensional Ising model, a phase transition occurs  at $`kT_c/J2.269`$. The magnetization is non-zero below this temperature and becomes zero above $`T_c`$. In addition, in the limit of large system ($`L\mathrm{}`$), heat capacity per spin and fluctuation of magnetization diverge. These intrinsic properties make computer simulation near critical point very difficult.

Cluster algorithms  overcome this difficulty successfully. For example, for the two-dimensional Ising model, the dynamical critical exponent defined in $`\tau L^z`$ is reduced from 2.17 for the single-spin flip  to much small value for the cluster algorithm . It turns out that a precise characterization of the Swendsen-Wang dynamical critical exponent is very difficult, due to weak size dependence of the decorrelation time. Table 1 represents a recent extensive calculation, based on the definition Eq. (13) for the total energy, rather than the usual method of extracting information from time-dependent auto-correlation functions. In this calculation, the variance of the sum of $`M`$ consecutive energies (c.f. Eq. (5)) are computed explicitly. Good convergence to the limiting value is already achieved typically for $`M10^2`$. An $`1/M`$ extrapolation is used to get more accurate estimates of the limit. From this calculation, we find for the two-dimensional Ising model, the convergence to asymptotics are slow. It appears that the divergence is slightly faster than $`\mathrm{log}L`$. If we fit to power law $`L^{z_{sw}}`$, using two consecutive sizes of $`L`$ and $`2L`$, the exponent $`z_{sw}`$ decreases from 0.35 to 0.21. It is not clear whether $`z_{sw}`$ will convergence to a finite value or continue decreasing to 0 in the limit of $`L\mathrm{}`$.

The nonlocal cluster algorithm introduced by Swendsen and Wang for the the Ising model (and more generally the Potts model) goes as follows: (1) go over each nearest neighbor pair and create a bond with probability

$$1\mathrm{exp}\left(2J\delta _{\sigma _i,\sigma _j}/(kT)\right).$$
(14)
That is, if the two nearest neighbor spins are the same, a bond is created between them with probability $`1\mathrm{exp}(2J/kT)`$; if spin values are different, there will be no bond. (2) Identify clusters as a set of sites connected by zero or more bonds (i.e., connected component of a graph). Relabel each cluster with a fresh new value $`+1`$ or $`1`$ at random.

We note that each Monte Carlo step per spin still takes $`O(1)`$ in computational cost. The method is applicable to models containing Ising symmetry, i.e., the energy is the same when $`\sigma _i`$ is changed to $`\sigma _i`$ globally.

The algorithm is based on a mapping of Ising model to a random cluster model of percolation. Specifically, we have

$`Z`$ $`=`$ $`{\displaystyle \underset{\{\sigma \}}{}}\mathrm{exp}\left(K{\displaystyle \underset{i,j}{}}(\delta _{\sigma _i,\sigma _j}1)\right)`$ (15)
$`=`$ $`{\displaystyle \underset{\{\sigma \}}{}}{\displaystyle \underset{\{n\}}{}}{\displaystyle \underset{i,j}{}}\left[p\delta _{\sigma _i,\sigma _j}\delta _{n_{ij},1}+(1p)\delta _{n_{ij},0}\right]`$ (16)
$`=`$ $`{\displaystyle \underset{\{n\}}{}}p^b(1p)^{Ndb}2^{N_\mathrm{c}},`$ (17)
where $`K=2J/kT`$, $`p=1\mathrm{exp}(K)`$, $`\delta _{i,j}`$ is the Kronecker delta, and $`b`$ is the number of bonds, $`d`$ is the dimension of a simple hypercubic lattice, and $`N_\mathrm{c}`$ is the number of clusters. $`\sigma _i`$ is spin on the site and $`n_{ij}=0,1`$ is the bond variable between the sites $`i`$ and $`j`$. It is evident that the moves in the Swendsen-Wang algorithm preserves configuration probability of the augmented model containing both the spins and bonds.

A single cluster variant due to Wolff  is very easy to program. The following C code generates and flip one cluster:

```
      #define Z  4
      #define L  8
      double p = 0.7;
      int s[L*L];

      void flip(int i, int s0)
      {
         int j, nn[Z];

         s[i] = - s0;
         neighbor(i, nn);
         for(j = 0; j < Z; ++j)
            if(s0 == s[nn[j]] && drand48() < p)
               flip(nn[j], s0);
      }

```

In this single cluster version, a site $`i`$ is selected at random. The value of the spin before flip is $`s_0`$. It flips the spin and looks for its neighbors. If the value of the neighbor spins are the same as $`s_0`$, with probability $`p`$ a neighbor site becomes part of the cluster. This is performed recursively. It turns out that single cluster is somewhat more efficient than the original Swendsen-Wang, particularly in high dimensions. For dimensions greater than or equal to 4, Swendsen-Wang dynamics gives the dynamic critical exponent $`z=1`$, while Wolff single cluster is $`z=0`$ .

It is easy to see why the single cluster algorithm of the above works. Let us consider a general cluster flip algorithm with two bond probabilities: $`P_s`$ will be the probability of connecting two parallel spin sites; $`P_d`$ the probability of connecting anti-parallel sites. Consider the transition between two configurations $`A`$ and $`B`$, characterized by flipping a cluster $`C`$. A cluster is growing from a site $`i`$ until the perimeter of the cluster is not connected to the outside. The transition probabilities can be written down as

$`W(AB)=P_I{\displaystyle \underset{C()}{}}(1P_s){\displaystyle \underset{C()}{}}(1P_d),`$ (18)
$`W(BA)=P_I{\displaystyle \underset{C()}{}}(1P_d){\displaystyle \underset{C()}{}}(1P_s).`$ (19)
Note that the bond configuration probabilities are the same in the interior of the clusters. The difference occurs at the boundary $`C`$, where parallel spins ($``$) in $`A`$ becomes anti-parallel spins ($``$) in $`B`$, or vice versa. Detailed balance requires that

$$\frac{W(AB)}{W(BA)}=\frac{P(B)}{P(A)}=(1P_s)^{N_1N_2}(1P_d)^{N_2N_1}=e^{\frac{\mathrm{\Delta }E}{kT}},$$
(20)
where $`N_1`$ is the number of parallel spins on boundary of $`C`$ in configuration $`A`$; and $`N_2`$ is the number of anti-parallel spins on boundary of $`C`$ in configuration $`A`$. Since we have $`\mathrm{\Delta }E=E_BE_A=2(N_1N_2)J`$, we obtain

$$\frac{1P_s}{1P_d}=\mathrm{exp}\left(\frac{2J}{kT}\right).$$
(21)
Although the algorithm is valid for any $`0P_d<1`$, it is most efficient at $`P_d=0`$, the Coniglio-Klein  bond probability value.

A quite large number of statistical models can be treated with cluster algorithms, with varied success. Excellent performance has been obtained for Ising model, Potts models, and antiferromagnetic Potts models , XY model and general $`O(n)`$ models , $`\varphi ^4`$ field-theoretic model , some regularly frustrated models , six-vertex model , etc. Cluster algorithms are proposed for hard sphere fluid systems , quantum systems , microcanonical ensembles , conserved order parameters , etc. Invaded cluster algorithm  and other proposal  are excellent methods for locating critical points. The cluster algorithms are also used in image processing .

The cluster algorithms do not help much in temperature-driven first-order phase transition – the slow convergence has been shown rigorously . Models with frustration, spin glass being the archetype, do not have efficient cluster algorithms, although there are attempts  with limited success. Breakthrough in this area will have a major impact on the simulation methods.

## 4 Reweighting Methods

In this and the following sections, we discuss a class of Monte Carlo simulation approaches that aim at an efficient use of data collected, and sampling methods that enhance rare events.

The computation of free energy, Eq. (2), poses a difficult problem for Monte Carlo method. A traditional method is to use thermodynamic integration, e.g.,

$$\frac{F(T_2)}{T_2}\frac{F(T_1)}{T_1}=_{T_1}^{T_2}\frac{E}{T^2}𝑑T,$$
(22)
based on the relation $`E=\mathrm{ln}Z/\beta `$, where $`\beta =1/(kT)`$.

If we can estimate the density of states (the number of states with a given energy $`E`$ for discrete energy models), then we can compute free energy, as well as thermodynamic averages. The result is obtained as a function of temperature $`T`$, rather than a single datum point for a specific value of $`T`$, as in standard Monte Carlo simulation.

This idea has been pursued over the last decade by Ferrenberg and Swendsen , Berg et al , Lee , Oliveira et al , and Wang . Consider the following decomposition of summations over the states

$$\underset{\{\sigma \}}{}A(\sigma )e^{\left(\frac{E(\sigma )}{kT}\right)}=\underset{E}{}e^{\frac{E}{kT}}\underset{E(\sigma )=E}{}A(\sigma )=\underset{E}{}e^{\frac{E}{kT}}n(E)A(\sigma )_E,$$
(23)
where $`A_E`$ is the microcanonical ensemble average,

$$A_E=\frac{1}{n(E)}\underset{E(\sigma )=E}{}A(\sigma ).$$
(24)
Since the state space is exponentially large ($`2^N`$ for the Ising model with $`N`$ spins), and the range of $`E`$ is typically of order $`N`$, if $`n(E)`$ can be computed accurately, the task is done. The canonical average of $`A`$ is related to the microcanonical average through

$$A_T=\frac{\underset{E}{}A_E\mathrm{exp}\left(\frac{E}{kT}\right)n(E)}{_E\mathrm{exp}\left(\frac{E}{kT}\right)n(E)},$$
(25)
and free energy is computed as

$$F=kT\mathrm{ln}\underset{E}{}\mathrm{exp}\left(\frac{E}{kT}\right)n(E).$$
(26)
### 4.1 Histogram method

Ferrenberg and Swendsen  popularized a method which in a sense is to compute the density of states (up to a multiplicative constant) in a range close to a given simulation temperature. This method is generalized as multiple histogram method to combine simulations at differential temperatures, to get the whole energy range . We discuss here only the single histogram method for its simplicity.

During a normal canonical simulation at fixed temperature $`T^{}`$, we collect the histogram of energy, $`H(E)`$, which is proportional to probability distribution of energy,

$$H(E)=cn(E)\mathrm{exp}\left(E/(kT^{})\right).$$
(27)
The constant $`c`$ is related to the partition function, $`c=M/Z(T^{})`$, where $`M`$ is the total number of samples collected. From the above equation, we find $`n(E)H(E)\mathrm{exp}\left(E/(kT^{})\right)`$. With this information, we can compute the free energy difference between temperature $`T^{}`$ and a nearby temperature $`T`$. Similarly, moments of energy can be computed after the simulation, through histogram reweighting,

$$E^n_T=\frac{\underset{E}{}E^nH(E)\mathrm{exp}\left(E/(kT)+E/(kT^{})\right)}{_EH(E)\mathrm{exp}\left(E/(kT)+E/(kT^{})\right)}.$$
(28)
The range of $`E`$ that the histogram data can be collected at a fixed temperature is limited by the energy distribution, which for the canonical distribution away from critical point, is of order of $`\sqrt{N}`$. The whole range of energy $`E`$ is of order $`N`$. This limit the usefulness of single histogram method.

### 4.2 Multicanonical Monte Carlo

The multicanonical Monte Carlo method has been shown to be very effective to overcome supercritical slowing down, reducing the relaxation time from exponential divergence with respect to system size to a power, at the first-order phase transitions . Multicanonical ensemble flattens out the energy distribution, so that the computation of the density of states $`n(E)`$ can be done for all values of $`E`$. A multicanonical ensemble is defined to following the probability density for the states as

$$p(\sigma )=f\left(E(\sigma )\right)\frac{1}{n\left(E(\sigma )\right)},$$
(29)
such that the energy histogram $`H(E)`$ ($`n(E)f(E)`$) is a constant. From the histogram samples obtained by a simulation with the weight of state at energy $`E`$ as $`f(E)`$, the density of state can be computed  from $`n(E)=H(E)/\left(f(E)c\right)`$. However, unlike canonical simulation where $`f(E)=e^{E/kT}`$ is given, in a multicanonical simulation, $`f(E)`$ is unknown to start with.

Berg proposed an iterative method to compute the weight in a parametrized form $`f(E)=\mathrm{exp}\left(\beta (E)E+\alpha (E)\right)`$, starting with no information, $`f_0(E)=\mathrm{const}`$. A new estimate at iteration $`n`$ is then based on the results of all previous iterations. We refer to references  for details.

## 5 Transition Matrix Monte Carlo and Flat Histogram Algorithm

The flat histogram algorithm offers an efficient bootstrap to realize the multicanonical ensemble, while transition matrix Monte Carlo utilizes more data that can be collected in a simulation to improve statistics.

### 5.1 Transition matrix

We start from the detailed balance equation for some given dynamics:

$$p(\sigma )W(\sigma \sigma ^{})=p(\sigma ^{})W(\sigma ^{}\sigma ).$$
(30)
By summation over the states $`\sigma `$ of fixed energy $`E`$, and $`\sigma ^{}`$ of fixed energy $`E^{}`$, and assuming that the probability of the state is a function of energy only, $`p(\sigma )f\left(E(\sigma )\right)`$, we get

$$n(E)f(E)T(EE^{})=n(E^{})f(E^{})T(E^{}E),$$
(31)
where the transition matrix in the space of energy is defined as

$$T(EE^{})=\frac{1}{n(E)}\underset{E(\sigma )=E}{}\underset{E(\sigma ^{})=E^{}}{}W(\sigma \sigma ^{}).$$
(32)
The matrix $`T`$ has a number of interesting properties: it is a stochastic matrix in the sense of $`T(EE^{})0`$ and $`_E^{}T(EE^{})=1`$; the stationary solution of $`T`$ is the energy distribution $`n(E)f(E)`$; the dynamics associated with $`T`$ is considerably faster than that of $`W`$ .

We specialize to the case of single-spin-flip dynamics for the Ising model. The transition matrix $`W`$ for the spin states consists of a product of two factors, the probability of choosing a spin to flip $`S(\sigma \sigma ^{})`$, and the flip rate $`a(EE^{})=\mathrm{min}(1,f(E^{})/f(E))`$. We have $`S(\sigma \sigma ^{})=0`$ less the two configurations $`\sigma `$ and $`\sigma ^{}`$ differ by one spin, in this case, the value of $`S`$ is $`1/N`$, where $`N`$ is the number of spins in the system. Using these results, we can rewrite the transition matrix as

$$T(EE^{})=\frac{1}{N}N(\sigma ,E^{}E)_Ea(EE^{}),EE^{}.$$
(33)
The diagonal elements are determined by normalization. Substituting Eq.(33) into Eq.(31), using the relation between $`a(EE^{})`$ and $`f(E)`$, we obtain

$$n(E)N(\sigma ,E^{}E)_E=n(E^{})N(\sigma ^{},EE^{})_E^{}.$$
(34)
This is known as broad histogram equation  which forms the basis for the flat histogram algorithm presented below. Additionally, this equation also gives us a way of computing the density of states $`n(E)`$ by the quantity $`N(\sigma ,E^{}E)_E`$ obtained from spin configurations generated from any distribution $`f\left(E(\sigma )\right)`$. The quantity $`N(\sigma ,\mathrm{\Delta }E)`$ is the number of ways that the system goes to a state with energy $`E+\mathrm{\Delta }E`$, by a single spin flip from state $`\sigma `$. The angular brackets indicate a microcanonical average:

$$N(\sigma ,\mathrm{\Delta }E)_E=\frac{1}{n(E)}\underset{E(\sigma )=E}{}N(\sigma ,\mathrm{\Delta }E).$$
(35)
### 5.2 Flat histogram algorithm

The following algorithm  generates flat histogram in energy and realizes the multicanonical ensemble.

1. Pick a site at random.
2. Flip the spin with probability

$$a(EE^{})=\mathrm{min}(1,\frac{N(\sigma ^{},\mathrm{\Delta }E)_{E+\mathrm{\Delta }E}}{N(\sigma ,\mathrm{\Delta }E)_E}),$$
(36)
where the current state $`\sigma `$ has energy $`E`$, the new state $`\sigma ^{}`$ has energy $`E^{}=E+\mathrm{\Delta }E`$.
3. Accumulate statistics for $`N(\sigma ,\mathrm{\Delta }E)_E`$.
4. Go to 1.

We note by virtue of Eq. (34), the flip rate is the same as that in multicanonical simulation with a weight $`1/n(E)`$ and Metropolis acceptance rate. While in multicanonical sampling, the weight is obtained through several simulations iteratively, the quantities $`N(\sigma ,\mathrm{\Delta }E)_E`$ is much easier to obtain, through a single simulation. This quantity serves a dual purpose—it is used to construct a Monte Carlo algorithm (used as input), and at the same time, it is used to compute the density of states (output of the simulation). Clearly, this is circular unless approximation is made. We have considered replacing the exact microcanonical average by an accumulative average, over the history of simulation generated so far, i.e.,

$$N(\sigma ,\mathrm{\Delta }E)_E\frac{1}{H(E)}\underset{i=1}{\overset{M}{}}\delta _{E(\sigma ^i),E}N(\sigma ^i,\mathrm{\Delta }E),$$
(37)
where $`\{\sigma ^i,i=1,2,\mathrm{}\}`$ is the sequence of states generated with the algorithm given above; $`H(E)`$ is the number of samples accumulated at the energy bin $`E`$. In case the data for computing the flip rate is not available, we simply accept the move to explore the new state. A more rigorous way of doing simulation is to iterate the simulation with constant flip rate. For example, after the first simulation, we compute a first estimate to the density of states. In a second simulation, we perform multicanonical simulation á la Lee . The data collected in the second run for $`N(\sigma ,\mathrm{\Delta }E)_E`$ will be unbiased. It is found that even with a single simulation, the results converge to the exact values (for $`N(\sigma ,\mathrm{\Delta }E)_E`$ and $`n(E)`$) for sufficiently long runs, even though a rigorous mathematical proof of the convergence is lacking.

Wang and Landau  proposed recently a new algorithm that works directly with the density of states $`n(E)`$. The simulation proceeds with the flip rate $`\mathrm{min}(1,n(E)/n(E^{}))`$, but the value of the density of states is updated after every move by $`n(E)n(E)f`$ and letting $`f1`$ for convergence. Excellent results were obtained. A careful comparison with flat histogram method is needed.

### 5.3 N-fold way (rejection-free moves)

In Metropolis algorithm, moves are sometimes rejected. This rejection is important for realizing the correct stationary distribution. In 1975 Bortz, Kalos, and Lebowitz  proposed a rejection-free algorithm. It is still based on Metropolis flip rate, but the waiting due to rejection is taking into account by considering all possible moves. The Bortz-Kalos-Lebowitz N-fold way algorithm for the Ising model goes as follows:

1. Compute the acceptance probability for one attempt of a move

$$A=\underset{\mathrm{\Delta }E}{}\frac{N(\sigma ,\mathrm{\Delta }E)}{N}a(EE+\mathrm{\Delta }E).$$
(38)
2. Pick an energy change $`\mathrm{\Delta }E`$ according to probability

$$p_{\mathrm{\Delta }E}=\frac{N(\sigma ,\mathrm{\Delta }E)}{AN}a(EE+\mathrm{\Delta }E).$$
(39)
flip a site belonging to $`\mathrm{\Delta }E`$ with probability 1. The site is choosing from the $`N(\sigma ,\mathrm{\Delta }E)`$ sites with equal probability.
3. One N-fold-way move is equivalent to $`1/A`$ moves in the original dynamics. Thus thermodynamic averages are weighted by $`1/A`$, i.e., $`g=_t(g_t/A)/_t1/A`$, where summation is over every move $`t`$.

In order to implement step 2 efficiently, additional data structure is needed so that picking a spin in a given class characterized by $`\mathrm{\Delta }E`$ is done in $`O(1)`$ in computer time.

### 5.4 Equal-hit algorithm

Combining N-fold way and flat histogram algorithm is easy, since the important quantity $`N(\sigma ,\mathrm{\Delta }E)`$ is already computed in flat histogram algorithm. The flip rate $`a`$ is given by formula (36). In the flat histogram algorithm, the probability that the energy of the system is $`E`$ is a constant, i.e.,

$`H(E)`$ $`=`$ $`\delta _{E(\sigma ),E}n(E)f(E)=\text{const}`$ (40)
$`=`$ $`\delta _{E(\sigma ),E}{\displaystyle \frac{1}{A(\sigma )}}_N/{\displaystyle \frac{1}{A(\sigma )}}_N.`$
The averages in the second line of the above equation refer to samples generated in an N-fold way simulation.

In equal-hit algorithm (ensemble) , we require that the number of “fresh configurations” generated at each energy is a constant. More precisely equal-hit ensemble is defined by

$$u(E)=\delta _{E(\sigma ),E}_N=\text{const}.$$
(41)
One possible choice of the flip rate is

$$a(EE^{})=\mathrm{min}(1,\frac{\frac{1}{A}_{E^{},N}N(\sigma ^{},EE^{})_E^{}}{\frac{1}{A}_{E,N}N(\sigma ,E^{}E)_E}),$$
(42)
where

$$\frac{1}{A}_{E,N}=\frac{\delta _{E(\sigma ),E}\frac{1}{A(\sigma )}_N}{\delta _{E(\sigma ),E}_N}$$
(43)
is the inverse total acceptance rate $`1/A`$ arithmetic averaged over the N-fold way samples at energy $`E`$.

The histogram $`H(E)`$ generated in the equal-hit algorithm depends on the precise dynamics (the rate $`a`$) used. Since there are many possible choices of the rate, such “equal-hit ensemble” is not unique.

### 5.5 Determination of the density of states

While Eq. (34) gives us a way of obtaining the density of states, there are more equations than unknowns. We consider two optimization methods. The first method is based on the transition matrix itself. We define $`T_{E,\mathrm{\Delta }E}=N(\sigma ,\mathrm{\Delta }E)_E/N`$. Symbols with hat being unknown, and $`T_{E,\mathrm{\Delta }E}`$ the Monte Carlo estimate, consider

$$\mathrm{Minimize}\underset{E,\mathrm{\Delta }E}{}\sigma _{E,\mathrm{\Delta }E}^2\left(\widehat{T}_{E,\mathrm{\Delta }E}T_{E,\mathrm{\Delta }E}\right)^2,$$
(44)
subject to $`0\widehat{T}_{E,\mathrm{\Delta }E}1`$, $`_{\mathrm{\Delta }E}\widehat{T}_{E,\mathrm{\Delta }E}=1`$, and $`\widehat{T}_{E,1}\widehat{T}_{E+1,1}\widehat{T}_{E+2,2}=\widehat{T}_{E,2}\widehat{T}_{E+2,1}\widehat{T}_{E+1,1}`$. The last constraint needs more explanation. We assume that the energy level is equally spaced (as in the Ising model). Consider three energy levels, $`E`$, $`E+1`$, $`E+2`$. If we write down three equations of type (34), for transitions from $`E`$ to $`E+1`$, $`E+1`$ to $`E+2`$, and $`E+2`$ to $`E`$, we can cancel the density of states by multiplying the three equations together. This leaves the last equation above, and it is known as TTT identity. It can be shown that multiple T identities (four or more) are not independent, and they need not put in as constraints. For Ising model there is also one additional symmetry constraint, $`T_{E,i}=T_{E,i}`$.

When the solution for $`T`$ is found, we can use any of the energy detailed balance equation to find density of states $`n(E)`$. The TTT identity guarantees that the answer is unique whichever detailed balance equation is used.

The second method is based on optimization directly with variable $`n(E)`$, actually $`S(E)=\mathrm{ln}n(E)`$, by

$$\mathrm{minimize}\underset{E,\mathrm{\Delta }E}{}\sigma _{E,\mathrm{\Delta }E}^2\left(\widehat{S}(E+\mathrm{\Delta }E)\widehat{S}(E)\mathrm{ln}\frac{T_{E,\mathrm{\Delta }E}}{T_{E+\mathrm{\Delta }E,\mathrm{\Delta }E}}\right)^2,$$
(45)
subject to, for the Ising model, $`_En(E)=2^N`$, where $`N`$ is the total number of spins in the system. In addition, we can put in the known fact that the ground states are doubly degenerate, $`n(0)=2`$.

Figure 1 shows one of the simulation results of the density of states, using the second method. The errors by comparison with exact known values  are presented in the insert of the figure. The density of states $`n(E)`$ is determined to an accuracy of better than 2 percents in a matter of few minutes of computer time.

The flat-histogram dynamics is used to study spin glasses . The dynamic characteristics is quite similar to multi-canonical method of Berg. The study of lattice polymer and protein folding is under way. For related ideas and approaches, see refs..

## 6 Conclusion

The phenomenon of critical slowing down can be effectively dealt with by cluster algorithms for a large class of statistical mechanics models. We reported new and very accurate results for the decorrelation times of the Swendsen-Wang dynamics. Large size asymptotic behavior is analyzed. For super-critical slowing down occurring in first-order phase transitions, multicanonical ensemble simulation and flat-histogram or equal-hit algorithms are very effective. Since the latter algorithms and associated transition matrix method is efficient in computing the density of states, this method can also be useful for general counting problems by Monte Carlo method.

## Acknowledgements

The author thanks Prof. R. H. Swendsen for many of the work discussed here. He also thank Z. F. Zhan, T. K. Tay, and L. W. Lee for collaborations. This work is supported by a NUS research grant R151-000-009-112 and Singapore-MIT Alliance.